Arndt - Algorithms for Programmers (523138), страница 51

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 51)

To illustrate this let f have aproot of multiplicity p at r: f (x) = (x − r) h(x) with h(r) 6= 0. Thenf 0 (x) ==10 Thisword intentionally used twice.p−1p (x − r)p−1(x − r)³ph(x) + (x − r) h0 (x)´p h(x) + (x − r) h (x)0(13.137)(13.138)CHAPTER 13. ARITHMETICAL ALGORITHMS297andF (x) = f (x)/f 0 (x) =(x − r)h(x)p h(x) + (x − r) h0 (x)(13.139)The fraction on the right hand side does not vanish at the root r.With Householder’s formula (13.118) we get (iterations for F denoted by Φ%k ):Φ2 (x) =x−Φ%2 (x) =x−Φ3 (x) =x−Φ%3 (x) =x+Φ4 (x) =x+Φ%4 (x) =x+Φ5 (x) =x+ff0(13.140a)ff0f 02 − f f 002f f 02f 02 − f f 002f 2 f 00 − 2f f 02032f − 3f f 0 f 00 + f 2 f 0003f 2 f 00 − 6f f 02036f − 6f f 0 f 00 + f 2 f 0006f f 03 + 3f 3 f 000 − 9f 2 f 0 f 0030000f f − 6f 04 + 12f f 02 f 00 − 4f 2 f 0 f 000 − 3f 2 f 00224f f 03 + 4f 3 f 000 − 24f 2 f 0 f 0030000f f − 24f 04 + 36f f 02 f 00 − 8f 2 f 0 f 000 − 6f 2 f 002(13.140b)(13.140c)(13.140d)(13.140e)(13.140f)(13.140g)The terms in the numerators and denominators of Φ%k and Φk+1 are identical up to the integral constants.Schröder’s formula (13.128), when inserting f /f 0 , becomes:¡¢f 2 f 0 f f 0 f 000 − 2f f 002 + f 02 f 00ff0%Φ (x) = x +−−(13.141)3(f f 00 − f 02 )2 (f f 00 − f 02 )¡¢¡¢f 3 f 0 2f f 03 f 00 f 000 ± .

. . − 3f 2 f 02 f 0002f 4 f 0 3f 08 f 0000 ± . . . − 36f 3 f 02 f 002 f 0002−−−576 (f f 00 − f 02 )24 (f f 00 − f 02 )f k f 0 (. . .)−... −2k−1k! (f f 00 − f 02 )Checking convergence with the above example: the iteration is:x2 − ddx= 2x2 + dx +d√√Convergence is second order (independent of p): Φ%d(1 − ²2 /2 + O(²3 )).2 ( d(1 + ²)) =Φ%2 (x)= x+xSimilar to formula 13.42 we have for the error√ 1−e)Φ%k ( d·1+e=√d·1 − ek1 + ek(13.142)(13.143)if Householder’s iteration is used (one gets d divided by the k-th iteration obtained from f ).0Using√ the Schröder’s 3rd order formula for f /f with f as above we get a nice 4th (!) order iterationfor d:d − x2(d − x2 )2+ xd2d+x(d + x2 )3√ 1 + 3e2 − 3e4 − e6√ 1−e) =Φ%d3 ( d1+e1 + 3e2 + 3e4 + e6√ 1−ce2 + 3=where c = e4 2d1+c3e + 1Φ%3 (x) = x + x(13.144)(13.145)(13.146)CHAPTER 13.

ARITHMETICAL ALGORITHMS298√In general, the (1 + a k)-th order Schröder iteration for 1/ a d obtained through f /f 0 has an order ofconvergence that exceeds the expected one by one.For f (x) = 1 − d x2 one gets (3rd order Schröder):Φ%3 (x)= x+x1 − dx2(1 − dx2 )2+x21 + dx(1 + dx2 )3(13.147)√1which is even slightly nicer, also has 4th order convergence and the error expression Φ%3 ( din 13.145.13.8.41−e1+e )is asA general schemeStarting point is the Taylor series of a function f around x0 :f (x) =∞X1 (k)f (x0 ) (x − x0 )kk!(13.148)11f (x0 ) + f 0 (x0 ) (x − x0 ) + f 00 (x0 ) (x − x0 )2 + f 000 (x0 ) (x − x0 )3 + .

. .26(13.149)k=0=Now let f (x0 ) = y0 and r be the zero (f (0) = r). We then happily expand the inverse g = f −1 around y0g(0) =∞X1 (k)g (y0 ) (0 − y0 )kk!(13.150)k=011= g(y0 ) + g 0 (y0 ) (0 − y0 ) + g 00 (y0 ) (0 − y0 )2 + g 000 (y0 ) (0 − y0 )3 + . .

.26Using x0 = g(y0 ) and r = g(0) we get(13.151)11x0 − g 0 (y0 ) f (x0 ) + g 00 (y0 ) f (x0 )2 − g 000 (y0 ) f (x0 )3 + . . .(13.152)26Remains to express the derivatives of the inverse g in terms of (derivatives of) f . Not a difficult task,note thatr=f ◦ g = idthat is:f (g(x)) = x(13.153)and derive (chain rule) to get g 0 (f (x)) f 0 (x) = 1, so g 0 (y) = f 01(x) . Derive f (g(x)) − x multiple times andset the expressions to zero (arguments y of g and x of f are omitted for readability):1 =0 =0 =0 =f 0 g00 00(13.154a)0 2 00g f +f g(13.154b)+ 3f f g + f g0 00000 00 000gfgf0 00 000 3 0000 000+ 4f g f+ 3f00 2 00(13.154c)0 2 00 0000 4 0000g + 6f f g + f gThis system of linear equations in the derivatives of g is trivially solved because it is alreadyWe obtain:1g0 =f0f 00g 00 = − 0 3f´³100 20 0003f−ffg 000 =f 05´1 ³0 00 00000 30 2 0000g 0000 =10fff−15f−fff 07´1 ³432222g 00000 =105f 00 − f 0 f 00000 − 105f 0 f 00 f 000 + 15f 0 f 00 f 0000 + 10f 0 f 0009f0(13.154d)triangular.(13.155a)(13.155b)(13.155c)(13.155d)(13.155e)CHAPTER 13.

ARITHMETICAL ALGORITHMS299Thereby equation 13.152 can be written as (omitting arguments x of f everywhere)µ¶µ´¶11f 0011 ³ 00 20 000r = x− 0 f +− 03 f 2 −3f−fff3 + . . .f2f6 f 05¢ff2f 3 ¡ 002= x−−· f 00 −· 3f − f 0 f 000 − . . .003051! f2! f3! f(13.156)which is Schröder’s iteration (equation 13.129).The [i, j]-th Padé approximant of r in f gives an iteration of order i + j + 1. Write P[i,j] for an iteration(of order i + j + 1) that is obtained using the Padé approximant [i, j].µ¶−1f0xfxff(13.157)P[0,1] (x) = x2=x−=x−=x1+0f + x f0f + x f0x f0(x f )The iterations obtained in general involve powers of x in the ‘increment’ P − x:³´23xf 2f f 0 + xf f 00 + 2xf 02x3 f 0¢P[0,2] (x) ==x− ¡(13.158)2f 2 f 0 + 2xf f 0 2 + xf 2 f 00 + 2x2 f 0 32f 2 f 0 + 2xf f 0 2 + xf 2 f 00 + 2x2 f 0 3Ã!−1µ¶−1ff2f 2 f 00ff 2 (x2 f 00 )f2= x 1++ 2 02 +=x 1++3 + (x f 0 )2xf 0(x f 0 )x f2xf 0 32(x f 0 )The [i, j]-th Padé approximant of r − x in f gives an iteration of order i + j + 2 where the incrementconsists of products of powers f, f 0 , f 00 , .

. . . Write Φ[i,j] for an iteration that is obtained using the Padéapproximant [i, j]:Φ[0,1] (x)= x−2f f 0− f f 002f 02(13.159)happens to be Householder’s third order iteration.Fourth order iterations areΦ[1,1] (x)´³222 f f 0 f 000 − 3 f f 00 + 6 f 0 f 00¢= x− 0 ¡f 2 f f 0 f 000 − 6 f f 00 2 + 6 f 0 2 f 00f(13.160)3Φ[0,2] (x)=Φ[2,0] (x)12 f f 0¢12 f 0 4 − 6 f f 0 2 f 00 + 2 f 2 f 0 f 000 − 3 f 2 f 00 2³´ −10 00000 2−2fff3f000ffx− −−f2 f012 f 0 3³´422f 6 f 0 + 3 f f 0 f 00 − f 2 f 0 f 000 + 3 f 2 f 00x−6 f 05¢f2f 3 ¡ 002f· f 00 −· 3f − f 0 f 000x− 0 −0305f2f6f´³2f 3 f 00 − f 0 f 000f2 f0f 00x − 02++ff2 f06 f 03= x− ¡===(13.161)(13.162)The last one is of course just Schröder’s fourth order iteration.In general one gets n − 2 alternative forms of iterations using the approximants [0, n − 2], [1, n − 3],.

. . [n − 3, 1]. It is not difficult to find alternative Padé based iterations like the fourth order³´3222f 6 f f 0 f 00 f 000 − 9 f f 00 + 2 f 2 f 000 + 18 f 0 f 00¢(13.163)Φ4 = x − 0 ¡f 6 f f 0 f 00 f 000 − 18 f f 00 3 + 2 f 2 f 000 2 + 18 f 0 2 f 00 2CHAPTER 13. ARITHMETICAL ALGORITHMS300that can be obtained by the [2, 2]-approximant in f /f 0 and neglecting terms containing the fourth derivative of f .The iterations obtained are fractions with numerator and denominator consisting only of terms that areproducts of integral constants and f, f 0 , f 00 , .

. . , f (n−1) . There are obviously other forms of iterations, e.g.the third order iterationsÃ!µ¶q0001fff20Φ3 (x) = x − 00 f ± f 0 − 2f f 00 = x − 00 1 ± 1 − 2 02(13.164)fffthat stems from directly solving the truncated Taylor expansion of f (r) = 0 around x.f (r) =1f (x) + f 0 (x) (r − x) + f 00 (x) (r − x)22(13.165)(For f (x) = ax2 + bx + c it gives the two solutions of the quadratic equation f (x) = 0; for other functionsone gets an iterated square root expression for the roots.)Alternative rational forms can also be obtained in a way that generalizes the the method used for multipleroots: if we emphasize the so far notationally omitted dependency from the function f as Φ {f }. Theiteration Φ {f } has fixed points where f has a root r, so x − Φ {f } (x) again has a root at r. Hence we canbuild more iterations that will converge to those roots as Φn {x − Φ {f }}o (x).

For dealing with multiplecan only be expected to haveroots we used Φ {x − Φ {f }2 }k = Φ {f /f 0 }. An iteration Φ x − Φ {f }jka kth order convergence.Similarly, one can derive alternative iterations of given order by using functions that have roots where fhas them11 . For exampleg(x):=1−11 − αf (x)where α ∈ C, α 6= 0(13.166)leads to the second order iterationΦ {g}2=x−f (x)(1 + αf (x))f 0 (x)(13.167)Using g := xf (x) leads to the alternative second order iteration.Moreover, one could use a function g and its inverse ḡ := g −1 and the corresponding iteration for f (g(x))and finally apply g to get the root: (Let r0 be the zero of f (g(x)): f (g(r0 )) = 0 if g(r0 ) = r.

r0 is whatwe get from Φ {f ◦ g}.) A simple example is g(x) = ḡ(x) = 1/x, with f = xa − d and Schröder’s formulaone gets the division-less iterations for the (inverse) a-th root of d. g subject to reasonable conditions: itmust be invertible near the root r of f .13.8.5Improvements by the delta squared processGiven a sequence of partial sums xk the so called delta squared process computes a new sequence x∗k ofextrapolated sums:2x∗k= xk+2 −(xk+2 − xk+1 )xk+2 − 2 xk+1 + xk(13.168)The method is due to Aitken. The name delta squared is due to the fact that the formula can be writtensymbolically asx∗ = x −11 Itdoes not do any harm if g has additional roots.(∆x)2(∆2 x)(13.169)CHAPTER 13.

Характеристики

Список файлов книги