Shampine, Allen, Pruess - Fundamentals of Numerical Computing (523185), страница 20
Текст из файла (страница 20)
In particular, it is not convenient when we do notknow in advance what degree is appropriate. This is the most common situation when94CHAPTER 3 INTERPOLATIONapproximating data, so an alternative form due to Newton is preferred for practicalinterpolation by polynomials.
Polynomial interpolation underlies two kinds of methods widely used for the numerical solution of ordinary differential equations, Adamsmethods and backward differentiation formulas (Gear’s methods). At each step of thenumerical solution of an initial value problem, the codes attempt to find the most appropriate degree for the underlying polynomial interpolant. For this reason such codesuse either the Newton form of the polynomial or a closely related form. Although themachinery that must be developed for the Newton form may seem formidable at first,the calculations are easy to learn.A basic tactic of numerical analysis is to estimate the error in a quantity by comparing it to a quantity believed to be more accurate. If PN(x) interpolates at thenodes {x1, .
. . ,xN} and PN+1(x) interpolates at the same nodes plus xN+l, then in suitable circumstances the latter is a better approximation to f(x) and f(x) - PN(x)PN+1(x) - PN(x). If we do not know what degree is appropriate, this suggests a wayto proceed. Start with the constant polynomial P1(x) = f1 interpolating at x1. Havingcomputed PN(x), compute PN+l(x) and use it to estimate the error of PN(x). If the estimated error is too big, increase the degree by interpolating at another node and repeat.This process is the basis of the Newton form of the interpolating polynomial.For each n, the interpolant Pn(x) is constructed as a “correction” to Pn-l(x).
Because Pn-l(x) is of degree less than n - 1 and Pn(x) is of degree at most n - 1, theirdifference must be a polynomial Qn(x) of degree at most n - 1:Pn (x) = Pn-l (x) + Q n (x).(3.21)The polynomial Pn (x) interpolates at xl ,. . . ,xn-1 just as Pn-1 (x) does, so for j =1, .
. . ,n-1,This implies that the xl,. . . ,xn-l are roots of Qn(x). Because its degree is at mostn - 1, Qn(x) must have the formQn(x) = cn (x - xl ) (x - x2 ) ··· (x - xn-l )for some constant cn. The polynomial Pn(x) also interpolates at xn: fn = Pn(xn) =Because the nodes are distinct, noneof the factors (xn - xj) can vanish, and(3.22)The relations (3.21) and (3.22) along with P1(x1) = f1 provide the Newton formof the interpolating polynomial.
The coefficient cn is called the (n - 1) st divided difference off over the points x1,. . . ,xn. A number of notations are seen. A commonone is3.3 NEWTON DIVIDED DIFFERENCE FORM95In this notationthe Newton divided difference form is(3.23)It is clear from (3.23) that the leading coefficient (the coefficient of the highestdegree term) of PN(x) is f[x 1 , . . . ,xN]. Some authors take this as the definition of the( N - 1)st divided difference.Before working some examples we present a theorem that relates an nth orderdivided difference to a pair of (n - 1)st order divided differences.
The relation leads toan algorithm for computing the cn that is computationally more convenient than (3.22).Theorem 3.3.For distinct nodes {xj} and any k > i,(3.24)andProof Let R1(x) be the polynomial of degree less than k - i that interpolates f(x)on xi+l, . . . , xk and let R2(x) be the polynomial of degree less than k - i that interpolateson xi , . .
. ,xk-l. The polynomial(3.25)has a degree at most one higher than the degrees of R1(x) and R2(x). Accordingly, itsdegree is less than k - i + l. For j = i + l, . . ., k - 1,so S(x) interpolates f(x) on xi+1 , . . . ,xk-l. Moreover, S(xi ) = fi and S(xk) = fk. ByTheorem 3.1, S(x) is the interpolating polynomial of degree less than k - i + 1 thatinterpolates f(x) on all the data. The result (3.24) simply expresses the fact that theleading coefficient of the left-hand side of (3.25) equals the leading coefficient of theright-hand side.nTo illustrate the use of this theorem, we construct a divided difference table. Suppose that three rows and columns of differences have already been computed and written in a lower triangular array as follows:96CHAPTER 3INTERPOLATIONTo add another row corresponding to the node x4, start with the data f [x 4 ] = f4.
ThenNotice the pattern in these calculations:In general, the first column of the divided difference table is xj, the second is fj,the next is the first divided difference, and so on. The table provides a convenientdevice for constructing the required divided differences: the coefficients of the interpolating polynomial are the quantities along the diagonal.Example 3.4.For the data from Example 3.2, first form the difference table:Then according to (3.23),P4(x) = 0.0 + 1.91(x - 1.82) - 0.19(x - 1.82)(x - 2.50)-0.83( x - 1.82)(x - 2.50)(x - 3.65).For computational efficiency this should be evaluated in the nested formP4 (x) = (x - 1.82){1.91 + x - 2.50)[-0.19 - 0.83(x - 3.65)]}.Of course, if you expand this out, you should get the same (except for roundoff) as thenLagrange form.There are two parts to an algorithm for calculating the Newton divided difference form of the interpolating polynomial.
The first computes the divided differencesneeded for the coefficients of PN(x). It is not necessary to save the whole table as wecan use a vector cj to save the entry in the current row j as long as we compute onediagonal at a time (and do the calculations in the correct order):3.3 NEWTON DIVIDED DIFFERENCE FORM97cN := fNfor j = N - l, . . . ,l begincj = fjfor k = j + l, . . . , N beginck := (c k - c k-l) / (xk - xj)end kend j.Once these coefficients are available, the second part of the algorithm is to evaluatePN(x) for a given x:PN := cNfor k = N - l, . . . ,l beginPN := PN * (x - xk ) + ckend k.Divided differences can be related to the derivatives of f(x) using Theorem 3.2.Application of the theorem to Pn-l(x) with x = xn leads towheremin (x1 , .
. . ,xn) << max (x1, . . . ,xn ).However, we also haveEquating the two expressions shows that(3.26)for a pointin the span of the data xl,. . . ,xn. With the efficient way of computingdivided differences just presented and (3.26), we have a way to approximate derivativesof a function f(x) known only from its values at certain points.This last result gives us a better understanding of the error estimate we used tomotivate the approach. According to Theorem 3.2,We have just seen that98CHAPTER 3INTERPOLATIONComparison of these two expressions shows that if f(N) does not vary much over thespan of the data, the error of PN(x) can be estimated by comparing it to PN+1(x).It should be appreciated that the Newton form of the interpolating polynomial wasused to obtain this result, but it is true regardless of the form used for computing thepolynomial.The Newton form (3.23) is closely related to the Taylor series of f(x) about thepoint xl :As a consequence of (3.26), the Newton form of the interpolating PN(x) becomes theTaylor polynomial of degree N - 1 when the nodes x2,.
. . ,xn all tend to xl.EXERCISES(c) in the Newton divided difference form (3.23).3.15 For the data3.17 Compute the divided difference table and P3(x) for thedata from Example 3.1. Verify that this polynomial isthe same as the one in Lagrange form.calculate P4(x)(a) in the Lagrange form (3.2),(b) using the matrix method discussed in Exercise 3.4 3.18 What is the operation count for the evaluation of thecoefficients in the Newton divided difference form of(the linear system here is small enough to be solved bythe interpolating polynomial? What is the operationhand), andcount for each evaluation of PN(x)? How does this(c) in the Newton divided difference form (3.23).compare to the Lagrange form?3.16 For the data3.19 Implement the algorithms for the Newton divided difference form.
See if you can reproduce the graphs incalculate P 5 (x)Figures 3.4 and 3.5. Try your algorithm on some of(a) in the Lagrange form (3.2),the other data sets in the exercises to test the conjecture that high degree polynomial interpolation often(b) using the matrix method discussed in Exercise 3.4results in approximations whose graphs show unde(the linear system here is small enough to be solved bysirable behavior.hand), and3.4 ASSESSING ACCURACYHow do we know when we have a good approximation? We have already seen acouple of possibilities.
One is to use (3.4), that is,Since wN is a polynomial, it is easily evaluated at any x. The derivative factor presentsproblems, however, since we certainly do not knowand probably do not even know3.4 ASSESSING ACCURACY99Figure 3.4 P6(x) from Example 3.5.f(N). Another possibility is to compare the result of interpolating on one set of nodesto that of a result of higher degree obtained by interpolating on the same set plus onemore node. A variation on this is to compare results of the same degree obtained byinterpolating on different sets of nodes.
Often the best approach is to reserve somedata and evaluate the exact error f(x) - PN(x) at these nodes. A realistic appraisalmay require that a lot of data be held in reserve, and it is far from clear how to decidewhich nodes to use for interpolation and which to keep in reserve. Usually we havesome idea about the behavior of the underlying function. A graph of the data andthe interpolant is then a great help in deciding whether the interpolant reproduces thisbehavior adequately.It is illuminating to see some examples of polynomial interpolation. The programsused to calculate the interpolating polynomials below are straightforward implementations of the Newton divided difference form. We do not provide the codes becausethey are easy to write and, as will be seen, interpolating with high degree polynomialsis generally not a good idea.Example 3.5.