f5-5

172 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).Ofcourseyoushouldcalculaterepeatedsubexpressions,like c/dord/c,onlyonce. Complex square root is even more complicated, since we must both guard intermediate results, and also enforce a chosen branch cut (here taken to be the negative real axis). To take the square root of c+id, first compute w≡  0 c=d=0 /radicalbig |c|/radicalBigg 1+/radicalbig 1+(d/c )2 2|c|≥|d| /radicalbig |d|/radicalBigg |c/d|+/radicalbig 1+(c/d )2 2|c|<|d|(5.4.6 ) Then the answer is √ c+id=  0 w=0 w+i/parenleftbiggd 2w/parenrightbigg w/negationslash=0,c≥0 |d| 2w+iw w /negationslash=0,c< 0,d≥0 |d| 2w−iw w /negationslash=0,c< 0,d< 0(5.4.7 ) CITED REFERENCES AND FURTHER READING: Midy, P., andYakovlev, Y. 1991, Mathematics and Computers inSimulation , vol. 33, pp. 33–49. Knuth,D.E.1981, SeminumericalAlgorithms ,2nded.,vol.2of TheArtofComputerProgramming (Reading, MA: Addison-Wesley) [see solutions to exercises 4.2.1.16 and 4.6.4.41]. 5.5 Recurrence Relations and Clenshaw’s Recurrence Formula Many useful functions satisfy recurrence relations, e.g., (n+1 )Pn+1(x)=( 2n+1 )xPn(x)−nPn−1(x) (5.5.1) Jn+1(x)=2n xJn(x)−Jn−1(x) (5.5.2) nEn+1(x)=e−x−xEn(x) (5.5.3) cosnθ=2 c o sθcos(n−1)θ−cos(n−2)θ (5.5.4) sinnθ=2 c o sθsin(n−1)θ−sin(n−2)θ (5.5.5) wherethefirstthreefunctionsareLegendrepolynomials,Besselfunctionsofthefirst kind, and exponential integrals, respectively. (For notation see [1].) These relations 5.5RecurrenceRelationsandClenshaw’sRecurrenceFormula 173Sample 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).are useful for extending computationalmethods from two successive values of nto other values, either larger or smaller. Equations(5.5.4)and(5.5.5)motivateustosayafewwordsabouttrigonometric functions. If yourprogram’srunningtime is dominatedbyevaluatingtrigonometric functions,youareprobablydoingsomethingwrong. Trigfunctionswhoseargumentsform a linear sequence θ=θ 0+nδ,n=0,1,2,..., are efficiently calculated by the following recurrence, cos(θ+δ)=c o sθ−[αcosθ+βsinθ] sin(θ+δ)=s i nθ−[αsinθ−βcosθ](5.5.6 ) whereαandβare the precomputed coefficients α≡2s i n2/parenleftbiggδ 2/parenrightbigg β≡sinδ (5.5.7 ) The reason for doing things this way, rather than with the standard (and equivalent) identities for sums of angles, is that here αandβdo not lose significance if the incrementalδis small. Likewise, the adds in equation (5.5.6) should be done in the order indicated by square brackets. We will use (5.5.6) repeatedly in Chapter12, when we deal with Fourier transforms. Another trick, occasionally useful, is to note that both sinθand cosθcan be calculated via a single call to tan: t≡tan/parenleftbiggθ 2/parenrightbigg cosθ=1−t2 1+t2sinθ=2t 1+t2(5.5.8 ) The cost of getting both sinand cos, if you need them, is thus the cost of tanplus 2 multiplies, 2 divides, and 2 adds. On machines with slow trig functions, this can be a savings. However, note that special treatmentis requiredif θ→±π. And also note that many modern machines have very fast trig functions; so you should not assume that equation (5.5.8) is faster without testing. Stability of Recurrences You need to be aware that recurrence relations are not necessarily stable against roundofferrorin the directionthat youproposeto go (eitherincreasing nor decreasingn). A three-term linear recurrence relation yn+1+anyn+bnyn−1=0,n =1,2,... (5.5.9 ) hastwolinearlyindependentsolutions, fnandgnsay. Onlyoneofthesecorresponds to the sequence of functions fnthat you are trying to generate. The other one gn maybe exponentiallygrowingin the directionthat you want to go, or exponentially damped,orexponentiallyneutral(growingordyingassomepowerlaw,forexample).If it is exponentiallygrowing, then the recurrence relation is of little or no practical use in that direction. This is the case, e.g., for (5.5.2) in the direction of increasing n, whenx<n. You cannot generate Bessel functions of high nby forward recurrence on (5.5.2). 174 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).To state things a bit more formally, if fn/gn→0asn→∞ (5.5.10 ) thenfniscalledthe minimalsolutionoftherecurrencerelation(5.5.9). Nonminimal solutionslikegnarecalled dominant solutions. Theminimalsolutionis unique,ifit exists, but dominant solutions are not — you can add an arbitrary multiple of fnto ag i v e ngn. You can evaluate any dominant solution by forward recurrence, but not the minimal solution . (Unfortunatelyit is sometimes the one you want.) AbramowitzandStegun(intheirIntroduction) [1]givea list ofrecurrencesthat are stable in the increasing or decreasing directions. That list does not contain all possible formulas, of course. Given a recurrence relation for some function fn(x) you can test it yourself with about five minutes of (human) labor: For a fixed x in your range of interest, start the recurrence not with true values of fj(x)and fj+1(x), but (first) with the values 1 and 0, respectively, and then (second) with 0 and 1, respectively. Generate 10 or 20 terms of the recursive sequences in the direction that you want to go (increasing or decreasing from j), for each of the two starting conditions. Look at the difference between the corresponding members ofthe two sequences. If the differences stay of order unity (absolute value less than 10, say), then the recurrence is stable. If they increase slowly, then the recurrence maybemildlyunstablebutquitetolerablyso. Iftheyincreasecatastrophically,thenthere is an exponentially growing solution of the recurrence. If you know that the function that you want actually corresponds to the growing solution, then you can keep the recurrence formula anyway e.g., the case of the Bessel function Y n(x)for increasingn, see§6.5; if you don’tknow which solution your functioncorresponds to, you must at this point reject the recurrence formula. Notice that you can do thistestbeforeyou go to the trouble of finding a numerical method for computing the two starting functions f j(x)andfj+1(x): stability is a property of the recurrence, not of the starting values. An alternative heuristic procedure for testing stability is to replace the recur- rencerelationbyasimilar onethatis linearwith constantcoefficients. Forexample, the relation (5.5.2) becomes yn+1−2γyn+yn−1=0 ( 5.5.11 ) whereγ≡n/xis treated as a constant. You solve such recurrence relations by trying solutions of the form yn=an. Substituting into the above recur- rence gives a2−2γa+1=0 ora=γ±/radicalbig γ2−1( 5.5.12 ) The recurrenceis stable if |a|≤ 1forall solutions a. This holds (as youcan verify) if|γ|≤ 1orn≤x. The recurrence(5.5.2)thus cannotbe used, startingwith J0(x) andJ1(x), to computeJn(x)for largen. Possibly you would at this point like the security of some real theorems on this subject (although we ourselves always follow one of the heuristic procedures).Here are two theorems, due to Perron [2]: TheoremA. Ifin(5.5.9)an∼anα,bn∼bnβasn→∞,andβ< 2α, then gn+1/gn∼−anα,f n+1/fn∼− (b/a )nβ−α(5.5.13 ) 5.5RecurrenceRelationsandClenshaw’sRecurrenceFormula 175Sample 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).andfnis the minimal solution to (5.5.9). Theorem B. Under the same conditions as Theorem A, but with β=2α, consider the characteristic polynomial t2+at+b=0 ( 5.5.14 ) If the rootst1andt2of (5.5.14)have distinct moduli, |t1|>|t2|say, then gn+1/gn∼t1nα,f n+1/fn∼t2nα(5.5.15 ) andfnis again the minimal solution to (5.5.9). Cases other than those in these two theorems are inconclusive for the existence of minimal solutions. (For more on the stability of recurrences, see [3].) Howdoyouproceedifthesolutionthatyoudesire istheminimalsolution? The answer lies in that old aphorism,that everycloudhas a silverlining: If a recurrence relation is catastrophically unstable in one direction, then that (undesired) solutionwill decrease very rapidly in the reverse direction. This means that you can start withanyseed values for the consecutive f jandfj+1and (when you have gone enoughsteps in the stable direction) you will convergeto the sequence of functionsthat you want, times an unknown normalization factor. If there is some other way to normalize the sequence (e.g., by a formula for the sum of the f n’s), then this can be a practical means of function evaluation. The method is called Miller’s algorithm . An exampleoftengiven [1,4]uses equation(5.5.2)in just this way,along with the normalization formula 1=J0(x)+2J2(x)+2J4(x)+2J6(x)+··· (5.5.16 ) Incidentally, there is an important relation between three-term recurrence relations and continuedfractions . Rewrite the recurrencerelation (5.5.9) as yn yn−1=−bn an+yn+1/yn(5.5.17 ) Iterating this equation, starting with n,g i v e s yn yn−1=−bn an−bn+1 an+1−··· (5.5.18 ) Pincherle’s Theorem [2]tells us that (5.5.18) converges if and only if (5.5.9) has a minimalsolution fn,inwhichcase itconvergesto fn/fn−1. Thisresult,usuallyfor thecasen=1andcombinedwithsomewaytodetermine f0,underliesmanyofthe practicalmethodsforcomputingspecial functionsthat we givein the nextchapter. 176 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).Clenshaw’sRecurrence Formula Clenshaw’s recurrence formula [5]is an elegant and efficient way to evaluate a sum of coefficients times functions that obey a recurrence formula, e.g., f(θ)=N/summationdisplay k=0ckcoskθorf(x)=N/summationdisplay k=0ckPk(x) Here is how it works: Suppose that the desired sum is f(x)=N/summationdisplay k=0ckFk(x)( 5.5.19 ) and thatFkobeys the recurrence relation Fn+1(x)=α(n,x )Fn(x)+β(n,x )Fn−1(x)( 5.5.20 ) for some functions α(n,x )andβ(n,x ). Now define the quantities yk(k= N,N−1,..., 1)by the following recurrence: yN+2=yN+1=0 yk=α(k,x )yk+1+β(k+1,x)yk+2+ck (k=N,N−1,..., 1)(5.5.21 ) If you solve equation (5.5.21) for ckon the left, and then write out explicitly the sum (5.5.19), it will look (in part) like this: f(x)=··· +[y8−α(8,x)y9−β(9,x)y10]F8(x) +[y7−α(7,x)y8−β(8,x)y9]F7(x) +[y6−α(6,x)y7−β(7,x)y8]F6(x) +[y5−α(5,x)y6−β(6,x)y7]F5(x) +··· +[y2−α(2,x)y3−β(3,x)y4]F2(x) +[y1−α(1,x)y2−β(2,x)y3]F1(x) +[c0+β(1,x)y2−β(1,x)y2]F0(x)(5.5.22 ) Notice that we have added and subtracted β(1,x)y2in the last line. If you examine the terms containinga factorof y8in (5.5.22),youwill findthat theysum to zeroas a consequence of the recurrence relation (5.5.20); similarly all the other yk’s down throughy2. The only surviving terms in (5.5.22) are f(x)=β(1,x)F0(x)y2+F1(x)y1+F0(x)c0 (5.5.23 ) 5.5RecurrenceRelationsandClenshaw’sRecurrenceFormula 177Sample 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).Equations (5.5.21) and (5.5.23) are Clenshaw’s recurrence formula for doing the sum (5.5.19): You make one pass down through the yk’s using (5.5.21); when you have reachedy2andy1you apply (5.5.23) to get the desired answer. Clenshaw’s recurrence as written above incorporates the coefficients ckin a downward order, with kdecreasing. At each stage, the effect of all previous ck’s is “remembered” as two coefficients which multiply the functions Fk+1andFk (ultimatelyF0andF1). If the functions Fkare small when kis large,andif the coefficientsckare small when kissmall, then the sum can be dominated by small Fk’s. In this case the remembered coefficients will involve a delicate cancellation and there can be a catastrophic loss of significance. An example would be to sumthe trivial series J 15(1) = 0 ×J0(1) + 0 ×J1(1) +...+0×J14(1) + 1 ×J15(1) (5.5.24 ) HereJ15, which is tiny, ends up represented as a canceling linear combination of J0andJ1, which are of order unity. The solution in such cases is to use an alternative Clenshaw recurrence that incorporatesck’s in an upward direction. The relevant equations are y−2=y−1=0 ( 5.5.25 ) yk=1 β(k+1,x)[yk−2−α(k,x )yk−1−ck], (k=0,1,...,N −1) ( 5.5.26 ) f(x)=cNFN(x)−β(N,x )FN−1(x)yN−1−FN(x)yN−2 (5.5.27 ) The rare case where equations (5.5.25)–(5.5.27) should be used instead of equations (5.5.21) and (5.5.23) can be detected automatically by testing whether the operands in the first sum in (5.5.23) are opposite in sign and nearly equal in magnitude. Other than in this special case, Clenshaw’s recurrence is always stable,independent of whether the recurrence for the functions F kis stable in the upward or downward direction. 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), pp. xiii, 697. [1] Gautschi, W. 1967, SIAM Review , vol. 9, pp. 24–82. [2] Lakshmikantham,V.,andTrigiante,D.1988, TheoryofDifferenceEquations:NumericalMethods and Applications (San Diego: Academic Press). [3] Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), pp. 20ff. [4] Clenshaw, C.W. 1962, Mathematical Tables ,vol. 5, NationalPhysical Laboratory (London: H.M. Stationery Office). [5] Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), §4.4.3, p. 111. Goodwin, E.T. (ed.) 1961, Modern Computing Methods , 2nd ed. (New York: Philosophical Li- brary), p. 76. 178 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).5.6 Quadratic and Cubic Equations Therootsofsimplealgebraicequationscanbeviewedasbeingfunctionsofthe equations’ coefficients. We are taught these functions in elementary algebra. Yet, surprisingly many people don’t know the right way to solve a quadratic equation with two real roots, or to obtain the roots of a cubic equation. There are two ways to write the solution of the quadratic equation ax2+bx+c=0 ( 5.6.1 ) with real coefficients a,b,c, namely x=−b±√ b2−4ac 2a(5.6.2 ) and x=2c −b±√ b2−4ac(5.6.3 ) If you use either(5.6.2)or(5.6.3) to get the two roots, you are asking for trouble: If eitheraorc(or both) are small, then one of the roots will involvethe subtraction ofbfroma verynearlyequal quantity(the discriminant);you will get that rootvery inaccurately. The correct way to compute the roots is q≡−1 2/bracketleftBig b+sgn(b)/radicalbig b2−4ac/bracketrightBig (5.6.4 ) Then the two roots are x1=q aandx2=c q(5.6.5 ) If the coefficients a,b,c, are complexrather than real, then the aboveformulas still hold, except that in equation (5.6.4) the sign of the square root should be chosen so as to make Re(b*/radicalbig b2−4ac)≥0( 5.6.6 ) where Re denotes the real part and asterisk denotes complex conjugation. Apropos of quadratic equations, this seems a convenient place to recall that the inverse hyperbolic functions sinh−1and cosh−1are in fact just logarithms of solutions to such equations, sinh−1(x)= l n/parenleftbig x+/radicalbig x2+1/parenrightbig (5.6.7 ) cosh−1(x)=±ln/parenleftbig x+/radicalbig x2−1/parenrightbig (5.6.8 ) Equation (5.6.7)is numericallyrobust for x≥0. For negativex, use the symmetry sinh−1(−x)=−sinh−1(x). Equation (5.6.8) is of course valid only for x≥1. Since FORTRAN mysteriously omits the inverse hyperbolic functions from its list of intrinsic functions, equations (5.6.7)–(5.6.8)are sometimes quite essential.