Chapter 5.Evaluation of Functions 5.0 Introduction NUMERICAL The purpose of this chapter is to acquaint you with a selection of the techniques that are frequently used in evaluating functions.In Chapter 6,we will apply and Cambridge 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 there is 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 、鱼 compu Press. C:THEA 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 the SCIENTIFIC standard program libraries 6 CITED REFERENCES AND FURTHER READING: Fike,C.T.1968,Computer Evaluation of Mathematical Functions(Englewood Cliffs,NJ:Prentice- Hall). r Numerical Lanczos,C.1956,Applied Analysis;reprinted 1988(New York:Dover),Chapter 7. Further Numerical (ISBN 0-521- Recipes 491069 5.1 Series and Their Convergence North Software. Everybody knows that an analytic function can be expanded in the neighborhood of a point zo in a power series, America). f)=ak(-xo) (5.1.1) k=0 Such series are straightforward to evaluate.You don't,of course,evaluate the kth power ofxr-zo ab initio for each term;rather you keep the k-1st power and update it with a multiply.Similarly,the form of the coefficients a is often such as to make use of previous work:Terms like k!or(2k)!can be updated in a multiply or two. 165
Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Chapter 5. Evaluation of Functions 5.0 Introduction The purpose of this chapter is to acquaint you with a selection of the techniques 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 there is 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 the standard program libraries. CITED REFERENCES AND FURTHER READING: Fike, C.T. 1968, Computer Evaluation of Mathematical Functions (Englewood Cliffs, NJ: PrenticeHall). Lanczos, C. 1956, Applied Analysis; reprinted 1988 (New York: Dover), Chapter 7. 5.1 Series and Their Convergence Everybody knows that an analytic function can be expanded in the neighborhood of a point x0 in a power series, f(x) = ∞ k=0 ak(x − x0) k (5.1.1) Such series are straightforward to evaluate. You don’t, of course, evaluate the kth power of x−x0 ab initio for each term; rather you keep the k −1st power and update it with a multiply. Similarly, the form of the coefficients a is often such as to make use of previous work: Terms like k! or (2k)! can be updated in a multiply or two. 165
166 Chapter 5.Evaluation of Functions 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 technique to 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 e times the magnitude of the sum thus far accumulated.(But watch out if isolated instances of ak=0 are possible!). A weakness of a power series representation is that it is guaranteed not to converge farther than that distance from xo at 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 generally also 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 function and the Bessel function of the first kind, 茶 sin (-1) +可 2k+1 (5.1.2) k=0 （-z2)* 9 =()”】 (5.1.3) Both of these series converge for all x.But both don't even start to converge until;before this,their terms are increasing.This makes these series useless for large x. IENTIFIC Accelerating the Convergence of Series 6 There are several tricks for accelerating the rate of convergence of a series(or, equivalently,of a sequence of partial sums).These tricks will not generally help in cases like (5.1.2)or(5.1.3)while the size of the terms is still increasing.For series with terms of decreasing magnitude,however,some accelerating methods can be startlingly good.Aitken's 62-process is simply a formula for extrapolating the partial 10.621 sums of a series whose convergence is approximately geometric.If Sn1,Sn,Sn+1 % are three successive partial sums,then an improved estimate is Numerical Recipes 43106 (outside S%三Sn+1- (Sn+1-Sn)2 Sn+1-2Sn Sn-1 (5.1.4) North You can also use (5.1.4)with n +1 and n-1 replaced by n+p and n-p respectively,for any integer p.If you form the sequence of Ss,you can apply (5.1.4)a second time to that sequence,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 much more susceptible to roundoff error. For alternating series (where the terms in the sum alternate in sign),Euler's transformation can be a powerful tool.Generally it is advisable to do a small
166 Chapter 5. Evaluation of Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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 technique to 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 times the magnitude of the sum thus far accumulated. (But watch out if isolated instances of ak = 0 are possible!). A weakness of a power series representation is that it is guaranteed not to converge farther than that distance from x0 at 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 generally also 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 function and the Bessel function of the first kind, sin x = ∞ k=0 (−1)k (2k + 1)!x2k+1 (5.1.2) Jn(x) = x 2 n∞ 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 |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 convergence of a series (or, equivalently, of a sequence of partial sums). These tricks will not generally help in cases like (5.1.2) or (5.1.3) while the size of the terms is still increasing. For series with terms of decreasing magnitude, however, some accelerating methods can be startlingly good. Aitken’s δ2-process is simply a formula for extrapolating the partial sums of a series whose convergence is approximately geometric. If S n−1, Sn, Sn+1 are three successive partial sums, then an improved estimate is S n ≡ Sn+1 − (Sn+1 − Sn)2 Sn+1 − 2Sn + Sn−1 (5.1.4) You can also use (5.1.4) with n + 1 and n − 1 replaced by n + p and n − p respectively, for any integer p. If you form the sequence of S i’s, you can apply (5.1.4) a second time to that sequence, 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 much more susceptible to roundoff error. For alternating series (where the terms in the sum alternate in sign), Euler’s transformation can be a powerful tool. Generally it is advisable to do a small
5.1 Series and Their Convergence 167 number of terms directly,through term n-1 say,then apply the transformation to the rest of the series beginning with term n.The formula(for n even)is ∑(-1r,=o-4+g-u-1+∑Aud (5.1.5) s= Here A is the forward difference operator,i.e., △un三un+1-ln △2un≡un+2-2un+1+un (5.1.6) △3un≡un+3-3un+2+3un+1-un etc. Of course you don't actually do the infinite sum on the right-hand side of(5.1.5). but only the first,say,p terms,thus requiring the first p differences(5.1.6)obtained from the terms starting at un. Euler's transformation can be applied not only to convergent series.In some 令 cases it will produce accurate answers from the first terms of a series that is formally divergent.It is widely used in the summation of asymptotic series.In this case it is generally wise not to sum farther than where the terms start increasing in ，%4 Press. magnitude;and you should devise some independent numerical check that the results ART Programs are meaningful. 9 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 SCIENTIFIC one at a time,in order.For each incorporation it either increases p by 1,equivalent to computing one further difference (5.1.6),or else retroactively increases n by 1, without having to redo all the difference calculations based on the old n value!The decision as to which to increase,n or p,is taken in such a way as to make the convergence most rapid.Van Wijngaarden's technique requires only one vector of COMPUTING (ISBN 18881092 saved partial differences.Here is the algorithm: #include 10-621 void eulsum(float *sum,float term,int jterm,float wksp[]) idge.org FuurrgPgoglrion Numerical Recipes 43106 Incorporates into sum the jterm'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 call to this routine,with the first term in the series,should be with jterm=1.On the second call,term should be set to the second term of the series,with sign opposite to that of the first call,and (outside jterm should be 2.And so on.wksp is a workspace array provided by the calling program, Software. dimensioned at least as large as the maximum number of terms to be incorporated. North int j; static int nterm; float tmp,dum; if (jterm ==1){ Initialize: nterm=1; Number of saved differences in wksp. *sum=0.5*(wksp[1]=term); Return first estimate. else tmp=wksp[1]: wksp[1]=term; for (j=1;j<=nterm-1;j++){ Update saved quantities by van Wijn- dum=wksp[j+1]; gaarden's algorithm
5.1 Series and Their Convergence 167 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). number of terms directly, through term n − 1 say, then apply the transformation to the rest of the series beginning with term n. The formula (for n even) is ∞ s=0 (−1)s us = u0 − u1 + u2 ... − un−1 +∞ s=0 (−1)s 2s+1 [∆sun] (5.1.5) Here ∆ is the forward difference operator, i.e., ∆un ≡ un+1 − un ∆2un ≡ un+2 − 2un+1 + un ∆3un ≡ un+3 − 3un+2 + 3un+1 − un etc. (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, p terms, thus requiring the first p differences (5.1.6) obtained from the terms starting at un. Euler’s transformation can be applied not only to convergent series. In some cases it will produce accurate answers from the first terms of a series that is formally divergent. It is widely used in the summation of asymptotic series. In this case it is generally wise not to sum farther than where the terms start increasing in magnitude; and you should devise some independent numerical check that the results 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 incorporation it either increases p by 1, equivalent to computing one further difference (5.1.6), or else retroactively increases n by 1, without having to redo all the difference calculations based on the old n value! The decision as to which to increase, n or p, 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: #include void eulsum(float *sum, float term, int jterm, float wksp[]) Incorporates into sum the jterm’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 call to this routine, with the first term in the series, should be with jterm=1. On the second call, term should be set to the second term of the series, with sign opposite to that of the first call, and jterm should be 2. And so on. wksp is a workspace array provided by the calling program, dimensioned at least as large as the maximum number of terms to be incorporated. { int j; static int nterm; float tmp,dum; if (jterm == 1) { Initialize: nterm=1; Number of saved differences in wksp. *sum=0.5*(wksp[1]=term); Return first estimate. } else { tmp=wksp[1]; wksp[1]=term; for (j=1;j<=nterm-1;j++) { Update saved quantities by van Wijndum=wksp[j+1]; gaarden’s algorithm
168 Chapter 5.Evaluation of Functions wksp[j+1]=0.5*(wksp[j]+tmp); tmp=dum; wksp [nterm+1]=0.5*(wksp [nterm]+tmp); if (fabs(wksp[nterm+1])<=fabs(wksp[nterm])) Favorable to increase p, *sum +(0.5*wksp[++nterm]); and the table becomes longer. else Favorable to increase n, *sum +wksp[nterm+1]; the table doesn't become longer The powerful Euler technique is not directly applicable to a series of positive terms.Occasionally it is useful to converta series of positive terms into an alternating series,just so that the Euler transformation can be used!Van Wijngaarden has given a transformation for accomplishing this [1]: -1)r-1w (5.1.7 3 where RECIPES wr三Ur+22x+4U4r+8U8r+··· (5.1.8) 2 Equations(5.1.7)and(5.1.8)replace a simple sum by a two-dimensional sum,each term in(5.1.7)being itself an infinite sum(5.1.8).This may seem a strange way to 3 Press. ART save on work!Since,however,the indices in(5.1.8)increase tremendously rapidly, as powers of 2,it often requires only a few terms to converge(5.1.8)to extraordinary Progra accuracy.You do,however,need to be able to compute the v,'s efficiently for “random”values r.The standard“updating”tricks for sequential r's,mentioned OF SCIENTIFIC( above following equation(5.1.1),can't be used. Actually,Euler's transformation is a special case of a more general transforma- 6 tion of power series.Suppose that some known function g(z)has the series ga)=∑6n2n (5.1.9) 1920 COMPUTING (ISBN n=0 and that you want to sum the new,unknown,series Numerical 10521 Recipes 43108 (5.1.10) (outside Then it is not hard to show(see [2))that equation(5.1.10)can be written as North Software. fe)=∑ae90 n! (5.1.11) n=0 which often converges much more rapidly.Here A(m)co is the nth finite-difference operator (equation 5.1.6),with A(0)coco,and g(m)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(2)=1+z 1 =1-z+z2-z3+… (5.1.12)
168 Chapter 5. Evaluation of Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). wksp[j+1]=0.5*(wksp[j]+tmp); tmp=dum; } wksp[nterm+1]=0.5*(wksp[nterm]+tmp); if (fabs(wksp[nterm+1]) <= fabs(wksp[nterm])) Favorable to increase p, *sum += (0.5*wksp[++nterm]); and the table becomes longer. else Favorable to increase n, *sum += wksp[nterm+1]; the table doesn’t become longer. } } The powerful Euler technique is not directly applicable to a series of positive terms. Occasionally it is useful to convert a series of positive terms into an alternating series, just so that the Euler transformation can be used! Van Wijngaarden has given a transformation for accomplishing this [1]: ∞ r=1 vr = ∞ 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-dimensional sum, each term in (5.1.7) being itself an infinite sum (5.1.8). This may seem a strange way to save on work! Since, however, the indices in (5.1.8) increase tremendously rapidly, as powers of 2, it often requires only a few terms to converge (5.1.8) to extraordinary 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’s transformation is a special case of a more general transformation of power series. Suppose that some known function g(z) has the series g(z) = ∞ n=0 bnzn (5.1.9) and that you want to sum the new, unknown, series f(z) = ∞ n=0 cnbnzn (5.1.10) Then it is not hard to show (see [2]) that equation (5.1.10) can be written as f(z) = ∞ n=0 [∆(n) c0] g(n) n! zn (5.1.11) which often converges much more rapidly. Here ∆(n)c0 is 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.2 Evaluation of Continued Fractions 169 into equation (5.1.11),and then setting 2=1. Sometimes you will want to compute a function from a series representation even when the computation is not efficient.For example,you may be using the values obtained to fit the function to an approximating form that you will use subsequently (cf.$5.8).If you are summing very large numbers of slowly convergent terms,pay attention to roundoff errors!In floating-point representation it is more accurate to sum a list of numbers in the order starting with the smallest one,rather than starting with the largest one.It is even better to group terms pairwise,then in pairs of pairs, etc.,so that all additions involve operands of comparable magnitude. 81 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,Numerica/Methods (Englewood Cliffs,NJ:Prentice-Hall). ICAL 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 by RECIPES Dover 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),82.3.[2] g 令 Press. 5.2 Evaluation of Continued Fractions IENTIFIC Continued fractions are often powerful ways of evaluating functions that occur in scientific applications.A continued fraction looks like this: f(z)=6o+ (5.2.1) MPUTING 02 b1+ 1920 b2+ (ISBN 44 b3+ b4++ Numerica 10.621 Printers prefer to write this as 是 43108 a1 a2 a3 a4 a5 Recipes f()=b如+1+2+b朗+1+5+ (5.2.2) (outside In either(5.2.1)or(5.2.2),the a's and b's can themselves be functions of x,usually Software. 首 linear or quadratic monomials at worst (i.e.,constants times x or times z2).For example,the continued fraction representation of the tangent function is tamr=1-3-5-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
5.2 Evaluation of Continued Fractions 169 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). into equation (5.1.11), and then setting z = 1. Sometimes you will want to compute a function from a series representation even when the computation is not efficient. For example, you may be using the values obtained to fit the function to an approximating form that you will use subsequently (cf. §5.8). If you are summing very large numbers of slowly convergent terms, pay attention to roundoff errors! In floating-point representation it is more accurate to sum a list of numbers in the order starting with the smallest one, rather than starting with the largest one. It is even better to group terms 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 Library), 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 Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 by Dover 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 Continued fractions are often powerful ways of evaluating functions that 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 x or times x2). For example, the continued fraction representation of the tangent function is tan x = 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