f5-9

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 Chebyshevcoefficients that approximatea functionin a certain range (e.g., from chebftin§5.8), then it is a simple matter to transform them to Chebyshev coefficients 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-fitted ab initio. The relevant formulas are these: If ci,i =1 ,...,mare the coefficients that approximateafunction finequation(5.8.9), Ciarethecoefficientsthatapproximate theindefiniteintegralof f,and c/prime iarethecoefficientsthatapproximatethederivative 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 coefficient 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 190 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).con=0.25*(b-a) Factor that normalizes to the interval b-a. sum=0. Accumulates the constant of integration. fac=1. Will equal ±1. do11j=2,n-1 cint(j)=con*(c(j-1)-c(j+1))/(j-1) Equation (5.9.1). sum=sum+fac*cint(j) fac=-fac enddo 11 cint(n)=con*c(n-1)/(n-1) Special case of (5.9.1) for n. sum=sum+fac*cint(n) cint(1)=2.*sum Set the constant of integration. returnEND Clenshaw-CurtisQuadrature Since a smooth function’s Chebyshev coefficients cidecrease rapidly, generally expo- nentially, equation (5.9.1) is often quite efficient as the basis for a quadrature scheme. Theroutines chebftandchint, used in that order, can be followed by repeated calls to chebev if/integraltext x af(x)dxis required for many different values of xin the range a≤x≤b. If only the single definite integral/integraltextb af(x)dxis required, then chintandchebevare replaced by the simpler formula, derived from equation (5.9.1), /integraldisplayb af(x)dx=(b−a)/bracketleftbigg1 2c1−1 3c3−1 15c5−···−1 (2k+ 1)(2 k−1)c2k+1−···/bracketrightbigg (5.9.3 ) where the ci’s are as returned by chebft. The series can be truncated when c2k+1becomes negligible, and the first neglected term gives an error estimate. This scheme is known as Clenshaw-Curtis quadrature [1]. It is often combined with an adaptive choice of N, the number of Chebyshev coefficients calculated via equation (5.8.7), which is also the number of function evaluations of f(x). If a modest choice of Ndoes not give a sufficiently small c2k+1in equation (5.9.3), then a larger value is tried. In this adaptive case, it is even better to replace equation (5.8.7) by the so-called “trapezoidal” orGauss-Lobatto ( §4.5) variant, c j=2 NN/summationdisplay/prime/prime k=0f/bracketleftbigg cos/parenleftbiggπk N/parenrightbigg/bracketrightbigg cos/parenleftbiggπ(j−1)k N/parenrightbigg j=1,...,N (5.9.4 ) where (N.B.!) the two primes signify that the first and last terms in the sum are to be multiplied by 1/2.I f Nis doubled in equation (5.9.4), then half of the new function evaluationpointsareidenticaltotheoldones,allowingthepreviousfunctionevaluationstobereused. This feature, plus the analytic weights and abscissas (cosine functions in 5.9.4), giveClenshaw-Curtis quadrature an edge over high-order adaptive Gaussian quadrature (cf. §4.5), which the method otherwise resembles. Ifyourproblemforcesyoutolargevaluesof N,youshouldbeawarethatequation(5.9.4) canbeevaluated rapidly,andsimultaneouslyforallthevaluesof j,byafastcosinetransform. (See§12.3, especially equation 12.3.17.) (We already remarked that the nontrapezoidal form (5.8.7) can also be done by fast cosine methods, cf. equation 12.3.22.) CITED REFERENCES AND FURTHER READING: Goodwin, E.T. (ed.) 1961, Modern Computing Methods , 2nd ed. (New York: Philosophical Li- brary), pp. 78–79. Clenshaw, C.W., and Curtis, A.R. 1960, Numerische Mathematik , vol. 2, pp. 197–205. [1] 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