Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 25
Текст из файла (страница 25)
Analysis of Algorithm 4.1 andcompensated summation can also be found in Espelid [356, 1978].The earliest error analysis of summation is that of Wilkinson for recursivesummation in [1084, 1960], [1088, 1963].Pairwise summation was first discussed by McCracken and Dorn [743,1964, pp. 61-63], Babuska [35, 1969], and Linz [707, 1970]. Caprani [185, 1971]shows how to implement the method on a serial machine using temporarystorage of size [log2 n] + 1 (without overwriting the xi ).The use of compensated summation with a Runge-Kutta formula is described by Vitasek [1055, 1969]. See also Butcher [170, 1987, pp.
118-120]and the experiments of Linnainrnaa [703, 1974]. Davis and Rabinowitz [267,1984, §4.2.1] discuss pairwise summation and compensated summation in thecontext of quadrature.Problems4.1. Define and evaluate a condition number C(x) for the summation Sn (x) =When does the condition number take the value 1?4.2. (Wilkinson [1088, 1963, p. 19]) Show that the bounds (4.3) and (4.4) arenearly attainable for recursive summation. (Hint: assume u = 2- t , set n = 2r(r << t), and definex(1) = 1,x (2) = 1 - 2-t ,x (3:4) = 1 - 21-t ,x(5:8) = 1 - 22-t ,x(2r-1 + 1:2r ) = 1 - 2r - 1 -t .)PROBLEMS4.3. Let S n =order. Show that101xi be computed by recursive summation in the naturaland hence that En =satisfiesWhich ordering of the xi minimizes this bound?4.4. Let M be a floating point number so large that fl (10 + M) = M.
Whatare the possible values ofwhere= {1, 2, 3, 4, M, -M}and the sum is evaluated by recursive summation?4.5. The “±” method for computingis defined as follows: formthe sum of the positive numbers, S+, and the sum of the nonpositive numbers,S -, separately, by any method, and then form S n = S - + S+. Discuss thepros and cons of this met hod.4.6. Let {x i } be a convergent sequence with limit . Aitken’s ∆2-method(Aitken extrapolation) generates a transformed sequence {y i } defined byUnder suitable conditions (typically that {x i } is linearly convergent), the yiconverge tofaster than the xi . Which of the following expressions shouldbe used to evaluate the denominator in the formula for yi ?( a) (xi +2 - 2xi+1) + xi .( b) (x i+2 - xi+1) - (x i+1 - xi ).(c) (x i+2 + xi ) - 2x i+ l .4.7. Analyse the accuracy of the following method for evaluating4.8.
In numerical methods for quadrature and for solving ordinary differentialequation initial value problems it is often necessary to evaluate a function onan equally spaced grid of points on a range [a,b]: xi := a + ih, i = 0:n,where h = (b-a)/n. Compare the accuracy of the following ways to form xi .Assume that a and b, but not necessarily h, are floating point numbers.102SUMMATION(a) x i = x i-1 + h (x0 = a).(b) xi = a + ih.(c) xi = a(1 - i/n) + (i/n)b.Note that (a) is typically used without comment in, for example, a NewtonCotes quadrature rule or a Runge-Kutta method with fixed step size.4.9. (R ESEARCH P ROBLEM ) Priest [844, 1992 , pp.
61-62] has proved thatif |x1 | > |x2| > |x3| then compensated summation computes the sum x 1 +x2 + x3 with a relative error of order u (under reasonable assumptions on thearithmetic, such as the presence of a guard digit). He also gives the examplex 1 = 2t+ 1,x2=2 t+ 1 - 2 , x3= x4 = x5 = x 6 =-(2 t -1),for which the exact sum is 2 but compensated summation computes 0 inIEEE single precision arithmetic (t=24).
What is the smallest n for whichcompensated summation applied to x1,. . . , xn ordered by decreasing absolutevalue can produce a computed sum with large relative error?PreviousHomeChapter 5PolynomialsThe polynomial (z - 1)(z - 2) . . . (z - 20) is not a ‘difficult’ polynomial per se . . .The ‘difficulty’ with the polynomial Π(z - i) is that ofevaluating the explicit polynomial accurately.If one already knows the roots, then the polynomial can be evaluatedwithout any loss of accuracy.-J. H. WILKINSON, The Perfidious Polynomial (1984)I first used backward error analysis in connection withsimple programs for computing zeros of polynomialssoon after the PILOT ACE came into use.-J.
H. WILKINSON, The State of the Art in Error Analysis (1985)The Fundamental Theorem of Algebra asserts thatevery polynomial equation over the complex field has a root.It is almost beneath the dignity of such a majestic theoremto mention that in fact it has precisely n roots.-J. H. WILKINSON, The Perfidious Polynomial (1984)103Next104POLYNOMIALSTwo common tasks associated with polynomials are evaluation and interpolation: given the polynomial find its values at certain arguments, and given thevalues at certain arguments find the polynomial.
We consider Horner’s rulefor evaluation and the Newton divided difference polynomial for interpolation.A third task not considered here is finding the zeros of a polynomial. Much research was devoted to polynomial zero finding up until the late 1960s; indeed,Wilkinson devotes a quarter of Rounding Errors in Algebraic Processes [1088,19 6 3 ] to the topic.
Since the development of the QR algorithm for findingmatrix eigenvalues there has been less demand for polynomial zero finding,since the problem either arises as, or can be converted to (see §26.6 and [346,1995], [1007, 1994]), the matrix eigenvalue problem.5.1. Horner’s MethodThe standard met hod for evaluating a polynomialp(x) = a0 + a 1 x + .
. . + anxn(5.1)is Horner’s method (also known as Horner’s rule and nested multiplication),which consists of the following recurrence:qn(x) = anfor i = n - 1: -1:0qi(x) = xqi+ 1 (x) + aiendp(x) = q0(x)The cost is 2n flops, which is n less than the more obvious method of evaluationthat explicitly forms powers of x (see Problem 5.2).To analyse the rounding errors in Horner’s met hod it is convenient to usethe relative error counter notation <k> (see (3.9)).
We haveIt is easy to either guess or prove by induction that(5.2)1055.1 H ORNER'S M ETHODwhere we have used Lemma 3.1, and where |φ k | < ku/(1 - ku) =: γk . Thisresult shows that Horner’s method has a small backward error: the computedis the exact value at x of a polynomial obtained by making relativeperturbations of size at most γ2 n to the coefficients of p( x).A forward error bound is easily obtained: from (5.2) we have(5.3)whereThe relative error is bounded according toClearly, the factor ψ ( p , x ) can be arbitrarily large.
However, ψ(p,x) = 1 ifai > 0 for all i and x > 0, or if (-l)i ai > 0 for all i and x < 0.In a practical computation we may wish to compute an error bound alongwithThe bound (5.3) is entirely adequate for theoretical purposes and canitself be computed by Horner’s method. However, it lacks sharpness for tworeasons. First, the bound is the result of replacing each γ k by γ2 n . Second,and more importantly, it is an a priori bound and so takes no account ofthe actual rounding errors that occur. We can derive a sharper, a posterioribound by a running error analysis.For the ith step of Horner’s method we can write(5.4)where we have used both (2.4) and (2.5).
Defining=: qi + fi , we haveorHenceSince fn = 0, we have |fi | < uπi , whereWe can slightly reduce the cost of evaluating the majorizing sequence πiworking withwhich satisfies the recurrenceWe can now furnish Horner’s method with a running error bound.by106POLYNOMIALSAlgorithm 5.1. This algorithm evaluates y = fl(p(x)) by Horner’s method,It also evaluates a quantity µ such that |y - p(x)| < µ.where p(x) =y = anµ = |y|/2for i = n- 1 : - 1 : 0y = xy + aiµ = |x|µ + |y |endµ = µ( 2µ - |y|)Cost: 4n flops.It is worth commenting on the case where one or more of the ai and xis complex.
The analysis leading to Algorithm 5.1 is still valid for complexdata, but we need to remember that the error bounds for fl ( x op y) are notthe same as for real arithmetic. In view of Lemma 3.5. it suffices to replacethe last line of the algorithm by µ =An increase in speed ofthe algorithm, with only a slight worsening of the bound, can be obtained byreplacing |y| = (( Rey )2 + (Imy)2)½ by |Rey| + |Imy| (and, of course, |x|should be evaluated once and for all before entering the loop).One use of Algorithm 5.1 is to provide a stopping criterion for a polynomialzero-finder: if |fl(p ( x))| is of the same order as the error bound µ, then furtheriteration serves no purpose, for as far as we can tell, x could be an exact zero.As a numerical example, for the expanded form of p (x) = (x + 1)32 wefound in MATLAB thatand for p(x) the Chebyshev polynomial of degree 32,In these two cases. the running error bound is, respectively.
62 and 31 timessmaller than the a priori one.In another experiment we evaluated the expanded form of p(x) = (x - 2)3is simulated single precision in MATLAB (u 6 × 10-8) for 200 equally spacedpoints near x = 2. The polynomial values, the error, and the a priori andrunning error bounds are all plotted in Figure 5.1. The running error boundis about seven times smaller than the a priori one.5.2. Evaluating DerivativesSuppose now that we wish to evaluate derivatives of p. We could simplydifferentiate (5.1) as many times as necessary and apply Horner’s method to5.2 E VALUATING D ERIVATIVES107Figure 5.1. Computed polynomial values (top) and running and a priori bounds(bottom) for Horner’s method.each expression, but there is a more efficient way.
Observe that if we defineq (x) = q 1 + q 2 x + . . . + qnxn -1 r = q0,where the qi = qi (a) are generated by Horner’s method for p(a), thenp(x) = (x - a)q(x) + r.In other words, Horner’s method carries out the process of synthetic division.Clearly, p'(a) = q(a). If we repeat synthetic division recursively on q(x), wewill be evaluating the coefficients in the Taylor expansionand after a final scaling by factorials, we will obtain the derivatives of p at a .The resulting algorithm is quite short.Algorithm 5.2. This algorithm evaluates the polynomialand its first k derivatives at a, returning yi = p( i)(a), i = 0:k.y0 = any(1:k) = 0P OLYNOMIALS108for j = n- 1 : - 1 : 0for i = min(k, n - j): -1:1yi = ayi + yi- 1endy 0 = ay0 + ajendm = lfor j = 2:km = m*jyj = m*yjendcost: nk + 2(k + n) - k2/2 flops.How is the error bounded for the derivatives in Algorithm 5.2? To answerthis quest ion with the minimum of algebra, we express the algorithm in matrixnotation.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
