f5-8

184 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).appropriate setting is ld=1, and the value of the derivative is the accumulated sum divided by the sampling interval h. CITED REFERENCES AND FURTHER READING: Dennis,J.E., andSchnabel,R.B. 1983, NumericalMethods forUnconstrained Optimizationand Nonlinear Equations (Englewood Cliffs, NJ: Prentice-Hall), §§5.4–5.6. [1] Ridders, C.J.F. 1982, Advances in Engineering Software , vol. 4, no. 2, pp. 75–76. [2] 5.8 Chebyshev Approximation The Chebyshev polynomial of degree nis denoted Tn(x), and is given by the explicit formula Tn(x)=c o s ( narccos x)( 5.8.1 ) This may look trigonometric at first glance (and there is in fact a close relation between the Chebyshev polynomials and the discrete Fourier transform); however (5.8.1) can be combined with trigonometric identities to yield explicit expressionsforT n(x)(see Figure 5.8.1), T0(x)=1 T1(x)=x T2(x)=2 x2−1 T3(x)=4 x3−3x T4(x)=8 x4−8x2+1 ··· Tn+1(x)=2 xT n(x)−Tn−1(x)n≥1.(5.8.2 ) (There also exist inverse formulas for the powers of xin terms of the Tn’s — see equations 5.11.2-5.11.3.) TheChebyshevpolynomialsareorthogonalintheinterval [−1,1]overaweight (1−x2)−1/2. In particular, /integraldisplay1 −1Ti(x)Tj(x)√ 1−x2dx=/braceleftBigg0 i/negationslash=j π/2 i=j/negationslash=0 πi =j=0(5.8.3 ) The polynomial Tn(x)hasnzeros in the interval [−1,1], and they are located at the points x=c o s/parenleftbiggπ(k−1 2) n/parenrightbigg k=1,2,...,n (5.8.4 ) 5.8ChebyshevApproximation 185Sample 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).Chebyshev polynomials1 .5 0 −.5 −1 −.8−.6−.4−.2 0 x.2 .4 .6 .8 1 −1T1T0 T2 T3 T6 T5 T4 Figure 5.8.1. Chebyshev polynomials T0(x)through T6(x). Note that Tjhasjroots in the interval (−1,1)and that all the polynomials are bounded between ±1. Inthis same intervalthereare n+1extrema(maximaandminima),locatedat x=c o s/parenleftbiggπk n/parenrightbigg k=0,1,...,n (5.8.5 ) At all of the maxima Tn(x)=1, while at all of the minima Tn(x)=−1; it is precisely this property that makes the Chebyshev polynomials so useful in polynomial approximation of functions. The Chebyshev polynomialssatisfy a discrete orthogonalityrelation as well as the continuous one (5.8.3): If xk(k=1,...,m )are the mzeros of Tm(x)given by (5.8.4), and if i, j < m , then m/summationdisplay k=1Ti(xk)Tj(xk)=/braceleftBigg0 i/negationslash=j m/ 2 i=j/negationslash=0 mi =j=0(5.8.6 ) It is not too dif ficult to combine equations (5.8.1), (5.8.4),and (5.8.6)to prove the followingtheorem: If f(x)is an arbitraryfunctionin the interval [−1,1], and if Ncoefficients cj,j=1,...,N, are defined by cj=2 NN/summationdisplay k=1f(xk)Tj−1(xk) =2 NN/summationdisplay k=1f/bracketleftbigg cos/parenleftbiggπ(k−1 2) N/parenrightbigg/bracketrightbigg cos/parenleftbiggπ(j−1)(k−1 2) N/parenrightbigg(5.8.7 ) 186 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).then the approximation formula f(x)≈/bracketleftBiggN/summationdisplay k=1ckTk−1(x)/bracketrightBigg −1 2c1 (5.8.8 ) isexactforxequal to all of the Nzeros of TN(x). For afixedN, equation (5.8.8) is a polynomial in xwhich approximates the function f(x)intheinterval [−1,1](whereallthezerosof TN(x)arelocated). Why isthisparticularapproximatingpolynomialbetterthananyotherone,exactonsomeother set of Npoints? The answer is notthat (5.8.8) is necessarily more accurate thansomeotherapproximatingpolynomialofthesame order N(forsomespeci fied definition of“accurate”), but rather that (5.8.8) can be truncated to a polynomial of lowerdegree m/lessmuchNinaverygracefulway,onethat doesyieldthe“mostaccurate ” approximation of degree m(in a sense that can be made precise). Suppose Nis so large that (5.8.8) is virtually a perfect approximation of f(x). Now consider the truncated approximation f(x)≈/bracketleftBigg m/summationdisplay k=1ckTk−1(x)/bracketrightBigg −1 2c1 (5.8.9 ) with the same cj’s, computed from (5.8.7). Since the Tk(x)’s are all bounded between ±1, the difference between (5.8.9) and (5.8.8) can be no larger than the sum of the neglected ck’s(k=m+1,...,N). In fact, if the ck’s are rapidly decreasing (which is the typical case), then the error is dominated by cm+1Tm(x), an oscillatory function with m+1equal extrema distributed smoothly over the interval [−1,1]. Thissmoothspreadingoutoftheerroris averyimportantproperty: The Chebyshev approximation (5.8.9) is very nearly the same polynomial as thatholygrailofapproximatingpolynomialsthe minimaxpolynomial ,which(amongall polynomialsof the same degree)has the smallest maximumdeviation fromthe true function f(x). The minimax polynomial is very dif ficult tofind; the Chebyshev approximatingpolynomialis almost identical and is very easy to compute! So, given some (perhaps dif ficult) means of computing the function f(x),w e now need algorithms for implementing (5.8.7)and (after inspection of the resulting c k’s and choiceofa truncatingvalue m) evaluating(5.8.9). Thelatter equationthen becomes an easy way of computing f(x)for all subsequent time. Thefirst of these tasks is straightforward. A generalizationof equation (5.8.7) that is here implemented is to allow the range of approximation to be between two arbitrarylimits aandb,insteadofjust −1to1. Thisiseffectedbyachangeofvariable y≡x−1 2(b+a) 1 2(b−a)(5.8.10 ) and by the approximation of f(x)by a Chebyshev polynomial in y. SUBROUTINE chebft(a,b,c,n,func) INTEGER n,NMAXREAL a,b,c(n),funcDOUBLE PRECISION PI EXTERNAL func PARAMETER (NMAX=50, PI=3.141592653589793d0) Chebyshev fit: Given a function func, lower and upper limits of the interval [ a,b], and a maximum degree n, this routine computes the ncoefficients cksuch that func (x)≈ 5.8ChebyshevApproximation 187Sample 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).[/summationtextn k=1ckTk−1(y)]−c1/2,w h e r e yand xare related by (5.8.10). This routine is to be used with moderately large n(e.g., 30 or 50), the array of c’s subsequently to be truncated at the smaller value msuch that cm+1and subsequent elements are negligible. Parameters: Maximum expected value of n,a n d π. INTEGER j,k REAL bma,bpa,fac,y,f(NMAX) DOUBLE PRECISION sumbma=0.5*(b-a)bpa=0.5*(b+a) do 11k=1,n We evaluate the function at the npoints required by (5.8.7). y=cos(PI*(k-0.5d0)/n)f(k)=func(y*bma+bpa) enddo 11 fac=2./ndo 13j=1,n sum=0.d0 We will accumulate the sum in double precision, a nicety that you can ignore. do12k=1,n sum=sum+f(k)*cos((PI*(j-1))*((k-0.5d0)/n)) enddo 12 c(j)=fac*sum enddo 13 return END (Ifyoufindthattheexecutiontimeof chebftisdominatedbythecalculationof N2cosines,ratherthanbythe Nevaluationsofyourfunction,thenyoushouldlook aheadto §12.3,especially equation12.3.22,whichshows howfast cosine transform methods can be used to evaluate equation 5.8.7.) Nowthatwe havetheChebyshevcoef ficients,howdoweevaluatetheapproxi- mation? Onecouldusetherecurrencerelationofequation(5.8.2)togeneratevalues forTk(x)from T0=1,T1=x, while also accumulating the sum of (5.8.9). It is better to use Clenshaw ’s recurrence formula ( §5.5), effecting the two processes simultaneously. Applied to the Chebyshev series (5.8.9), the recurrence is dm+2≡dm+1≡0 dj=2xd j+1−dj+2+cj j=m, m−1,..., 2 f(x)≡d0=xd2−d3+1 2c1(5.8.11 ) FUNCTION chebev(a,b,c,m,x) INTEGER m REAL chebev,a,b,x,c(m) Chebyshev evaluation: All arguments are input. c(1:m) is an array of Chebyshev coeffi- cients, the first melements of coutput from chebft(which must have been called with the same aandb). The Chebyshev polynomial/summationtextm k=1ckTk−1(y)−c1/2is evaluated at a point y=[x−(b+a)/2]/[(b−a)/2], and the result is returned as the function value. INTEGER j REAL d,dd,sv,y,y2 if ((x-a)*(x-b).gt.0.) pause ’x not in range in chebev’d=0.dd=0. y=(2.*x-a-b)/(b-a) Change of variable. y2=2.*ydo 11j=m,2,-1 Clenshaw’s recurrence. sv=d 188 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).d=y2*d-dd+c(j) dd=sv enddo 11 chebev=y*d-dd+0.5*c(1) Last step is different. return END If we are approximatingan evenfunction on the interval [−1,1], its expansion will involve only even Chebyshev polynomials. It is wasteful to call chebevwith all the odd coef ficients zero [1]. Instead, using the half-angle identity for the cosine in equation (5.8.1), we get the relation T2n(x)=Tn(2x2−1) ( 5.8.12 ) Thus we can evaluate a series of even Chebyshev polynomials by calling chebev with the evencoef ficientsstored consecutivelyin the array c, but with the argument xreplaced by 2x2−1. An odd function will have an expansion involving only odd Chebyshev poly- nomials. It is best to rewrite it as an expansion for the function f(x)/x, which involves only even Chebyshev polynomials. This will give accurate values for f(x)/xnearx=0. The coef ficients c/prime nforf(x)/xcan be found from those for f(x)by recurrence: c/prime N+1=0 c/prime n−1=2cn−c/prime n+1,n =N,N−2,...(5.8.13 ) Equation (5.8.13) follows from the recurrence relation in equation (5.8.2). IfyouinsistonevaluatinganoddChebyshevseries,theef ficientwayistoonce again use chebevwith xreplaced by y=2x2−1, and with the odd coef ficients stored consecutively in the array c. Now, however, you must also change the last formula in equation (5.8.11) to be f(x)=x[(2y−1)d2−d3+c1]( 5.8.14 ) and change the corresponding line in chebev. CITED REFERENCES AND FURTHER READING: Clenshaw,C.W.1962, MathematicalTables ,vol.5,NationalPhysicalLaboratory, (London:H.M. Stationery Office). [1] Goodwin, E.T. (ed.) 1961, Modern Computing Methods , 2nd ed. (New York: Philosophical Li- brary), Chapter 8. Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), §4.4.1, p. 104. Johnson, L.W., and Riess, R.D. 1982, Numerical Analysis , 2nd ed. (Reading, MA: Addison- Wesley), §6.5.2, p. 334. Carnahan, B., Luther, H.A., and Wilkes, J.O. 1969, Applied Numerical Methods (New York: Wiley), §1.10, p. 39. 5.9DerivativesorIntegralsofa Chebyshev-approximatedFunction 189Sample 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.9 Derivatives or Integrals of a Chebyshev-approximated Function If youhaveobtainedthe Chebyshevcoef ficients that approximatea functionin a certain range (e.g., from chebftin§5.8), then it is a simple matter to transform them to Chebyshev coef ficients corresponding to the derivative or integral of the function. Having done this, you can evaluate the derivative or integral just as if it were a function that you had Chebyshev- fittedab initio. The relevant formulas are these: If ci,i =1,...,mare the coef ficients that approximateafunction finequation(5.8.9), Ciarethecoef ficientsthatapproximate theindefiniteintegralof f,and c/prime iarethecoef ficientsthatapproximatethederivative off, then Ci=ci−1−ci+1 2(i−1)(i> 1) ( 5.9.1 ) c/prime i−1=c/prime i+1+2 (i−1)ci (i=m, m−1,..., 2) ( 5.9.2 ) Equation(5.9.1)isaugmentedbyanarbitrarychoiceof C1,correspondingtoan arbitrary constant of integration. Equation (5.9.2), which is a recurrence, is started withthevalues c/prime m=c/prime m+1=0,correspondingtonoinformationaboutthe m+1st Chebyshev coef ficient of the original function f. Here are routines for implementing equations (5.9.1) and (5.9.2). SUBROUTINE chder(a,b,c,cder,n) INTEGER nREAL a,b,c(n),cder(n) Given a,b,c(1:n) , as output from routine chebft §5.8,a n dg i v e n n, the desired degree of approximation (length of cto be used), this routine returns the array cder(1:n) ,t h e Chebyshev coefficients of the derivative of the function whose coefficients are c(1:n). INTEGER j REAL concder(n)=0. n andn-1are special cases. cder(n-1)=2*(n-1)*c(n) do 11j=n-2,1,-1 cder(j)=cder(j+2)+2*j*c(j+1) Equation (5.9.2). enddo 11 con=2./(b-a) do12j=1,n Normalize to the interval b-a. cder(j)=cder(j)*con enddo 12 returnEND SUBROUTINE chint(a,b,c,cint,n) INTEGER n REAL a,b,c(n),cint(n) Given a,b,c(1:n) , as output from routine chebft §5.8,a n dg i v e n n, the desired degree of approximation (length of cto be used), this routine returns the array cint(1:n) ,t h e Chebyshev coefficients of the integral of the function whose coefficients are c. The constant of integration is set so that the integral vanishes at a. INTEGER j REAL con,fac,sum