Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 50
Текст из файла (страница 50)
. . . XI(N+1) STRICTLY INCREASING SEQUENCE OF BREAKPOINTS.CC(l,I), C(2,I), VALUE AND FIRST DERIVATIVE AT XI (I), I=1 ,... ,N+l,COF THE PIECEWISE CUBIC FUNCTION.C****** O U T P U T ******C C(l,I), C(2,1), C(3,I), C(4,I) POLYNOMIAL COEFFICIENTS OF THE FUNCCTION ON THE INTERVAL (XI (I), XI(I+1)) , I=l,...,N .CDO 10 I=l,NDX = XI(I+1) - XI (I)DIVDFl = (C(l,I+l) - C(l,I))/DXDIVDF3 = C(2,I) + C(2,I+l) - 2.*DIVCFlC(3,I) = (DIVDFl - C(2,I) - DIVDF3)/DX10Ct4,I) = DIVDF3/ (DX*DX)RETURNENDExample 6.15 Solve the interpolation problem of Example 2.4 using piecewise-cubicHermite interpolation; i.e., for N = 2, 4, . .
. , 16, chooseand interpolatef(x) = (1 + x 2 ) - 1at these points, estimating as before the maximum interpolation error in [-5, 5].The following FORTRAN program solves this problem:C PROGRAM FOR EXAMPLE 6.15 .INTEGER I,J,K,NREAL C(4,l7) ,ERRMAX,H,X(l7) ,YCPIECEWISE CUBIC HERMTTE INTERPOLATION AT EQUALLY SPACED POINTSTO THE FUNCTIONCF(Y) = l./(l. + Y*Y)CPRINT 660600 FORMAT('1 N',5X,'MAXIMUM ERROR')DO 40 N=2,16,2H = l0./FLOAT(N)DO 10 I=l,N+IX(I) = FLOAT(I-l)*H - 5.C(l,I) = F(X(I))CC(2,I) = F'(X(1))10C(2,I) = -2.*X(I)*C(1,I)**2CALL CALCCF ( X, C, N )ESTIMATE MAXIMUM INTERPOLATION ERROR ON (-5,5).CERRMAX = 0.DO 30 I=1,101Y =.1*I - 5.ERRMAX = MAX(ERRMAX, ABS(F(Y)-PCUBIC(Y,X,C,N)))30CONTINUE40PRINT 640, N,ERRMAX640 FORMAT(I5,E18.7)STOPEND288APPROXIMATIONCOMPUTER OUTPUT FOR EXAMPLE 6.15N246810121416MAXIMUM ERROR4.9188219E - 012.1947326E - 019.1281965E - 023.512825OE - 021.2705882E - 024.0849234E - 031.6011164E - 031.6953134E - 03In contrast to polynomial interpolation (see Example 2.4) the maximum error nowdecreases quite nicely as N increases.The error in piecewise-cubic Hermite interpolation is easily estimated.Since, forwhere Pi (x) interpolates f(x) at xi ,xi , xi+1, xi+1, it follows from (2.37) that, forf(x) - g 3 (x) = f[x,i xi , xi+l, xi+l, x](x - xi ) 2 (x - xi+1 ) 2provided f(x) is four times continuously differentiable.
Further,ThereforeFor a < x < b:(6.78)Piecewise-cubic Hermite interpolation requires knowledge of f´(x). Inpractice, it is often difficult, if not impossible, to acquire the needednumbers f´(xi ), i = 1, . . . , N + 1. In such a case, one uses for si somereasonable approximation to f´(xi ), i = 1, . . . , N + 1. Thus, in piecewisecubic Bessel interpolation, one uses(6.79)instead of si = f´(xi ), but proceeds otherwise as before, determining thecoefficients cj,i for the cubic pieces by (6.77). Note that (6.79) requires thetwo additional points x 0 , x N+2 to give some number for the boundaryderivatives s1, sn+1 , of g3(x).
One chooses these points somehow, e.g.,x0 = x3xN+2 = xN-l6.7PIECEWISE-POLYNOMIAL APPROXIMATION289Or, corresponding to the choice x0 = a, xN+2 = b, one uses(6.80)s N+1 = f´(b)if these numbers are available. Yet another possibility is to choose s1 andsN+1 in such a way that g3(x) satisfies the “free-end” conditionss l = f´(a)(6.81)g´´3(a) = g´´3(b) = 0If we continue to use fi = f(xi ), i + 1, . . .
, N + 1, in (6.77), thenregardless of the particular choice of numbers si , i = 1, . . . , N + 1, theresulting piecewise-cubic function g3(x) interpolates f(x) at x1, . . . , xN+1.Further, g 3 (x) is not only continuous, but also continuously differentiableon [a,b], since (6.77) implies thati = 2, . . . , NP´i-1 (xi ) = si = P´i (xi )As we now show, it is always possible to determine the numberss1, . . .
, sN+1 in such a way that the resulting g3(x) is even twice continuously differentiable. This method of determining g3(x) is known as cubicspline interpolation. The name “spline” has been given to the interpolantg 3 (x) in this case, since its graph approximates the position which adraftman’s spline (i.e., a thin flexible rod) would occupy if it were constrained to pass through the points {xi ,fi }, i = 1, . . . , N + 1.The requirement that g 3 (x) be twice continuously differentiable isequivalent to the condition thati=2,...,NP´´i-1(xi ) = P´´i (xi )or with (6.73),2 c3,i-1 + 6c4,i-1∆xi-1 = 2c 3 , ii = 2,. .
.,NHence, with (6.77) we wanti=2, . . ., NIf we use (6.77) to express c4,i-1 and c4,i in terms of the fj 's and sj's, andsimplify, we get(6.82)This is a system of N - 1 linear equations in the N + 1 unknownss1, . . . , sN+1. If we somehow choose s1 and sN+1, for example, by (6.79) or(6.80), we can solve (6.82) for s2, . . . , sN by Gauss elimination (see Chap.4). The coefficient matrix of (6.82) is then strictly row diagonally dominant, hence (see Exercise 4.6-3) invertible, so that (6.82) has then a uniquesolution.
Once we obtain the solution s2 , . . . , sN of the linear system290APPROXIMATION(6.82), we use it, together with the boundary slopes s 1 and s N+1 , inCALCCF to construct the local polynomial coefficients of the interpolating cubic spline.It can be shown (see, e.g., de Boor [40; V(6)]) that the error in thecubic spline interpolant satisfiesFor a < x < b:(6.83)This error bound is only 5 times as big as the error bound (6.78) for cubicHermite interpolation, even though cubic Hermite interpolation uses twiceas much information about the function f(x), viz., the values f´(x i ),i = 2, .
. . , N in addition to the function values. This suggests that theslopes g´3(xi ) of the interpolating spline must be good approximations to thecorresponding slopes f´(xi ) of f(x). One can show (see, e.g., de Boor [40;V(11)-(12))) thatFor a < x < b:(6.84)while, in case of a uniform point sequence, xi = x0 + ih, all i, one even hasFor i = 2, . . . , N:(6.85)This has made cubic-spline interpolation popular as a means for numericaldifferentiation (see Chap. 7).The FORTRAN subprogram SPLINE below uses Gauss eliminationadapted to take advantage of the tridiagonal character of the coefficientmatrix of (6.82) (see Algorithm 4.3) to calculate c2,i = si , i = 2, .
. . , N, asthe solution of (6.82), given the numbers c1,i = fi , i = 1, . . . , N + 1, andc 2,1 = s 1 , c 2,N+l = s N+1 .SUBROUTINE SPLINE ( XI, C, N )PARAMETER NPlMAX=50INTEGER N, MD(NPlMAX),DIAG(NPlMAX),GREAL C(4,N+l) ,XI(N+l),c****** I N P U T ******C XI(l), . . . . XI(N+l) STRICTLY INCREASING SEQUENCE OF BREAKPOINTSC C(l,I), C(2,I), VALUE AND FIRST DERIVATIVE AT XI(I), I=l,...,N+l,OF THE CUBIC SPLINE.CC****** O U T P U T ******C C(l,I), C(2,I), C(3,I), C(4,I) POLYNOMIAL COEFFICIENTS OF THE SPLINEON THE INTERVAL (XI (I), XI(I+l)) , I=l,...,N .CDATA DIAG(l),D(l) /l.,0./DO 10 M=2,N+lD (M) = XI(M) - XI(M-1)l0DIAG(M) = (C(l,M) - C(l,M-l))/D(M)DO 20 M=2,NC(2,M) = 3.* (D(M)*DIAG(M+l) + D(M+l)*DIAG(M))6.7PIECEWISE-POLYNOMIAL APPROXIMATION29120DIAG(M) = 2.*(D(M) + D(M+l))DO 30 M=2,NG = -D(M+l)/DIAG (M-1)DIAG(M) = DIAG(M) + G*D(M-1)30C(2,M) = C(2,M) + G*C(2,M-1)DO 45 M=N,2,-140C(2,M) = (C(2,M) - D(M)*C(2,M+l))/DIAG(M)RETURNENDExample 6.16: approximating a design curve by a cubic spline We are given a designcurve, a cross section of part of a car door, say, as pictured in Fig.
6.10 a. The curve hasa slope discontinuity at x = 6.1. Measurements have been taken and end slopes havebeen estimated graphically, as indicated in Fig. 6.10a and c. The problem is to find afunction s(x) which fits the data and “looks smooth.”A solution to this problem is easily provided by cubic spline interpolation to thegiven data, using two cubic splines which join continuously, but with differing slopes, atFigure 6.10 Cubic spline approximation to a design curve.292APPROXIMATIONx = 6.1. The following FORTRAN program accomplishes this, using the subprogramsSPLINE and CALCCF discussed earlier. The program reads in the data up to x = 6.1,including the two given end slopes, and stores the calculated polynomial coefficients ofthe first six polynomial pieces inC(J,I), J = l, . .
. , 4I = 1, . . . , 6Then thc data from x = 6.1 to x = 18 are read in, together with the two end slopes, andusing SPLINE and CALCCF once again, the coefficientsC(J,I), J = l, . . . , 4I = 7, . . . , 16of the remaining 10 polynomial pieces are found. Finally, the calculated piecewise-cubicfunction s(x) is evaluated, using PCUBIC, for various values of x; some of these valuesare plotted in Fig. 6.10b. Even without the slope discontinuity, polynomial interpolationto these data would produce an “unsmooth,” i.e., oscillatory, approximation because theregion of relatively high curvature near 6.1 is followed by a rather flat and enigmaticsection (see Exercise 6.7-2).FORTRAN PROGRAM FOR CUBIC SPLINE INTERPOLATION(EXAMPLE 6.16)CPROGRAM FOR EXAMPLE 6.10PARAMETER NPlMAX = 50INTEGER I, IEND, N, Nl, N2REAL C(4, NPlMAX), FX, X, XI (NPlMAX)READ 500, Nl500 FORMAT(I2)READ 501, (XI (I),C(l,I)I=1,N1),C(2,1),C(2,Nl)501 FORMAT (2E10.3)N = N1 - 1CALL SPLINE(XI,C,N)CALL CALCCF(XI,C,N)CCREAD 500, N2IEND = N + N2READ 501, (XT(T) ,C(1,I),T=Nl,IEND),C(2,Nl),C(2,IEND)N = N2 - 1CALLSPLINE(XI(N1),C(l,N1),N)CALLCALCCF(XI(Nl),C(l,Nl),N)N = IEND - 1X = XI(1)DO 12 I=1,40FX = PCUBIC(X,XI,C,N)PRINT 600, I,X,FXFORMAT(15,Fl0.l,E20.9)600X = X + .510STOPENDWe have given here only a short introduction to piecewise-polynomialapproximation.
For more detail, see, e.g., de Boor [40].Polynomial approximation and piecewise-cubic approximation differin several important aspects which become already apparent when oneconsiders interpolation. If data are given at equally spaced points, thenpolynomial interpolation becomes increasingly poor as the number ofpoints increases, as we saw in Example 2.4. There are no such difficultieseven in cubic-spline interpolation. (Note that there are also no difficulties6.7PIECEWISE-POLYNOMIAL APPROXIMATION293in trigonometric polynomial interpolation.) Also, as the number of pointsincreases, the polynomial (and the trigonometric polynomial) becomesmore and more complex in the sense that it becomes more costly toevaluate it. Also, because of the illcondition of the power form, one has touse double precision or write the polynomial in some other form, e.g., interms of Chebyshev polynomials, when the degree exceeds 10 or so.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
