Thompson - Computing for Scientists and Engineers (523188), страница 31
Текст из файла (страница 31)
Detailed discussions ofthe interpolation function, SplineInterp, and of the integration function, SplineInt, are deferred to the sections in which they are used, namely Sections 5.4 and5.5, respectively.The main program, Cubic SplinesThe structure of Cubic Splines is as follows. The control loop in the main program is controlled by n, the number of points in the fit.
Since at least 4 points areneeded to describe a cubic, we input n < 1 to signal program termination.5.3PROGRAM FOR SPLINE FITTING167If n is greater than the number of points assigned for the arrays (starting at 1 tofollow Fortran usage) or if 0 < n < 4, then a new value of n is requested. Otherwise, it’s time to input the n data values to be splined. Then the second-derivativesat the endpoints, sd2 [1] and sd2 [n] , are input. If you want natural splines, justgive both of these derivatives as zero. The final input quantities are the minimumand maximum x values, Since the program already knows the number of steps andthey are uniform, the stepsize (h) is readily computed from xmin and xmax.Exercise 5.5(a) Code and compile program Cubic Splines. At this stage you may wish tobypass coding or use of the integration function SplineInt.
If so, replace it bya “stub,” which will either just print a message that the function was entered orwill return a value of zero, or both.(b) Run the program in order to test a variety of combinations of the controlloops, such as the allowed and terminating values of n, the allowed range of xfor interpolation (which is, by definition, in the interval [ xmin , xmax] ), and(for the same reason) the range allowed for integration.
In addition, in order tosimplify programming, the upper limit on an integral, xU, is not allowed to beless than the lower limit, XL, even though this is mathematically permissible andcould be programmed, at some increase in coding complexity. nThe function SplineFitThe function SplineFit, which does the main work, is organized very much as thespline algorithm subsection in Section 5.1 describes. The only aspect of coding thealgorithm that may be unusual to Fortran programmers is in computing the secondderivatives of the spline.
As (5.20) and (5.21) show, at this stage in the algorithmone makes a downward recurrence on j . Since the C language allows for loopsto decrement, it is most natural to use the decrementing operation (coded as j - -).The required results can also be achieved, with great loss of clarity and naturalness,by using an incrementing loop-control variable (a DO loop in Fortran), then computing the decrementing array index from this variable. It’s messy and error-prone, believe me.We now describe the control structures for using the function Spline Interp.In Section 5.4 we discuss how to program this function. Because splines provide ameans of interpolating, it is reasonable to control the interpolation loop by terminating its execution if the x value at which the spline is to be interpolated is beyond theinput-data range [ xmin , xmax] , which is how the while loop over x in themain program is controlled.
At any chosen x within this range the functionSplineInterp returns the interpolated value, and also the interpolated first, second, and third derivatives. They are then printed by the main program before it requests input of another value of x.168FITTING CURVES THROUGH DATAExercise 5.6Modify the Cubic Splines program, Program 5.1, to be convenient foryour computing environment. In particular, the simple but general-purpose input given there is not convenient for many interactive systems. A practical alternative is to input the data from files or from a spreadsheet.
Similarly, modify theoutput so that it is convenient for your system and graphics software. nNow that we have an algorithm and an understanding of the main parts of Cubic Splines, we are ready to learn how to use it to interpolate, differentiate, andintegrate — techniques that we develop and practice in the next three sections.5.4 INTERPOLATING BY SPLINESSpline fitting is a very useful technique, because it allows us to make a fit to the datawhich is locally smooth and that has decreasing influence from points further away.It is therefore particularly suitable for interpolation between data before using othertechniques for data massaging and analysis.Interpolating values and derivativesTo estimate values at a point x in the range x1 to xn, locate the appropriate interpolation index i by(5.29)Next, find the remainder,(5.30)The interpolated value is then simply obtained from(5.31)The first derivative from the spline interpolation is obtained from (5.31) (noting thatthe derivative with respect to x is the same as the derivative with respect to e)(5.32)The second derivative estimate is thence obtained from (5.32) as(5.33)5.4 INTERPOLATING BY SPLINES169Finally (for a cubic), the third derivative estimate is merely a constant in each interval, so that(5.34)and it is usually discontinuous across intervals.Exercise 5.7(a) Verify the correctness of each of the derivative formulas (5.32) — (5.34)starting with the spline equation (5.31).(b) Suppose that y (x) is exactly a cubic in the range of spline fitting.
Verifythat the spline derivatives and the derivatives of y are the same, at least analytically. nThe C function SplineInterpThe function SplineInterp is coded to compute formulas (5.29) through (5.34).Both the formulas and the code are written in the Horner-polynomial form (Section 4.1) for efficiency. This function also illustrates the use in C of dereferencing,so that the values of the interpolated derivative estimates, * sd1 , * sd2 , * sd3, arepassed back through the function to the program using them. Here are some suggestions for testing the interpolation function.Exercise 5.8(a) Patch in some code in function yw just above the return statement.
Thiscode should look like that immediately above, but you should use formulas for afunction that is a simple cubic or lower-order polynomial whose values and derivatives you can check easily and accurately by hand calculations. Then compile and run this version of Cubic Spline, including the interpolation option.The interpolated and input values should agree within computer roundoff errors.Explain why.(b) A nice way to see the input values at desired x values is to give a list of zeros as input to be fitted. Then the error in the interpolated fit at the knots willjust be the values of your cubic (or lower-order) polynomial. Clever, eh?(c) If your test runs fail to give agreement with a cubic, drop off the cubic termcontributions from the value and derivative formulas, then try again. If discrepancies beyond computer roundoff error persist, keep decreasing the order of thetest polynomial.
The coding bug will be in the derivative one order higher thanthe lowest successful order. Note carefully that you must use exact endpointconditions for each order polynomial. nOnce the interpolation sections of Cubic Spline are running correctly, it is time totest them on our working function presented in Section 4.1.170FITTING CURVES THROUGH DATAInterpolating working-function values and derivativesAs our first test of spline interpolation, consider the working function (4.1), whichwas also considered extensively in Sections 4.1 and 4.4. Program 5.1, CubicSplines, contains the function yw that evaluates the working function and its firstthree derivatives and passes them to the main program.For testing interpolation by cubic splines, it is clear from Figures 5.1 and 5.2,which used only five knots, that reasonable accuracy of values and derivatives forsuch a wiggly (sixth-order) polynomial requires considerably more points, even ifexact endpoint second derivatives (which are usually not known) are used, as Figure 5.2 shows.
The choice of n = 9 knots is a compromise between ease of inputand accuracy of interpolated values and some derivatives of the working function.Theknotsarechosenat x = -0.4 bystepsof h = 0.1 up to x = 0.4. We interpolate values to points midway between the knots. Either exact (subscript E onthe error) or natural-spline (Section 5.2, subscript N on the error) boundary conditions are used.Table 5.1 shows the results numerically, and they are displayed graphically inFigure 5.3.TABLE 5.1 Nine-point cubic spline input yw, errors in interpolated values using exact endpointsecond derivatives (eE) and natural-spline endpoints (eN).Notice from Table 5.1 and Figure 5.3 that the errors in fitting the function values are largest near the endpoints, where there is least information about the function.
Thus, whenever reliable second derivatives are available near the endpointsthey should be used. For our pathological working function and small number ofspline knots, using exact endpoint conditions (eE in Figure 5.3) results in almostone order of magnitude better fit to the function than obtained using natural endpoints (eN).5.4 INTERPOLATING BY SPLINES171FIGURE 5.3 Spline interpolation using nine points of the working function, yw.
Scaled errorsare shown for the natural-spline fit (dashed curve, eN ) and for the exact-endpoint spline fit (dottedcurve, eE ).Exercise 5.9Run the program Cubic Splines using the values of the working function atthe nine knots, as given in Table 5.1. First make a run with exact second derivatives at the endpoints, then make a second run with zero endpoint derivativesfor the natural splines. Second derivatives can be obtained as suggested in Exercise 5.8 (b).
Verify the two sets of numerical results for the errors, as given inTable 5.1. nNow that we have some idea how spline fitting works for function values, let uslearn how to estimate function derivatives by splines. Spline fitting provides a particularly simple, convenient, and accurate way of estimating derivatives of functionsor of data, because the method is derived from considerations of matching derivatives at the knots, as discussed in Section 5.1. Notice that the interpolated third derivative will jump between adjacent knots, often quite quickly, as for our workingfunction, since it has nonzero derivatives up through sixth order.As our first, and very demanding, test of spline derivatives we use the workingfunction,(4.1) or (4.2), a sixth-order polynomial whose properties are discussed extensively in Section 4.1.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
