Conte, de Boor - Elementary Numerical Analysis. An Algorithmic Approach (523140), страница 9
Текст из файла (страница 9)
For if we use the abbreviationthen the calculations of Algorithm 2.2 readIn particular, the number d i,k-1 is not used any further once dik has beencalculated, so that we can safely store d ik over d i,k-1 .Algorithm 2.3: Calculation of the coefficients for the Newton formulaGiven the n + 1 distinct points x0, . . . , xn, and, correspondingly, thenumbers f(x0), . . . , f(xn), with f(xi ) stored in di , i = 0, . . .
, n.ThenExample 2.4Let f(x) = (1 + x2)-1. For n = 2, 4, . . . , 16, calculate the polynomialPn(x) of degree < n which interpolates f(x) at the n + 1 equally spaced pointsThen estimate the maximum interpolation erroron the interval [-5, 5] by computing2.3THE DIVIDED-DIFFERENCE TABLE45whereThe FORTRAN program below uses Algorithms 2.1 and 2.3 to solve this problem.FORTRAN PROGRAM FOR EXAMPLE 2.4C PROGRAM FOR EXAMPLE 2.4INTEGER I,J,K,N,NP1REALD(17),ERRMAX,H,PNOFY,X(17),YC POLYNOMIAL INTERPOLATION AT EQUALLY SPACED POINTS TO THE FUNCTIONF(Y) = l./(l. + Y*Y)CPRINT 600N',5X,'MAXIMUM ERROR')600 FORMAT('1DO 40 N=2,16,2NP1 = N+1H = 10./FLOAT(N)DO 10 I=1,NP1X(I) = FLOAT(I-1)*H - 5.D(I) = Fix(I))10CONTINUECCALCULATE DIVIDED DIFFERENCES BY ALGORITHM 2.3DO 20 K=1,NDO 20 I=1,NP1-RD(I) = (D(I+1) - D(I))/(X(I+K) - X(I))CONTINUE20ESTIMATE MAXIMUM INTERPOLATION ERROR ON (-5,5)CERRMAX = 0.DO 30 J=1,101Y = FLOAT(J-1)/10. - 5.CCALCULATE PN(Y) BY ALGORITHM 2.1PNOFY = D(1)DO 29 K=2,NP1PNOFY = D(K) + (Y - X(K))*PNOFY'29CONTINUEERRMAX = MAX(ABS(F(Y) - PNOFY) , ERRMAX)CONTINUE30PRINT 630, N,ERRMAXFORMAT(I5,El8.7)63040 CONTINUE .STOPENDCOMPUTER OUTPUT FOR EXAMPLE 2.4N246810121416MAXIMUM ERROR6.4615385E - 014.3813387E - 016.1666759E - 011.0451739E + 001.9156431E + 003.6052745E + 007.192008OE + 0014051542E + 01Note how the interpolation error soon increases with increasing degree even though we usemore and more information about the function f(x) in our interpolation process.
This isbecause we have used uniformly spaced interpolation points; see Exercise 6.1-12 and Eq.(6.20).46INTERPOLATION BY POLYNOMIALSEXERCISES2.3-l From a table of logarithms we obtain the following values of log x at the indicatedtabular points.xlog x1.01.52.03.03.54.00.00.176090.301030.477 120.544070.60206Form a divided-difference table based on these values.2.3-2 Using the divided-difference table in Exercise 2.3-1, interpolate for the followingvalues: log 2.5, log 1.25, log 3.25. Use a third-degree interpolating polynomial in its Newtonform.2.3-3 Estimate the error in the result obtained for log 2.5 in Exercise 2.3-2 by computing thenext term in the interpolating polynomial.
Also estimate it by comparing the approximationfor log 2.5 with the sum of log 2 and the approximation for log 1.25.2.3-4 Derive the formulaThen use it to interpret the Nested Multiplication Algorithm 2.1, applied to the polynomial(2.10), as a way to calculate p[z, x0, . . . , xn-1], p[z, x0, . .
. , xn-2], . . . , p[z, x0] and P[z], i.e.,as a way to get another diagonal in the divided difference table for p(x).2.3-5 By Exercise 2.2-3, the polynomial of degree < k which interpolates a function f(x) atx0, . . . , xk is f(x) itself if f(x) is a polynomial of degree < k. This fact may be used to checkthe accuracy of the computed interpolating polynomial. Adapt the FORTRAN program givenin Example 2.4 to carry out such a check as follows: For n = 4, 8, 12, .
. . , 32, find thepolynomial pn(x) of degree < n which interpolates the functionat0,1,2, . . . ,n. Then estimatewherethe yi's are a suitably large number of points in [0, n] .2.3-6 Prove that the first derivative p'2(x) of the parabola interpolating f(x) at x0 < xl < x2 isequal to the straight line which takes on the value f[xi-1, xi] at the point (xi-1 + xi) /2, fori = 1, 2.
Generalize this to describe p'n(x) as the interpolant to dataforin case pn(x) interpolates f(x) at x0 < x1 < · · · < xn.appropriate*2.4 INTERPOLATION AT AN INCREASING NUMBER OFINTERPOLATION POINTSConsider now the problem of estimating f(x) at a pointusingpolynomial interpolation at distinct points x0, x1, x2, . . . .
With pk(x) thepolynomial of degree < k which interpolates f(x) at x0, . . . , xk, we calcuuntil, so we hope, the differencelate successivelybetweenandis sufficiently small. The Newton form for the*2.4INTERPOLATION AT AN INCREASING NUMBER OF INTERPOLATION POINTS47interpolating polynomialwithis expressly designed for such calculations. If we knowthen we can calculateandAlgorithm 2.4: Interpolation using an increasing number of interpolationpoints Given distinct points x 0 , x 1 , x 2 , . .
. and the value!f(x0), f(x1), f(x2), . . . of a function f(x) at these points. Also, given apointFor k = 0, 1, 2, . . . , until satisfied, do:This algorithm generates the entries of the divided-difference table forf(x) at x0 , x1 , x2 , . . . a diagonal at a time. During the calculation ofthe upward diagonal emanating from f[xk+1] is calculated up toand including the number f[x0, .
. . , xk+1], using the number f[xk+1] =f(xk+1) and the previously calculated entries f[xk], f[xk-1, xk],. . . , f[x0, . . . , xk] in the preceding diagonal. Hence, even if only themost recently calculated diagonal is saved (in a FORTRAN program, say),the algorithm provides incidentally the requisite coefficients for the Newton form for pk+1(x) with centers xk+1, . . . , x1:(2.14)Example 2.5 We apply Algorithm 2.4 to the problem of Examples 2.2 and 2.3, usingx0 = 1, x1 = 4, x2 = 6, and in addition, x3 = 0.
For this example,We getNext, with K[x1] = 1.5727, we get0.0006, and withwe get1.5724.48INTERPOLATION BY POLYNOMIALSAdding the point x 2 = 6, we have K[x2 ] = 1.5751; hence K[x1 , x 2 ] = 0.0012,K[x0, x1, x2] = 0.00012; therefore, asthe number calculated earlier in Example 2.3. To check the error for this approximationto K(3.5), we add the point x 3 = 0. With K[x3 ] = 1.5708, we compute K[x2 , x 3 ] =0.000717, K[x1, x2, x3] = 0.000121, K[x0, x1, x2, x3] = - 0.000001, and get, with= (-2.5)(-1.25) = 3.125, thatindicating that 1.5722 or 1.5723 is probably the value of K(3.5) to within the accuracy ofthe given values of K(x).These calculations, if done by hand, are conveniently arranged in a table as shownin Fig.
2.2, which also shows how Algorithm 2.4 gradually builds up the divided-difference table.We have listed below a FORTRAN FUNCTION, called TABLE,which uses Algorithm 2.4 to interpolate in a given table of abscissas andordinates X(I), F(I), I = 1, . . . , NTABLE, with F(I) = f(X(I)), and X(1)< X(2) < · · · , in order to find a good approximation to f(x) at x =XBAR.
The program generates p0 (XBAR), p1 (XBAR), . . . , untilwhere TOL is a given e r r o r r e q u i r e m e n t , o r u n t i l k + 1 =min(20, NTABLE), and then returns the number pk (XBAR). The sequencex0, x1, x2, . . . of points of interpolation is chosen from the tabular pointsX(1), X(2), . . . , X(NTABLE) as follows: If X(I) < XBAR < X(I + 1),then x0 = X(I + 1), x1 = X(I), x2 = X(I + 2), x3 = X(I - 1), . . . , exceptnear the beginning or the end of the given table, where eventually onlypoints to the right or to the left of XBAR are used. To protect the program(and the user!) against an unreasonable choice for TOL, the programshould be modified so as to terminate also if and when the successivedifferences |p k+1 (XBAR) - p k(XBAR)| begin to increase as k increases.(See also Exercise 2.4-1.)Figure 2.2*2.4INTERPOLATION AT AN INCREASING NUMBER OF INTERPOLATION POINTS49FORTRAN SUBPROGRAM FOR INTERPOLATION IN AFUNCTION TABLEREAL FUNCTION TABLE (XBAR, X, F, NTABLE, TOL, I'FLAG )C RETURNS AN INTERPOLATED VALUE TABLE AT XBAR FOR THE FUNCTIONC TABULATED AS (X(I),F(I)), I=l,...,NTABLE.INTEGER IFLAG,NTABLE,J,NEXT,NEXTL,NEXTRREALF(NTABLE),TOL,X(NTABLE),XBAR,A(20),ERROR,PSIK,XK(20)C****** I N P U T ******C XBAR POINT AT WHICH TO INTERPOLATE .C X(I), F(I), I=1 ,...,NTABLE CONTAINS THE FUNCTION TABLE .CA S S U M P T I O N ...
X IS ASSUMED TO BE INCREASING.)C NTABLE NUMBER OF ENTRIES IN FUNCTION TABLE.C TOL DESIRED ERROR BOUND .C****** O U T P U T ******C TABLE THE INTERPOLATED FUNCTION VALUE .C IFLAG AN INTEGER,C=l , SUCCESSFUL EXECUTION ,C=2 , UNABLE TO ACHIEVE DESIRED ERROR IN 20 STEPS,C=3 , XBAR LIES OUTSIDE OF TABLE RANGE. CONSTANT EXTRAPOLATION ISCUSED.C****** M E T H O D ******C A SEQUENCE OF POLYNOMIAL INTERPOLANTS OF INCREASING DEGREE IS FORMEDC USING TABLE ENTRIES ALWAYS AS CLOSE TO XBAR AS POSSIBLE. EACH INC TERPOLATED VALUE IS OBTAINED FROM THE PRECEDING ONE BY ADDITION OF AC CORRECTION TERM (AS IN THE NEWTON FORMULA). THE PROCESS TERMINATESC WHEN THIS CORRECTION IS LESS THAN TOL OR, ELSE, AFTER 20 STEPS.CCLOCATE XBAR IN THE X-ARRAY.IF (XBAR .GE. X(l) .AND.
XBAR .LE. X(NTABLE)) THENDO 10 NEXT=2,NTABLEIF (XBAR .LE. X(NEXT))GO TO 12CONTINUE10END IFIF (XBAR .LT. X(1)) THENTABLE = F(1)ELSETABLE = F(NTABLE)END IFPRINT 610,XBAR610 FORMAT(E16.7,' NOT IN TABLE RANGE.')IFLAG = 3RETURN12 XK(1) = X(NEXT)NEXTL = NEXT-lNEXTR = NEXT+1A(1) = F(NEXT)TABLE = A(1)PSIK = 1.USE ALGORITHM 2.4, WITH THE NEXT XK ALWAYS THE TABLECCENTRY NEARESTXBAR OF THOSE NOT YET USED.KP1MAX = MIN(20,NTABLE)DO 20 KP1=2,KP1MAXIF (NEXTL .EQ.
0) THENNEXT = NEXTRNEXTR = NEXTR+1ELSE IF (NEXTR .GT. NTABLE) THENNEXT = NEXTLNEXTL = NEXTL-1ELSE IF (XBAR - X(NEXTL) .GT. X(NEXTR) - XBAR) THENNEXT = NEXTRNEXTR = NEXTR+1ELSENEXT = NEXTLNEXTL = NEXTL-1END IFXK(KP1) = X(NEXT)A(KP1) - F(NEXT)DO 13 J=KP1-1,1,-lA(J) = (A(J+l) - A(J))/(XK(KP1) - XK(J))13CONTINUE50INTERPOLATION BY POLYNOMIALSFOR I=1 ,...,KP1, A(I) NOW CONTAINS THE DIV.DIFF. OFF(X) OF ORDER K-I AT XK(I) ,...,XK(KP1).PSIK = PSIK*(XBAR - XK(KP1-1))ERROR = A(1)+PSIKTEMPORARY PRINTOUTCPRINT613,KP1,XK(KP1),TABLE,ERRORFORMAT(110,3El7.7)613TABLE = TABLE + ERRORIF (ABS(ERROR) .LE.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
