f5-1

Sample 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).Chapter 5. Evaluation of Functions 5.0 Introduction Thepurposeofthischapteristoacquaintyouwithaselectionofthetechniques that are frequently used in evaluating functions. In Chapter 6, we will apply and illustrate these techniques by giving routines for a variety of specific functions. The purposes of this chapter and the next are thus mostly in harmony, but thereis nevertheless some tension between them: Routines that are clearest and most illustrative of the general techniques of this chapter are not always the methods of choice for a particular special function. By comparing this chapter to the next one, you should get some idea of the balance between “general” and “special” methods that occurs in practice. Insofar as that balance favors general methods, this chapter should give you ideas about how to write your own routine for the evaluation of a function which, while “special” to you, is not so special as to be included in Chapter 6 or thestandard program libraries. CITED REFERENCES AND FURTHER READING: Fike,C.T.1968, ComputerEvaluationofMathematicalFunctions (EnglewoodCliffs,NJ:Prentice- Hall). Lanczos, C. 1956, Applied Analysis ; reprinted 1988 (New York: Dover), Chapter 7. 5.1 Series and Their Convergence Everybodyknowsthatananalyticfunctioncanbeexpandedintheneighborhood of a point x0in a power series, f(x)=∞/summationdisplay k=0ak(x−x0)k(5.1.1 ) Such series are straightforward to evaluate. You don’t, of course, evaluate the kth powerof x−x0abinitioforeachterm;ratheryoukeepthe k−1stpowerandupdate it with a multiply. Similarly, the form of the coefficients ais often such as to make use of previouswork: Terms like k!or(2k)!can be updatedin a multiplyor two. 159 160 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).How do you know when you have summed enough terms? In practice, the terms had better be getting small fast, otherwise the series is not a good techniqueto use in the first place. While not mathematically rigorous in all cases, standard practice is to quit when the term you have just added is smaller in magnitude than some small /epsilon1times the magnitude of the sum thus far accumulated. (But watch out if isolated instances of a k=0are possible!). A weakness of a power series representation is that it is guaranteed notto converge farther than that distance from x0at which a singularity is encountered in the complex plane . This catastrophe is not usually unexpected: When you find a power series in a book (or when you work one out yourself), you will generallyalso know the radius of convergence. An insidious problem occurs with series that converge everywhere (in the mathematical sense), but almost nowhere fast enough to be useful in a numerical method. Two familiar examples are the sine functionand the Bessel function of the first kind, sinx= ∞/summationdisplay k=0(−1)k (2k+1 ) !x2k+1(5.1.2 ) Jn(x)=/parenleftBigx 2/parenrightBign∞/summationdisplay k=0(−1 4x2)k k!(k+n)!(5.1.3 ) Both of these series converge for all x. But both don’t even start to converge until k/greatermuch|x|; before this, their terms are increasing. This makes these series useless for large x. Accelerating the Convergence of Series There are several tricks for accelerating the rate of convergenceof a series (or, equivalently, of a sequence of partial sums). These tricks will notgenerally help in cases like (5.1.2) or (5.1.3)while the size of the terms is still increasing. For serieswith terms of decreasing magnitude, however, some accelerating methods can be startlinglygood. Aitken’s δ 2-processissimplyaformulaforextrapolatingthepartial sums of a series whose convergenceis approximatelygeometric. If Sn−1,S n,S n+1 are three successive partial sums, then an improved estimate is S/prime n≡Sn+1−(Sn+1−Sn)2 Sn+1−2Sn+Sn−1(5.1.4 ) You can also use (5.1.4) with n+1and n−1replaced by n+pand n−p respectively, for any integer p. If you form the sequence of S/prime i’s, you can apply (5.1.4) a second time to thatsequence, and so on. (In practice, this iteration will only rarely do much for you after the first stage.) Note that equation (5.1.4) should be computed as written; there exist algebraically equivalent forms that are muchmore susceptible to roundoff error. Foralternating series (where the terms in the sum alternate in sign), Euler’s transformation canbeapowerfultool. Generallyitisadvisabletodoasmallnumber 5.1SeriesandTheirConvergence 161Sample 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).oftermsdirectly,throughterm n−1say,thenapplythetransformationtotherestof the series beginning with term n. The formula (for neven) is ∞/summationdisplay s=0(−1)sus=u0−u1+u2...−un−1+∞/summationdisplay s=0(−1)s 2s+1[∆sun](5.1.5 ) Here ∆is theforward difference operator , i.e., ∆un≡un+1−un ∆2un≡un+2−2un+1+un ∆3un≡un+3−3un+2+3un+1−unetc.(5.1.6 ) Of course you don’t actually do the infinite sum on the right-hand side of (5.1.5), but only the first, say, pterms, thus requiringthe first pdifferences(5.1.6) obtained from the terms starting at un. Euler’s transformation can be applied not only to convergent series. In some casesitwillproduceaccurateanswersfromthefirsttermsofaseriesthatisformally divergent. It is widely used in the summation of asymptotic series. In this caseit is generally wise not to sum farther than where the terms start increasing in magnitude;andyoushoulddevisesomeindependentnumericalcheckthattheresults are meaningful. There is an elegant and subtle implementation of Euler’s transformation due to van Wijngaarden [1]: It incorporates the terms of the original alternating series one at a time, in order. For each incorporationit eitherincreases pby 1, equivalent to computing one further difference (5.1.6), or else retroactively increases nby 1, without havingto redoall the differencecalculations based on the old nvalue! The decision as to which to increase, norp, is taken in such a way as to make the convergence most rapid. Van Wijngaarden’s technique requires only one vector of saved partial differences. Here is the algorithm: SUBROUTINE eulsum(sum,term,jterm,wksp) INTEGER jtermREAL sum,term,wksp(jterm) Workspace, provided by the calling program. Incorporates into sumthejterm’th term, with value term, of an alternating series. sum is input as the previous partial sum, and is output as the new partial sum. The first callto this routine, with the first termin the series, should be with jterm=1. On the second call,termshould be set to the second term of the series, with sign opposite to that of the first call, and jtermshould be 2. And so on. INTEGER j,ntermREAL dum,tmpSAVE nterm if(jterm.eq.1)then Initialize: nterm=1 Number of saved differences in wksp. wksp(1)=term sum=0.5*term Return first estimate. else tmp=wksp(1)wksp(1)=term do 11j=1,nterm-1 Update saved quantities by van Wijngaarden’s algo- rithm. dum=wksp(j+1) wksp(j+1)=0.5*(wksp(j)+tmp) tmp=dum 162 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).enddo 11 wksp(nterm+1)=0.5*(wksp(nterm)+tmp) if(abs(wksp(nterm+1)).le.abs(wksp(nterm)))then Favorable to increase p, sum=sum+0.5*wksp(nterm+1)nterm=nterm+1 and the table becomes longer. else Favorable to increase n, sum=sum+wksp(nterm+1) the table doesn’t become longer. endif endif return END The powerful Euler technique is not directly applicable to a series of positive terms. Occasionallyitisusefultoconvertaseriesofpositivetermsintoanalternating series,justsothattheEulertransformationcanbeused! VanWijngaardenhasgiven a transformation for accomplishing this [1]: ∞/summationdisplay r=1vr=∞/summationdisplay r=1(−1)r−1wr (5.1.7 ) where wr≡vr+2v2r+4v4r+8v8r+··· (5.1.8 ) Equations (5.1.7)and (5.1.8)replace a simple sum by a two-dimensionalsum, each term in (5.1.7) being itself an infinite sum (5.1.8). This may seem a strange way tosave on work! Since, however,the indices in (5.1.8) increase tremendouslyrapidly, aspowersof2,itoftenrequiresonlyafewtermstoconverge(5.1.8)toextraordinary accuracy. You do, however, need to be able to compute the v r’s efficiently for “random” values r. The standard “updating” tricks for sequential r’s, mentioned above following equation (5.1.1), can’t be used. Actually,Euler’stransformationisa specialcaseofa moregeneraltransforma- tion of power series. Suppose that some known function g(z)has the series g(z)=∞/summationdisplay n=0bnzn(5.1.9 ) and that you want to sum the new, unknown, series f(z)=∞/summationdisplay n=0cnbnzn(5.1.10 ) Then it is not hard to show (see [2]) that equation (5.1.10)can be written as f(z)=∞/summationdisplay n=0[∆(n)c0]g(n) n!zn(5.1.11 ) which often convergesmuch more rapidly. Here ∆(n)c0is the nth finite-difference operator (equation 5.1.6), with ∆(0)c0≡c0, and g(n)is the nth derivative of g(z). The usual Euler transformation (equation 5.1.5 with n=0) can be obtained, for example, by substituting g(z)=1 1+z=1−z+z2−z3+··· (5.1.12 ) 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