Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 26
Текст из файла (страница 26)
Horner’s method for evaluating p(a) is equivalent to solution of thebidiagonal systemBy considering (5.4). we see thatHence(5.5)The recurrence for r0 = p'(a) can be expressed as Unr = q(1:n), wherer = r(0:n - 1), soHenceThis gives, using (5.5),(5.6)5.3 T HE N EWTON FORMANDP OLYNOMIAL I NTERPOLATION109NowBy looking at the form of r and q, we find from (5.6) that(5.7)This is essentially the same form of bound as for p(a) in (5.3). Analogousbounds hold for all derivatives.5.3. The Newton Form and Polynomial InterpolationAn alternative to the monomial representation of a polynomial is the Newtonform(5.8)which is commonly used for polynomial interpolation.
The interpolation problem is to choose p so that p(a i ) = fi, i = 0:n, and the numbers ci are knownP OLYNOMIALS110as divided differences. Assuming that the points a j are distinct, the divideddifferences may be computed from a standard recurrence:c(0)(0: n) = f( 0 :n )for k = 0:n - 1f o r j = n:-1:k +1endendc = c( n )Cost: 3n 2/2 flops.Two questions are of interest: how accurate are the computedandwhat is the effect of rounding errors on the polynomial values obtained byevaluating the Newton form? To answer the first question we express therecurrence in matrix-vector form:c(0)= f ,where Lk =c( k + l ) =L k c( k ),k=0:n-1,is lower bidiagonal, withDk = diag(ones(1:k + 1), ak+1- a 0 ,a k+2 - a1, .
. . ,an - an -k -1),The analysis that follows is based on the model (2.4), and so is valid only formachines with a guard digit. With the no-guard-digit model (2.6) the boundsbecome weaker and more complicated, because of the importance of termsf l(aj - aj-k-1) in the analysis.It is straightforward to show that(5.9)where Gk = diag(ones(1:k + 1), η k , k+2,. . .
,η k , n+1), where each η ij is of theform η ij = (1 + δ1)(1 + δ2)(1 + δ3), |δi | < u. Hence(5.10)From Lemma 3.7.(5.11)5.3 T HE N EWTON F ORMANDP OLYNOMIAL I NTERPOLATION111To interpret the bound, note first that merely rounding the data (f if i (1 + δi ), |δ i | < u) can cause an error ∆c as large as eround = u|L||f|,where L = Ln-1. . . L0, so errors of at least this size are inevitable.
Since|L n-1 | . . . |L0 | > |Ln- 1 . . .L 0 | = |L|, the error in the computed divided differences can be larger than eround only if there is much subtractive cancellationin the product L = Ln-1 . . . L0. If a 0 < a 1 < . . . < an , then each Li ispositive on the diagonal and nonpositive on the first subdiagonal; therefore|L n-1|. . . |L0| = |Ln-1. . . L0| = |L|, and we have the very satisfactory bound< ( ( 1-3u ) -n - 1)|L||f|. This same bound holds if the ai are arrangedin decreasing order.To examine how well the computed Newton form reproduces the f i we“unwind” the analysis above. From (5.9) we haveBy invoking Lemma 3.7 again, we obtain(5.12)If a 0 < a1 < ...
< an then> 0 for all i, and we obtain the verysatisfactory boundAgain, the samebound holds for points arranged in decreasing order.In practice it is found that even when the computed divided differences arevery inaccurate, the computed interpolating polynomial may still reproducethe original data well. The bounds (5.11) and (5.12) provide insight into thisobserved behaviour by showing thatandcan be large only whenthere is much cancellation in the products Ln-l .
. . L0 f andrespectively.The analysis has shown that the ordering a 0 < a 1 < . . . < an yields“optimal” error bounds for the divided differences and the residual, and somay be a good choice of ordering of interpolation points. However, if theaim is to minimize |p(x) - fl(p(x))| for a given xaj, then other orderingsneed to be considered. An ordering with some theoretical support is the Lejaordering, which is defined by the equations [863, 1990](5.13a)(5.13b)For a given set of n + 1 points ai , the Leja ordering can be computed in n 2flops (see Problem 5.4).We give a numerical example to illustrate the analysis.
Let n = 16 andlet a 0 < . . . < an be equally spaced points on [-1,1] . Working in simulated112POLYNOMIALSsingle precision with u = 2-246 × 10 -8 , we computed divided differencesfor two different vectors f. Error statistics were computed by regarding thesolutions computed in double precision as exact. We define the ratios(1) For f i from the normal N(0,1) distribution the divided differencesrange in magnitude from 1 to 105. and their relative errors range from 0 (thefirst divided difference, f 0 , is always exact) to 3 × 10-7. The ratio p1 = 16.3,so (5.11) provides a reasonably sharp bound for the error in .
The relativeerrors when f is reconstructed from the computed divided differences rangebetween 0 and 3 × 10-1 (it makes little difference whet her the reconstructionis done in single or double precision). Again, this is predicted by the analysis,in this case by (5.12), because p 2 = 2 × 107. For the Leja ordering, the divideddifferences are computed with about the same accuracy, but f is reconstructedmuch more accurately, with maximum relative error 7 × 10-6 ( p 1 = 1 × 103,p2 = 8 × 104).(2) For f i = exp(a i), the situation is reversed: we obtain inaccurate divided differences but an accurate reconstruction of f.
The divided differencesrange in magnitude from 10-4 to 10-1, and their relative errors are as large as1, but the relative errors in the reconstructed f are all less than 10-7. Again,the error bounds predict this behaviour: p 1 = 6 × 108, p2 = 1.02. The Lejaordering performs similarly.The natural way to evaluate the polynomial (5.8) for a given x is by ageneralization of Horner’s met hod:q n(x ) = cnfor i = n - 1:-1:0q i ( x) = ( x- a i )q i + 1 ( x) + ciendp (x) = q 0 (x)A straightforward analysis shows that (cf.
(5.2))Hence the computedis the exact value corresponding to a polynomial withslightly perturbed divided differences. The corresponding forward error boundis5.4 N OTESANDR EFERENCES1135.4. Notes and ReferencesBackward and forward error analysis for Horner’s rule was given by Wilkinson [1088, 1963, pp. 36-37, 49-50]; our results are simply Wilkinson’s presented in a different notation. The analysis has been redone by many other authors, sometimes without reference to Wilkinson’s results. Another early reference, which gives a forward error bound only, is McCracken and Dorn [743,1964, §3.5],For more on running error analysis see §3.3.Müller [782, 1983] gives a first-order error analysis for the evaluation of thedivided difference form of a polynomial.
Olver [809, 1986] derives a posteriorierror bounds for the Horner scheme with derivatives (Algorithm 5.2), phrasingthem in terms of his relative precision notation. Stewart [939, 1971] analysessynthetic division, using a matrix-oriented approach similar to that in §5.2.The relative merits of the monomial and Chebyshev representations ofa polynomial are investigated, with respect to accuracy of evaluation, byNewbery [793, 1974] and Schonfelder and Razaz [901, 1980]. Clenshaw [212,1955] showed how Horner’s method could be extended to evaluate a polynomial expressed in the Chebyshev form p(x) =where Ti is theChebyshev polynomial of degree i. Error analysis of Clenshaw’s method, andvariations of it, are given by Gentleman [433, 1969], Newbery [792, 1973], andOliver [805, 1977], [806, 1979].
Clenshaw’s scheme can be generalized to expansions in terms of arbitrary orthogonal polynomials; see Smith [927, 1965]and Algorithm 21.8.Running error bounds for Horner’s method were included in algorithms ofKahan and Farkas [637, 196 3], [638, 196 3] without explanation. Adams [5,1967] derives the bounds and extends them to evaluation of a real polynomialat a complex argument. Algorithm 5.1 is given in [5, 1967], and also in the classic paper by Peters and Wilkinson [827, 1971], which describes many aspectsof the solution of polynomial equations. Wilkinson’s paper “The PerfidiousPolynomial” [1103, 1984] (for which he was awarded the Chauvenet Prize) ishighly recommended as a beautifully written introduction to backward erroranalysis in general and error analysis for polynomials in particular.There seems to be little work on choosing the ordering of interpolationpoints to minimize the effect of rounding errors on the construction or evaluation of the interpolating polynomial.
Werner [1075, 1984] examines experimentally the effect of different orderings on the computed value of aninterpolating polynomial at a single point, for several forms of interpolatingpolynomial.The Leja ordering, which was proposed by Leja in a 1957 paper, is analysedin detail by Reichel [863, 1990 ].
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
