198 Chapter 5.Evaluation of Functions void pcshft(float a,float b,float d[],int n) Polynomial coefficient shift.Given a coefficient array d[..n-1],this routine generates a coefficient array[..1]such that where z and yare related by (5.8.10).i.e.,the interval-1 =j;k--) division is a miracle of high-school algebra.If you NUMERICAL d[k]-=cnst*d[k+1]; never learned it,go do so.You won't be sorry. RECIPES CITED REFERENCES AND FURTHER READING: Acton,F.S.1970,Numerica/Methods That Work;1990,corrected edition (Washington:Mathe- 号 Press. matical Association of America),pp.59.182-183 [synthetic division]. Programs 5.11 Economization of Power Series IENTIFIC One particular application of Chebyshev methods,the economization of power series,is an occasionally useful technique,with a flavor of getting something for nothing. Suppose that you are already computing a function by the use of a convergent power series,for example 9 (ISBN f)三1-+京-7+ (5.11.1) Numerica 10621 (This function is actually sin()/v,but pretend you don't know that.)You might be 43106 doing a problem that requires evaluating the series many times in some particular interval,say [0,(2)2].Everything is fine,except that the series requires a large number of terms before S Recipes its error (approximated by the first neglected term,say)is tolerable.In our example,with =(2)2,the first term smaller than 10-7 is/(27!).This then approximates the error of the finite series whose last term is 12/(25!). Software. Notice that because of the large exponent in the error is much smaller than 10-7 everywhere in the interval except at the very largest values of r.This is the feature that allows "economization":if we are willing to let the error elsewhere in the interval rise to about the same value that the first neglected term has at the extreme end of the interval,then we can replace the 13-term series by one that is significantly shorter. Here are the steps for doing so: 1.Change variables from x to y,as in equation (5.8.10),to map the z interval into -1≤y≤1. 2.Find the coefficients of the Chebyshev sum(like equation 5.8.8)that exactly equals your truncated power series (the one with enough terms for accuracy). 3.Truncate this Chebyshev series to a smaller number of terms,using the coefficient of the first neglected Chebyshev polynomial as an estimate of the error
198 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). void pcshft(float a, float b, float d[], int n) Polynomial coefficient shift. Given a coefficient array d[0..n-1], this routine generates a coefficient array g[0..n-1] such that n-1 k=0 dkyk = n-1 k=0 gkxk, where x and y are related by (5.8.10), i.e., the interval −1 =j;k--) d[k] -= cnst*d[k+1]; } CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), pp. 59, 182–183 [synthetic division]. 5.11 Economization of Power Series One particular application of Chebyshev methods, the economization of power series, is an occasionally useful technique, with a flavor of getting something for nothing. Suppose that you are already computing a function by the use of a convergent power series, for example f(x) ≡ 1 − x 3! + x2 5! − x3 7! + ··· (5.11.1) (This function is actually sin(√x)/ √x, but pretend you don’t know that.) You might be doing a problem that requires evaluating the series many times in some particular interval, say [0,(2π) 2]. Everything is fine, except that the series requires a large number of terms before its error (approximated by the first neglected term, say) is tolerable. In our example, with x = (2π) 2, the first term smaller than 10−7 is x13/(27!). This then approximates the error of the finite series whose last term is x12/(25!). Notice that because of the large exponent in x13, the error is much smaller than 10−7 everywhere in the interval except at the very largest values of x. This is the feature that allows “economization”: if we are willing to let the error elsewhere in the interval rise to about the same value that the first neglected term has at the extreme end of the interval, then we can replace the 13-term series by one that is significantly shorter. Here are the steps for doing so: 1. Change variables from x to y, as in equation (5.8.10), to map the x interval into −1 ≤ y ≤ 1. 2. Find the coefficients of the Chebyshev sum (like equation 5.8.8) that exactly equals your truncated power series (the one with enough terms for accuracy). 3. Truncate this Chebyshev series to a smaller number of terms, using the coefficient of the first neglected Chebyshev polynomial as an estimate of the error
5.11 Economization of Power Series 199 4.Convert back to a polynomial in y. 5.Change variables back to z. All of these steps can be done numerically,given the coefficients of the original power series expansion.The first step is exactly the inverse of the routine pcshft(85.10),which mapped a polynomial from y(in the interval [-1,1])to z(in the interval [a,b]).But since equation (5.8.10)is a linear relation between r and y,one can also use pcshft for the inverse.The inverse of pcshft(a,b,d,n) turns out to be (you can check this) pcshft -2-b-a2-b-a ’b-a ,d,n b-a 鱼 granted for 18881892 The second step requires the inverse operation to that done by the routine chebpc (which 1800 took Chebyshev coefficients into polynomial coefficients).The following routine,pccheb, accomplishes this,using the formula[1] from NUMERICAL RECIPESI (5.11.2) where the last term depends on whether k is even or odd America Press. k (k-1)/2 Ti(=)(k odd), + To()(k even) (5.113) IENTIFIC void pccheb(f1oatd☐，float c☐，intn) Inverse of routine chebpc:given an array of polynomial coefficients d[o..n-1],returns an equivalent array of Chebyshev coefficients c[O..n-1]. int j,jm,jp,k; float fac,pow; 1920 COMPUTING(ISBN Will be powers of 2. v@cam pow=1.0; c[0]=2.0*d[0]; 021 for (k=1;k=0;j-=2,jm--,jp++){ 膜 Increment this and lower orders of Chebyshev with the combinatorial coefficent times Software. d[k];see text for formula. c[i]+fac; fac *=((float)jm)/((float)jp); pow +pow; 2 The fourth and fifth steps are accomplished by the routines chebpc and pcshft, respectively.Here is how the procedure looks all together:
5.11 Economization of Power Series 199 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). 4. Convert back to a polynomial in y. 5. Change variables back to x. All of these steps can be done numerically, given the coefficients of the original power series expansion. The first step is exactly the inverse of the routine pcshft (§5.10), which mapped a polynomial from y (in the interval [−1, 1]) to x (in the interval [a, b]). But since equation (5.8.10) is a linear relation between x and y, one can also use pcshft for the inverse. The inverse of pcshft(a,b,d,n) turns out to be (you can check this) pcshft −2 − b − a b − a ,2 − b − a b − a ,d,n The second step requires the inverse operation to that done by the routine chebpc (which took Chebyshev coefficients into polynomial coefficients). The following routine, pccheb, accomplishes this, using the formula [1] xk = 1 2k−1 Tk(x) + k 1 Tk−2(x) + k 2 Tk−4(x) + ··· (5.11.2) where the last term depends on whether k is even or odd, ··· + k (k − 1)/2 T1(x) (k odd), ··· + 1 2 k k/2 T0(x) (k even). (5.11.3) void pccheb(float d[], float c[], int n) Inverse of routine chebpc: given an array of polynomial coefficients d[0..n-1], returns an equivalent array of Chebyshev coefficients c[0..n-1]. { int j,jm,jp,k; float fac,pow; pow=1.0; Will be powers of 2. c[0]=2.0*d[0]; for (k=1;k=0;j-=2,jm--,jp++) { Increment this and lower orders of Chebyshev with the combinatorial coefficent times d[k]; see text for formula. c[j] += fac; fac *= ((float)jm)/((float)jp); } pow += pow; } } The fourth and fifth steps are accomplished by the routines chebpc and pcshft, respectively. Here is how the procedure looks all together:
200 Chapter 5.Evaluation of Functions #define NFEW.. #define NMANY . float *c,*d,te,a,bi Economize NMANY power series coefficients e[0..NMANY-1]in the range (a,b)into NFEW coefficients d[o..NFEW-1]. c=vector(0,NMANY-1); d=vector(0,NFEW-1); e=vector(0,NMANY-1); pcshft((-2.0-b-a)/(b-a),(2.0-b-a)/(b-a),e,MANY); pccheb(e,c,NMANY); Here one would normally examine the Chebyshev coefficients c[0..NMANY-1]to decide how small NFEW can be. chebpc(c.d.NFEW): 鱼 18881892 pcshft(a,b,d,NFEW); 士/ 100 In our example,by the way,the 8th through 10th Chebyshev coefficients turn out to be on the order of-7 x 107 3 x 107',and-9 x 107,so reasonable truncations (for single precision calculations)are somewhere in this range,yielding a polynomial with 8- from NUMERICAL RECIPES I 10 terms instead of the original 13. Replacing a 13-term polynomial with a (say)10-term polynomial without any loss of (Nort server accuracy-that does seem to be getting something for nothing.Is there some magic in this technique?Not really.The 13-term polynomial defined a function f(x).Equivalent to computer, America e University Press. 令 THE economizing the series,we could instead have evaluated f()at enough points to construct its Chebyshev approximation in the interval of interest,by the methods of $5.8.We would ART have obtained just the same lower-order polynomial.The principal lesson is that the rate of convergence of Chebyshev coefficients has nothing to do with the rate of convergence of 9 Programs power series coefficients;and it is the former that dictates the number of terms needed in a polynomial approximation.A function might have a divergent power series in some region of interest,but if the function itself is well-behaved,it will have perfectly good polynomial approximations.These can be found by the methods of $5.8,but not by economization of 、三马h馆 series.There is slightly less to economization of series than meets the eye. to dir OF SCIENTIFIC COMPUTING(ISBN CITED REFERENCES AND FURTHER READING: 1988-1992 Acton,F.S.1970,Numerica/Methods That Work;1990,corrected edition (Washington:Mathe- matical Association of America),Chapter 12. Arfken,G.1970,Mathematical Methods for Physicists,2nd ed.(New York:Academic Press), p.631.[1) Numerical Recipes 10-621 43108 (outside 5.12 Pade Approximants Software. A Pade approximant,so called,is that rational function (of a specified order)whose power series expansion agrees with a given power series to the highest possible order.If the rational function is M R(x)三 k=0 N (5.12.1)
200 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). #define NFEW .. #define NMANY .. float *c,*d,*e,a,b; Economize NMANY power series coefficients e[0..NMANY-1] in the range (a, b) into NFEW coefficients d[0..NFEW-1]. c=vector(0,NMANY-1); d=vector(0,NFEW-1); e=vector(0,NMANY-1); pcshft((-2.0-b-a)/(b-a),(2.0-b-a)/(b-a),e,NMANY); pccheb(e,c,NMANY); ... Here one would normally examine the Chebyshev coefficients c[0..NMANY-1] to decide how small NFEW can be. chebpc(c,d,NFEW); pcshft(a,b,d,NFEW); In our example, by the way, the 8th through 10th Chebyshev coefficients turn out to be on the order of −7 × 10−6, 3 × 10−7, and −9 × 10−9, so reasonable truncations (for single precision calculations) are somewhere in this range, yielding a polynomial with 8 – 10 terms instead of the original 13. Replacing a 13-term polynomial with a (say) 10-term polynomial without any loss of accuracy — that does seem to be getting something for nothing. Is there some magic in this technique? Not really. The 13-term polynomial defined a function f(x). Equivalent to economizing the series, we could instead have evaluated f(x) at enough points to construct its Chebyshev approximation in the interval of interest, by the methods of §5.8. We would have obtained just the same lower-order polynomial. The principal lesson is that the rate of convergence of Chebyshev coefficients has nothing to do with the rate of convergence of power series coefficients; and it is the former that dictates the number of terms needed in a polynomial approximation. A function might have a divergent power series in some region of interest, but if the function itself is well-behaved, it will have perfectly good polynomial approximations. These can be found by the methods of §5.8, but not by economization of series. There is slightly less to economization of series than meets the eye. CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), Chapter 12. Arfken, G. 1970, Mathematical Methods for Physicists, 2nd ed. (New York: Academic Press), p. 631. [1] 5.12 Pade Approximants ´ A Pade approximant ´ , so called, is that rational function (of a specified order) whose power series expansion agrees with a given power series to the highest possible order. If the rational function is R(x) ≡ M k=0 akxk 1 + N k=1 bkxk (5.12.1)