Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 23
Текст из файла (страница 23)
Fortunately, the system is very special andit is practical to interpolate data involving thousands of nodes. First the tridiagonalsystem (3.42)–(3.44) must be solved. Since it is not necessary to do row interchangesin this case, the elimination formulas are very simple. For reinforcement, let us workthrough the details.To eliminate the first entry in row 2, multiply row 1 byand subtract.
The remaining equations have the same pattern, so at the kth stage multiply row k by thecurrentand subtract from row k + 1. The algorithm for elimination and modification of the right-hand side isfor k = 2, 3, . . . , N beginend k.3.5 SPLINE INTERPOLATION111Back substitution is also easy:The whole computation costs only 3N – 3 multiplications and 2N – 1 divisions. Oncethe c vector is known, vector d can be computed from (3.41), and vector b from (3.40).The storage required is a small multiple of N rather than the N* needed for a generalsystem of equations.We finish this section by discussing some of the mathematical properties of thecomplete cubic interpolator spline S(x).
The physical spline used by drafters can bemodeled using the theory of thin beams. In general, the curvature K(X) of a functionf(x) isand in this theory the expression is linearized to K(x) |f"(x)|. When (S'(x))2 << 1,the quantitycan be regarded as a measure of the curvature of the splineS(x). We prove now that in this measure, any smooth interpolating function satisfyingthe type (1) end conditions must have a curvature at least as large as that of S(x).
Thisis sometimes referred to as the minimum curvature property of the complete cubicspline. The same result is true for the natural cubic spline when the requirement thatthe interpolant satisfy the type (1) end conditions is dropped.Theorem 3.5.If g is any C2[x 1 ,x N ] function that interpolates f over {XI, . . . .x N} and satisfies the type (1) end conditions, thenwhere S(x) is the complete cubic interpolator spline. The same inequality holds forg that do not necessarily satisfy the type (1) end conditions when S(x) is the naturalcubic interpolator spline.Proof First observe thatIf it can be shown that the second integral on the right is zero, thensince the integral of a nonnegative function is always nonnegative, and we are finished.To establish that the desired integral is zero, note that112CHAPTER 3INTERPOLATIONand two integrations by parts giveSince S is a cubic on each [xn,xn+l ], it follows that S(4) 0, so the last integral is zero.Also, (g – S)= (fn+l - fn+l ) - (fn - fn) = 0 since both g and S interpolate f.Thus,which telescopes to (g' – S')S"|xn – (g' – S')S"|x1, and the type (1) end conditionsforce these terms to vanish.
The terms vanish without assuming that g satisfies type(1) end conditions when S(x) is the natural cubic spline because it satisfies the endnconditions S" (xl) = S" (xn) = 0 .While the minimum curvature property is nearly always desirable, there are circumstances in which it is a disadvantage. Note that f certainly interpolates itself, soTheorem 3.5 impliesIn examples where f has very largecurvature, there can be a considerable discrepancy between S and f unless there areenough knots (data points) in the region of large curvature that S can turn sufficientlyfast. Several illustrative examples are given in the next section.Convergence rates analogous to (3.17)-(3.20) for the Hermite cubic spline can beestablished in the complete cubic case. However, proofs are more difficult becauseS(x) is determined by all the data and it is not possible to treat each subinterval independently.
The following result is from [1 1].If f is inand S(x) is the complete cubic interpolatoryTheorem 3.6.spline for f with knots {xl < . . .< xN}, then for any x in [x1 ,xN]where M4 = maxIn contrast to polynomial interpolation, S(x) does converge to f(x) as Naslong as h0. The first and second derivatives of the spline also converge to thecorresponding derivatives of f. Because of this, the spline inherits the shape of fwhen h is small. For example, at a point t where f'(t) > 0, convergence implies thatfor all sufficiently small h, S’(t) >0.
Accordingly, the smooth spline inherits the3.5 SPLINE INTERPOLATION113monotonicity of f for small h, except possibly near the extrema of f. The shapepreserving spline is required only when the data are so sparse that we must imposedirectly the property of monotonicity on the spline. The same argument shows that forsmall h, the smooth spline is convex where f is, except possibly near inflection points.The complete cubic spline will converge when f has fewer than four continuousderivatives on [a,b], just not as fast. Experimentation with a physical spline showsthat the farther a node xk is from a given point t, the less the effect of the value of fkon S(t). This is also true of the mathematical spline, and a careful analysis, see [16],of convergence reveals that the rate of convergence at t depends only on how smoothf(x) is near t. In particular, the convergence rates of the theorem hold on subintervalsof [a,b] where f has four continuous derivatives.In practice, it is usually impossible to use the conclusions of Theorem 3.6 to estimate errors, given only discrete data, since M4 is not available.
As was suggested forpolynomial interpolation, it is wise to reserve some data as a check on the approximation. A graph of S can help in making judgments about the quality of the fit.ROUTINES FOR CUBIC SPLINE INTERPOLATIONTwo routines are provided for the calculation of the complete cubic interpolatory splineS. One, SPCOEF in FORTRAN, Spline_coeff in C, sets up the tridiagonal system(3.42)–(3.44) for {ci }, 1 solves it, and computes {di } and {bi } from (3.41) and (3.40).This routine should be called only once for a given set of data. The coefficients outputfrom this routine are then used in the evaluation routine, SVALUE in FORTRAN,Spline-value in C. It is called once for each point t where S(t) is to be evaluated.
Aroutine to compute the coefficients defining the shape-preserving interpolant is quiteuseful. It can be written easily by modifying SPCOEF or Spline_coeff so as to use theformulas of Section 3.5.2. Proceeding in this way, SVALUE or Spline_value can beused for the evaluation of both kinds of spline.Instead of using the slopes f´(x1) and f´(xN) needed for the end conditions of thecomplete cubic spline, the routines provided interpolate the four data points nearesteach end with cubits, and the slopes of these approximations are used in (3.43) and(3.44).
As h0, the resulting spline converges to the complete cubic spline. In practice this approximation works well enough if N is not too small. The approximationis not plausible for Example 3.8 because there are only eight data points, all of whichare used to approximate the derivatives at the end as well as the function throughout the interval.
When the data are this sparse, the shape-preserving spline is moreappropriate.A typical call in FORTRAN isCALL SPCOEF (N, X, F, B, C, D, FLAG)andflag = Spline_coeff(n, x, f, b, c, d);in the C and C++ versions. The input vectors X and F hold the data points (xi ,fi )to be interpolated and N is the number of such points. The output vectors B, C, and114CHAPTER 3INTERPOLATIOND contain the coefficients of the cubits. In normal circumstances the output variableFLAG is set to zero. However, if the input N < 2, then no calculations are performedand FLAG := - 1.
If the entries of X are not correctly ordered (so that some hj < 0),then FLAG := -2.To evaluate the spline the FORTRAN version SVALUE first finds an index i suchthat xi < t < xi+1 and then the ith cubic is evaluated to get S(t). A typical call in FORTRAN isCALL SVALUE (N, X, F, B, C, D, T, INTERV, S, FLAG)andflag = Spline_value(n, x, f, b, c, d, t, interval, s);in the C++ version. The last two parameters are output, so their addresses must explicitly be passed in C:flag = Spline-value(n, x, f, b, c, d, t, &interval, &s);As usual, arrays in the C and C++ versions are indexed starting at 0 rather than 1 asis typical of FORTRAN. The parameters N, X, F, B, C, and D have the same meaningfor SPCOEF and Spline_coeff.
The last three are input quantities that must have beenset by a prior call to SPCOEF or Spline_coeff. The variable T holds the point wherethe evaluation is to be made and the answer comes back in S. If the index i satisfyingxi < T < xi+1 is known, this can be input using the variable INTERV or interval, asthe case may be. However, it is not necessary to do this since the code will calculatethe correct value and assign it to INTERV or interval.
The normal value of FLAG (thereturn value in the C version) is zero. When N < 2, FLAG is returned with the value- 1. If T < x1, then FLAG is set to 1, and the cubic for [x 1 ,x 2 ] is used for S. If T > xN,then FLAG is set to 2, and the cubic for [xN-1, xN] is used for S.Example 3.1O. A sample driver is provided to interpolate sinx over { 0, 0.2, 0.4, 0.6,0.8 } ; the resulting S(x) is then tabulated at { 0.1, 0.3, 0.5, 0.7, 0.9 } to yield thefollowing.Note that only one call is made to SPCOEF or Spline_coeff even though the spline isnevaluated at five points. Why is FLAG = 2 in the last line?Example 3.11. A graph of the spline S(v) interpolating the 16 indicated data pointsfrom Example 3.6 appears in Figure 3.11.
It is a dramatically better approximationthan the polynomial of high degree appearing in Figure 3.2. Values of S at some ofthe reserved data points are S(0.29) = -2.88, S(1.11) = -1.57, S(5) = 5.50, S(8.5) =3.5 SPLINE INTERPOLATION115Figure 3.11 S(v) from Example 3.11.3.5 SPLINE INTERPOLATION111Back substitution is also easy:CN := YN/~Nfork=N-l,N-2,...,1 beginck := (Yk-Pk*Ck+l)/akend k.The whole computation costs only 3N - 3 multiplications and 2N - 1 divisions. Oncethe c vector is known, vector d can be computed from (3.41), and vector b from (3.40).The storage required is a small multiple of N rather than the N2 needed for a generalsystem of equations.We finish this section by discussing some of the mathematical properties of thecomplete cubic interpolatory spline S(x). The physical spline used by drafters can bemodeled using the theory of thin beams.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
