f5-10

5.10PolynomialApproximationfromChebyshevCoefficients 191Sample 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.10 Polynomial Approximation from Chebyshev Coefficients You may well ask after reading the preceding two sections, “Must I store and evaluate my Chebyshev approximation as an array of Chebyshev coefficients for a transformedvariable y? Can’t I convertthe ck’s into actual polynomialcoefficients in the original variable xand have an approximationof the followingform?” f(x)≈m/summationdisplay k=1gkxk−1(5.10.1 ) Yes, you can do this (and we will give you the algorithm to do it), but we cautionyouagainstit: Evaluatingequation(5.10.1),wherethecoefficient g’sreflect an underlying Chebyshev approximation, usually requires more significant figures than evaluation of the Chebyshev sum directly (as by chebev). This is because the Chebyshev polynomials themselves exhibit a rather delicate cancellation: The leading coefficient of Tn(x), for example, is 2n−1; other coefficients of Tn(x)are evenbigger;yettheyallmanagetocombineintoapolynomialthatliesbetween ±1. Onlywhen mis no larger than 7 or 8 should you contemplate writing a Chebyshev fit as a direct polynomial, and even in those cases you should be willing to tolerate two orso significantfiguresless accuracythanthe roundofflimit of yourmachine. You get the g’s in equation (5.10.1)from the c’s output from chebft(suitably truncatedatamodestvalueof m)bycallinginsequencethefollowingtwoprocedures: SUBROUTINE chebpc(c,d,n) INTEGER n,NMAX REAL c(n),d(n) PARAMETER (NMAX=50) Maximum anticipated value of n. Chebyshev polynomialcoefficients. Givenacoefficient array c(1:n)oflength n, this routine generates a coefficient array d(1:n)such that/summationtextn k=1dkyk−1=/summationtextn k=1ckTk−1(y)−c1/2. The method is Clenshaw’s recurrence (5.8.11), but now applied algebraically rather than arithmetically. INTEGER j,kREAL sv,dd(NMAX) do 11j=1,n d(j)=0.dd(j)=0. enddo 11 d(1)=c(n) do13j=n-1,2,-1 do12k=n-j+1,2,-1 sv=d(k) d(k)=2.*d(k-1)-dd(k)dd(k)=sv enddo 12 sv=d(1) d(1)=-dd(1)+c(j)dd(1)=sv enddo 13 do14j=n,2,-1 d(j)=d(j-1)-dd(j) enddo 14 d(1)=-dd(1)+0.5*c(1) returnEND 192 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).SUBROUTINE pcshft(a,b,d,n) INTEGER n REAL a,b,d(n) Polynomial coefficient shift. Given a coefficient array d(1:n), this routine generates a coefficient array g(1:n)such that/summationtextn k=1dkyk−1=/summationtextn k=1gkxk−1,w h e r e xand yare related by (5.8.10), i.e., the interval −1<y< 1is mapped to the interval a<x< b. The array gis returned in d. INTEGER j,kREAL const,fac const=2./(b-a) fac=constdo 11j=2,n First we rescale by the factor const... d(j)=d(j)*fac fac=fac*const enddo 11 const=0.5*(a+b) ...which is then redefined as the desired shift. do13j=1,n-1 We accomplish the shift by synthetic division. Synthetic division is a miracle of high-school algebra. If younever learned it, go do so. You won’t be sorry.do 12k=n-1,j,-1 d(k)=d(k)-const*d(k+1) enddo 12 enddo 13 returnEND CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), pp. 59, 182–183 [synthetic division]. 5.11 Economization of Power Series One particularapplication ofChebyshev methods, the economization ofpower series ,i s 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 doingaproblemthatrequiresevaluating theseriesmanytimesinsomeparticularinterval,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, withx=( 2 π) 2, the first term smaller than 10−7isx13/(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 intheintervalexceptattheverylargestvaluesof x. Thisisthefeaturethatallows “economization”: if we are willing to let the error elsewhere in the interval rise to about thesame value that the first neglected term has at the extreme end of the interval, then we canreplace the 13-term series by one that is significantly shorter. Here are the steps for doing so: 1. Change variables from xtoy, as in equation (5.8.10), to map the xinterval into −1≤y≤1. 2. FindthecoefficientsoftheChebyshev sum(likeequation5.8.8)thatexactlyequalsyour truncated power series (the one with enough terms for accuracy).