Thompson - Computing for Scientists and Engineers (523188), страница 32
Текст из файла (страница 32)
The coefficients for calculating its analytical derivativesare given in Table 4.1. In Section 4.4 we explored polynomial approximations forfirst and second derivatives of functions, including the working function. The results appear in Tables 4.4- 4.7. We can now use the spline methods to estimatethe same derivatives, so that we can judge the relative merits of the two techniques.With the same nine-point spline as used for interpolating values, one obtains thefirst derivatives at the knots summarized in Table 5.2 and Figure 5.4.172FITTING CURVES THROUGH DATATABLE 5.2 First derivatives of the working function estimated by a nine-point cubic-spline fitusing either exact endpoints (errors eE ) or natural-spline endpoints (errors eN ) .Notice in Table 5.2 and Figure 5.4 that by using exact endpoint second derivatives one obtains nearly a one order of magnitude improvement in the first derivatives at the knots, and that with both boundary conditions the accuracy of the firstderivatives is much better near the middle of the spline range than near the ends.FIGURE 5.4 Spline derivatives using nine points of the working function, yw.
The exact firstderivative is shown by the solid curve. Scaled errors are shown for the natural-spline derivatives(dashed curve) and for the exact-endpoint spline derivatives (dotted curve).Exercise 5.10Use the program Cubic Splines to make a nine-point cubic-spline fit to thedata in the second column of Table 5.2, which are the working-function valuesat the knots. Run two sets of data, one with exact endpoint derivatives and theother with natural-spline endpoints. Verify that the errors in the first derivativesare as given in Table 5.2 and compare these derivatives graphically with thoseshown in Figure 5.4.
n5.4 INTERPOLATING BY SPLINES173You may be pessimistic about the accuracy of spline fitting, on the basis of thelow accuracy of the spline descriptions of the working function (4.1). This is, however, a pathological function, with derivatives that keep on increasing through thesixth derivative. Its main use is to stress-test the methods and to make inaccuraciesof approximations glaringly obvious.
We now switch to a more harmonious function, the cosine, and explore how well it and its derivatives can be described byspline fitting.Interpolating cosine values and derivativesIn Project 4A (Section 4.5) we explored the numerical derivatives of the cosinefunction at x = 0, using methods based on power series around the point where thederivatives are required. We now repeat the analysis for an extended range of x,namely 0 to /2, and we compare the accuracy of the two methods.Exercise 5.11(a) Consider the cosine function in the range 0 to /2. Use exact endpoint conditions of second derivatives equal to -1 and zero, respectively. Give the exactvalues of the cosine at x = j /16, where j = 0, 1, .
. . . 8, so that there are 9knots. Compare the spline estimates of the first and second derivatives with theanalytical derivatives at the same points, -sin x and -cos x, respectively. Theresults are listed in Table 5.3 and are displayed in Figure 5.5.(b) Compare the spline estimates of derivatives with the results for numericaldifferentiation of the cosine at x = 0 obtained by the methods in Section 4.5.Those estimates of the first and second derivatives are given in Table 4.8. nTable 5.3 Cubic-spline 9-knot fit to the cosine function between x = 0 and x =/2 with exactendpoint conditions. The errors in first derivatives have superscript (1) and in second derivatives thesuperscript is (2).174FITTING CURVES THROUGH DATAFIGURE 5.5 Spline interpolation using nine points of the cosine function, as shown.
Scalederrors in first derivatives (dotted curve) and in second derivatives (dashed curves) are shown for theexact-endpoint spline fit.The standard deviation of the errors in the first derivatives in Table 5.3 is2 x 10-5, which is about the same as the error in the cosine first derivative from thethree-point central-derivative method at x = 0, given in Table 4.8 for a stepsizeh = 0.2, which is about the spacing between the data /16 = 0.196350) given inTable 5.3. Similarly for the errors in the second derivatives, the standard deviationof these in Table 5.3 is 2 x 10-3, comparable to that for the three-point method atx = 0. Overall, a cubic-spline fit of a well-behaved function such as a cosine (arapidly convergent series, as explored in Section 3.2) is three orders of magnitudemore accurate than a fit to our badly behaved working function considered in thepreceding subsection.We have also discovered that derivatives computed using cubic splines can easilybe of comparable accuracy to those obtained from the power-series methods in Sections 4.4 and 4.5.
Although the spline fitting initially requires more number crunching than do the power-series methods, once the spline coefficients have beenfound the computation of derivatives is very speedy.The unstable problem (Section 4.3) of differentiation solved by an unstablemethod involving subtractive cancellation has apparently been bypassed. This meritof splines is more apparent than real because the spline-algorithm formulas (5.15) —(5.18) show that we must take first and second differences of the yj ; therefore,subtractive-cancellation errors are possible. However, such errors are less likelywith splines because the spline procedure does not involve (either explicitly or implicitly) finding a limiting value of the differences as the spacing between the knotsdecreases to zero.5.5INTEGRATION METHODS USING SPLINES1755.5 INTEGRATION METHODS USING SPLINESEven though spline fitting usually emphasizes interpolation properties (as in Section 5.4), splines are also very useful for integration, which is especially straightforward for the cubic spline.
In this section we derive the algorithm for integrationby cubic splines with equal spacing of the knots. We then test the formulas and apply them to the working function and to the cosine. This allows us to compare therelative accuracy of splines with the trapezoid and Simpson formulas in Section 4.6.Deriving the integration algorithmEquation (5.31) shows that the spline fit in thejth interval, s(x), can be expressed as(5.35)where e is given by (5.30). Therefore, in order to integrate the spline fit from xL toxU we need only integrate s(x) over each subinterval, then sum the results. The indefinite integral of sj for any x in the range xj to xj+l is(5.36)In this formula we have expressed it efficiently in the Homer-polynomial algorithmform (Section 4.1), at the expense of loss of clarity in the typography.
An estimateof the total integral from xL to xU that involves no further approximation can therefore be obtained by summing the partial integrals given by (5.36).Note several improvements of cubic spline integration over the trapezoidal orSimpson rules in Section 4.6, namely:Within each interval the spline uses a higher-order polynomial — a cubic ratherthan linear over the interval (trapezoidal rule) or quadratic over double intervals(Simpson rule).Spline fitting provides smooth connections between intervals: there are no suchconstraints for the other methods.Integration with splines uses approximations that are consistent between interpolation and estimation of derivatives.
So, if these operations are also required,splines may be the method of choice.The cost of spline methods compared with the other integration methods is that theinitial fitting algorithm is more time-consuming. This may not be a penalty overall,since for a specified accuracy of the integral the stepsize h can usually be chosenmuch larger in this method than in the others, resulting in fewer yi values beingneeded. Therefore, if the yi are costly to obtain, either as experimental data or ascomputed values, spline methods will often be preferred.
Such methods have beenunder extensive development recently, as described in Chapter 2 of the compendiumon numerical-integration methods by Davis and Rabinowitz, which also gives references to appropriate computer programs.176FITTING CURVES THROUGH DATAThe C function for spline integrationNow that we have an algorithm for integration by using cubic splines, we describethe functions SplineInt,Horner4Poly, and ywInt in program Cubic Splines(Program 5.1). Equation (5.36) gives the integral from the start of the j th segmentto a distance e into the segment.
A direct way to combine the integrals from severalsegments is indicated in Figure 5.6.FIGURE 5.6 Segments and partial segments for the spline integration.The integration segments are combined as follows. In the notation used in Figure 5.6 and in function SplineInt, the integration is to be from xL to xu, and themain program has already checked that these lie within the input range[xmin, xmax] and that xu is not less than xL. The function therefore first computesthe integer knot positions jL and ju and the extra pieces eL and eU, both of which arepositive, because conversion to integers is done by truncation.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
