Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 22
Текст из файла (страница 22)
This was independently discovered by Fritsch and Carlson and published in the more accessible reference [10]. The argument involves studying H’(x)as a function of αn and β n. This is not too complicated since H’(x) is a quadratic on(x n ,x n + l ). A simple condition that guarantees monotonicity is preserved is that α nand β n lie in the interval [0,3]. There are many formulas for α n and β n that satisfy thisrestriction. One given in [9] that works pretty well is to use(3.3 1)with(3.32)for n = 2, 3, .
. . , N - 1. Ifthen the slopes change sign at xn. In sucha case we probably should not impose any requirements on the slope of H(x) at xn.Some people suggest setting bn = 0 when this happens. Others say to go ahead anduse (3.31) as long as there is no division by zero. The heuristic actually chosen canhave a significant impact on the behavior of the shape-preserving cubic spline nearthose regions whereAt the ends the simplest rule is to use b1 =andbN =A better choice is to use the end slope of the quadratic interpolatingthe three closest data points (assuming it satisfies the constraint on a and β); otherpossibilities are given in [9].
With (3.31) and the simple choice for b1 and bN it is easyto show that the sufficient conditions on αn and β n are satisfied. At the ends α1 = 1and β N-1 = l, which are certainly in [0,3]. For n = 2, 3, . . . , N - 1,so αn =andas desired.The algorithm for H(x) is very simple. Compute b1 by whatever formula youotherwise compute bn fromchoose; for n = 2, 3, .
. . , N - 1 take bn = 0 if(3.31), (3.32). Compute bN. The values cn and dn can be computed from (3.30) forn = l, . . . , N - 1 either at the same time the bn are computed or in a separate passover the data. Later in this chapter we provide a routine SVALUE/Spline_value for theevaluation of H(x).Examination of the algorithm shows that on the subinterval [x n ,x n + 1 ], the splinedepends on the data (xn-1, fn-1 ) , (xn, fn ) (xn+l, fn+l ) and (xn+2,fn+2).
It shouldbe no surprise that it depends on data from adjacent subintervals because the first106CHAPTER 3 INTERPOLATIONFigure 3.8 Data from Exercise 3.30 interpolated by the polynomial P8.derivatives of the polynomials in adjacent subintervals have to match at the nodes.Although not as local as the cubic splines of the preceding subsection, the constructionof the shape-preserving spline on a subinterval requires only data from the subintervalitself and the two adjacent subintervals.
As might be expected this H(x) works verywell on data that are always monotone but is less successful on oscillatory data. SeeHuynh [14] for some alternatives.This spline is not very accurate as h = max(xn+l - xn ) tends to zero, but that isnot its purpose. It should be used when the data are “sparse” and qualitative propertiesof the data are to be reproduced. It is a simple and effective automatic French curve.Example 3.8. Exercise 3.30 describes a situation in which a shape-preserving splineis particularly appropriate for the approximation of a function f(C). There are onlyeight data points and it is necessary to approximate the derivative f´(C). The concentration C and diffusion coefficient D(C) that make up the function f(C) = CD(C) tobe approximated are nonnegative.
As can be seen in Figure 3.8, a polynomial interpolant to f(C) is unsatisfactory because it fails to reproduce this fundamental property.Also, it does not reproduce the monotonicity of the data, casting doubt on its use forapproximating f´(C). The shape-preserving spline requires the monotonicity of theinterpolant to match that of the data; Figure 3.9 shows that the result is a much morenplausible fit.CONTINUOUS SECOND DERIVATIVEThe last spline considered has a historical origin. To draw smooth curves through datapoints, drafters once used thin flexible strips of plastic or wood called splines. The3.5 SPLINE INTERPOLATION107Figure 3.9 Data from Exercise 3.30 interpolated by the shape-preserving spline.data were plotted on graph paper and a spline was held on the paper with weights sothat it went over the data points.
The weights were constructed so that the spline wasfree to slip. The flexible spline straightened out as much as it could subject to theconstraint that it pass over the data points. The drafter then traced along the splineto get the interpolating curve. The smooth cubic spline presented here is the solutionof a linearized model of the physical spline. The physical analogy already points outsomething very different about this spline—its value at any point depends on all thedata.To construct the smooth cubic spline, we write once again(3.33)on each [x n ,x n+1 ], 1 < n < N - 1. There are 4(N - 1) free parameters to be determined.The interpolation conditions require that for 1 < n < N - 1(3.34)giving 2N - 2 conditions.
There remain 2N - 2 degrees of freedom that can be usedto make S(x) smooth on all of [x 1 ,x N ]. Note that (3.34) automatically ensures that S iscontinuous on [x l ,x N ]. For S´ to be continuous at interior knots,(3.35)This provides N - 2 conditions, so N degrees of freedom remain. For S´´ to be continuous at interior knots,(3.36)for another N - 2 conditions. Exactly 2 degrees of freedom are left.
This is not enoughto achieve a continuous S”’ (this is undesirable anyway since the resulting S wouldbe a cubic polynomial rather than a piecewise cubic polynomial). There are manypossibilities for the two additional constraints.108CHAPTER 3INTERPOLATIONFigure 3.10 Graphs of S(x) and S´´(X) from Example 3.9.Type 1. S´(x1 ) = f´(x1 ),S´(xN) = f´(xN).Type 2.
S´´(x1) = S´´(xN) = 0.Type 3. S´´(x1) = f´´´(xl),S´´(xN) = f´´´(xN).Type 4. S´´(x1) = S´(xN),S´´(x1) = S´´(xN).For obvious reasons these conditions are known as end conditions. The second condition is the one leading to a spline that approximates the physical spline. The physicalspline straightens out as much as possible past the last data point on each end, so it becomes a straight line with zero second derivative. In the form stated, the first and thirdconditions are useful only if extra information is available about f.
However, the exactslopes or curvatures needed here are often replaced by polynomial approximations inpractice. The last end condition is appropriate when f is periodic with period xN-x1because then f(x) and all its derivatives have the same values at x1 and xN.Example 3.9. LetIt is easily verified that S is in C2[0,2], and satisfies the interpolation conditions S(0) =2, S(1) = 1, S(2) = 4 and the end conditions S´(0) = 1, S´(2) = 13. Graphs of S andS´´ are shown in Figure 3.10. Note that the graph of S appears very smooth, while thatnfor S´´ has an obvious corner at the knot x = 1.Returning to the earlier characterization of S(x), we had 4N - 4 conditions on the4 N - 4 unknowns given by (3.33).
A matrix method is in order, but some preliminarymanipulations will simplify the task considerably. On each interval [x n ,xn +1 ](3.37)3.5 SPLINE INTERPOLATION109(3.38)The interpolation conditions immediately yield, from (3.33),an = fn, 1 < n < N - 1,and also fn+1 = an + bnhn + cnh2n+ dnh3n(3.39)which can be rewritten asbn = (fn+l - fn)/hn - cnhn-dnh2n, 1 < n < N -1.(3.40)This eliminates half of the unknowns. The continuity condition (3.36) on S´´ says that2cn = 2cn-l + 6d n-1 h n-1 or [with cN = S´´(xN)/2](3.41)Only formulas for c1, .
. . , cN remain. They are provided by the two end conditionsand the global continuity of S´. From (3.35) and (3.37) it follows that bn = bn-1 +2cn-1hn-l + 3dn-1h2n-1 for 2 < n < N - 1. Substitution in (3.40) and (3.41) givesfor 2 < n < N - 1 and a rearrangement yields(3.42)Only the first type of end conditions (prescribed slopes) is taken up here. From (3.33),(3.40), and (3.41),so2h l c 1 + h 1 c 2 = 3 {(f 2 - f 1 )/h1 - f´(x1)}.(3.43)Similarly, f´(xN) = S´(xN) leads to(344)Equations (3.41)–(3.44) provide a set of N equations in the N unknowns c1, c2, .
. ., cN.The coefficient matrix has the very special structure110CHAPTER 3INTERPOLATIONSuch a matrix is called tridiagonal (all the nonzero entries lie in three diagonal bands),symmetric (the entry in row i of column j equals the one in row j of column i ), anddiagonally dominant (in each column the magnitude of the entry on the main diagonalexceeds the sum of the magnitudes of the other entries in the column). We saw inSection 2.5.2 that for such matrices the system has a unique solution for each righthand side and the solution can be found accurately using Gaussian elimination withoutany row interchanges.Theorem 3.4.Given the knots x1 < x2 < ··· < xN and fn = f(xn), 1 < n < N,there exists one and only one function S(x) that satisfies each of the following:1.2.3.4.S(x) is a cubic polynomial in each [xn,xn+1], 1 < n < N - 1.S(x) is in C2[x 1 ,x N ].S(xn) = fn, 1 < n < N.S´(x1 ) = f´(x1 ),S´(xN) = f´(xN).For this choice of end conditions, S(x) is called the complete cubic spline.
Thecoefficient matrix has the same structure for end conditions of types (2) and (3) andsimilar results are true for them. With the choice of type (2), S(x) is called the naturalcubic spline. The matrix has a somewhat different form for the periodic end conditions,type (4), but similar results are true and the spline can be computed conveniently.Because the smooth cubic spline depends on all the data, a system of linear equations must be solved to construct it. Often large data sets are to be fit, and if thesolution of the linear equations were anything like as expensive as for a general system, the approach would be impractical.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
