f5-2

5.2EvaluationofContinuedFractions 163Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).into equation (5.1.11), and then setting z=1. Sometimes you will want to compute a function from a series representation evenwhenthecomputationis notefficient. Forexample,youmaybeusingthevalues obtainedto fit the functionto an approximatingformthat youwill use subsequently (cf.§5.8). If you are summing very large numbers of slowly convergentterms, pay attention to roundoff errors! In floating-point representation it is more accurate to sum a list ofnumbersin the orderstartingwith the smallest one,ratherthanstarting with the largest one. It is even better to groupterms pairwise, then in pairs of pairs, etc., so that all additions involve operands of comparable magnitude. CITED REFERENCES AND FURTHER READING: Goodwin, E.T. (ed.) 1961, Modern Computing Methods , 2nd ed. (New York: Philosophical Li- brary), Chapter 13 [van Wijngaarden’s transformations]. [1] Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), Chapter 3. Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), §3.6. Mathews, J., and Walker, R.L. 1970, Mathematical Methods of Physics , 2nd ed. (Reading, MA: W.A. Benjamin/Addison-Wesley), §2.3. [2] 5.2 Evaluation of Continued Fractions Continuedfractionsare oftenpowerfulways of evaluatingfunctionsthat occur in scientific applications. A continued fraction looks like this: f(x)=b0+a1 b1+a2 b2+a3 b3+a4 b4+a5 b5+···(5.2.1) Printers prefer to write this as f(x)=b0+a1 b1+a2 b2+a3 b3+a4 b4+a5 b5+··· (5.2.2) In either (5.2.1)or (5.2.2),the a’s and b’s can themselves be functions of x, usually linear or quadratic monomials at worst (i.e., constants times xor times x2). For example, the continued fraction representation of the tangent function is tanx=x 1−x2 3−x2 5−x2 7−··· (5.2.3) Continued fractions frequently converge much more rapidly than power series expansions, and in a much larger domain in the complex plane (not necessarily including the domain of convergence of the series, however). Sometimes the continued fraction converges best where the series does worst, although this is not 164 Chapter5. EvaluationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).a general rule. Blanch [1]gives a good review of the most useful convergencetests for continued fractions. Therearestandardtechniques,includingtheimportant quotient-differencealgo- rithm, for going back and forth between continued fraction approximations, power series approximations, and rational function approximations. Consult Acton [2]for an introductionto this subject, and Fike [3]for furtherdetails and references. How do you tell how far to go when evaluating a continued fraction? Unlike a series, you can’t just evaluate equation (5.2.1) from left to right, stopping when the change is small. Written in the form of (5.2.1), the only way to evaluate the continued fraction is from right to left, first (blindly!) guessing how far out tostart. This is not the right way. The right way is to use a result that relates continued fractions to rational approximations, and that gives a means of evaluating (5.2.1) or (5.2.2) from leftto right. Let f ndenote the result of evaluating (5.2.2) with coefficients through anandbn. Then fn=An Bn(5.2.4) where AnandBnare given by the following recurrence: A−1≡1 B−1≡0 A0≡b0 B0≡1 Aj=bjAj−1+ajAj−2 Bj=bjBj−1+ajBj−2 j=1,2,...,n (5.2.5) ThismethodwasinventedbyJ.Wallisin1655(!),andisdiscussedinhis Arithmetica Infinitorum [4]. You can easily prove it by induction. Inpractice,thisalgorithmhassomeunattractivefeatures: Therecurrence(5.2.5) frequently generates very large or very small values for the partial numerators and denominators AjandBj. There is thus the danger of overflow or underflow of the floating-pointrepresentation. However,therecurrence(5.2.5)islinearinthe A’sand B’s. At any point you can rescale the currently saved two levels of the recurrence, e.g., divide Aj,B j,A j−1,andBj−1all by Bj. This incidentally makes Aj=fj and is convenientfor testing whether you have gonefar enough: See if fjandfj−1 from the last iteration are as close as you would like them to be. (If Bjhappens to be zero, which can happen, just skip the renormalization for this cycle. A fancier level of optimization is to renormalize only when an overflow is imminent, saving the unnecessary divides. All this complicates the program logic.) Two newer algorithms have been proposed for evaluating continued fractions. Steed’smethod doesnotuse AjandBjexplicitly,butonlythe ratio Dj=Bj−1/B j. One calculates Djand∆fj=fj−fj−1recursively using Dj=1/(bj+ajDj−1)( 5.2.6) ∆fj=(bjDj−1)∆fj−1 (5.2.7) Steed’s method (see, e.g., [5]) avoids the need for rescaling of intermediate results. However, for certain continued fractions you can occasionally run into a situation 5.2EvaluationofContinuedFractions 165Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).where the denominator in (5.2.6) approaches zero, so that D jand∆fjare very large. The next ∆fj+1will typically cancel this large change, but with loss of accuracyinthenumericalrunningsumofthe fj’s. Itis awkwardtoprogramaround this, so Steed’s method can be recommended only for cases where you know in advance that no denominator can vanish. We will use it for a special purpose inthe routine bessik(§6.7). The best general method for evaluating continued fractions seems to be the modified Lentz’s method [6]. The need for rescaling intermediate results is avoided by using boththe ratios Cj=Aj/A j−1,D j=Bj−1/B j (5.2.8) and calculating fjby fj=fj−1CjDj (5.2.9) Fromequation(5.2.5),oneeasilyshowsthattheratiossatisfytherecurrencerelations Dj=1/(bj+ajDj−1),C j=bj+aj/C j−1 (5.2.10 ) In this algorithm there is the danger that the denominator in the expression for D j, orthe quantity Cjitself, mightapproachzero. Eitheroftheseconditionsinvalidates (5.2.10). However,ThompsonandBarnett [5]showhowtomodifyLentz’salgorithm to fix this: Just shift the offending term by a small amount, e.g., 10−30. If you work through a cycle of the algorithm with this prescription, you will see that fj+1 is accurately calculated. In detail, the modified Lentz’s algorithm is this: •Setf0=b0;i fb0=0setf0=tiny. •SetC0=f0. •SetD0=0. •Forj=1,2,... SetDj=bj+ajDj−1. IfDj=0, set Dj=tiny. SetCj=bj+aj/C j−1. IfCj=0setCj=tiny. SetDj=1/D j. Set∆j=CjDj. Setfj=fj−1∆j. If|∆j−1|<e p sthen exit. Here epsis your floating-point precision, say 10−7or10−15. The parameter tiny should be less than typical values of eps|bj|, say 10−30. The above algorithm assumes that you can terminate the evaluation of the continued fraction when |fj−fj−1|is sufficiently small. This is usually the case, but by no means guaranteed. Jones [7]gives a list of theorems that can be used to justify this termination criterion for various kinds of continued fractions. ThereisatpresentnorigorousanalysisoferrorpropagationinLentz’salgorithm. However, empirical tests suggest that it is at least as good as other methods. 166 Chapter5. EvaluationofFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).ManipulatingContinuedFractions Severalimportantpropertiesofcontinuedfractionscanbeusedtorewritethem in formsthat canspeedup numericalcomputation. An equivalencetransformation an→λa n,b n→λb n,a n+1→λa n+1 (5.2.11 ) leavesthevalueofacontinuedfractionunchanged. Byasuitablechoiceofthescale factor λyou can often simplify the form of the a’s and the b’s. Of course, you cancarryoutsuccessiveequivalencetransformations,possiblywithdifferent λ’s, on successive terms of the continued fraction. Theevenandoddparts of a continued fraction are continued fractions whose successive convergentsare f2nandf2n+1, respectively. Their main use is that they convergetwiceasfastastheoriginalcontinuedfraction,andsoiftheirtermsarenot much more complicated than the terms in the original there can be a big savings incomputation. The formula for the even part of (5.2.2) is f even=d0+c1 d1+c2 d2+··· (5.2.12 ) where in terms of intermediate variables α1=a1 b1 αn=an bnbn−1,n ≥2(5.2.13 ) we have d0=b0,c 1=α1,d 1=1+ α2 cn=−α2n−1α2n−2,d n=1+ α2n−1+α2n,n ≥2(5.2.14 ) You can find the similar formula for the odd part in the review by Blanch [1]. Often a combination of the transformations (5.2.14) and (5.2.11) is used to get the best form for numerical work. We will make frequent use of continued fractions in the next chapter. CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), §3.10. Blanch, G. 1964, SIAM Review , vol. 6, pp. 383–421. [1] Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), Chapter 11. [2] Cuyt, A., and Wuytack, L. 1987, Nonlinear Methods in Numerical Analysis (Amsterdam: North- Holland), Chapter 1. Fike,C.T.1968, ComputerEvaluationofMathematicalFunctions (EnglewoodCliffs,NJ:Prentice- Hall),§§8.2, 10.4, and 10.5. [3] Wallis,J.1695,in OperaMathematica ,vol.1,p.355,OxoniaeeTheatroShedoniano.Reprinted by Georg Olms Verlag, Hildeshein, New York (1972). [4] 5.3PolynomialsandRationalFunctions 167Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Thompson,I.J.,andBarnett,A.R.1986, JournalofComputationalPhysics ,vol.64,pp.490–509. [5] Lentz, W.J. 1976, Applied Optics , vol. 15, pp. 668–671. [6] Jones, W.B. 1973, in Pad´e Approximants and Their Applications , P.R. Graves-Morris, ed. (Lon- don: Academic Press), p. 125. [7] 5.3 Polynomials and Rational Functions A polynomial of degree N−1is represented numerically as a stored array of coefficients, c(j)with j=1,...,N. We will always take c(1)to be the constant term in the polynomial, c(N)the coefficientof xN−1; but of course other conventions are possible. There are two kinds of manipulations that you can do with a polynomial: numerical manipulations (such as evaluation), where you are given the numerical value of its argument, or algebraic manipulations, where you want totransformthe coefficientarrayin someway withoutchoosinganyparticular argument. Let’s start with the numerical. We assume that youknow enough neverto evaluatea polynomialthis way: p=c(1)+c(2)*x+c(3)*x**2+c(4)*x**3+c(5)*x**4 Come the (computer) revolution, all persons found guilty of such criminal behavior will be summarily executed, and their programs won’t be! It is a matter of taste, however, whether to write p=c(1)+x*(c(2)+x*(c(3)+x*(c(4)+x*c(5)))) or p=(((c(5)*x+c(4))*x+c(3))*x+c(2))*x+c(1) If the number of coefficients is a large number n, one writes p=c(n) do11j=n-1,1,-1 p=p*x+c(j) enddo 11 Another useful trick is for evaluating a polynomial P(x)and its derivative dP(x)/dxsimultaneously: p=c(n) dp=0. do11j=n-1,1,-1 dp=dp*x+pp=p*x+c(j) enddo 11 which returns the polynomial as pand its derivative as dp. The above trick, which is basically synthetic division [1,2], generalizes to the evaluation of the polynomial and nd-1of its derivatives simultaneously: