f4-2

130 Chapter4. Integrationof FunctionsSample 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).N = 1 234 (total after N = 4) Figure 4.2.1. Sequential calls to the routine trapzdincorporate the information fromprevious calls and evaluate the integrand only at those new points necessary to refine the grid. The bottom line shows thetotality of function evaluations after the fourth call. The routine qsimp, by weighting the intermediate results, transforms the trapezoid rule into Simpson’s rule with essentially no additional overhead. There are also formulas of higher order for this situation, but we will refrain from giving them. Thesemi-openformulas arejusttheobviouscombinationsofequations(4.1.11)– (4.1.14) with (4.1.15)–(4.1.18), respectively. At the closed end of the integration,use the weights from the former equations; at the open end use the weights from the latter equations. One example should give the idea, the formulawith error term decreasing as 1/N 3which is closed on the right and open on the left: /integraldisplayxN x1f(x)dx=h/bracketleftbigg23 12f2+7 12f3+f4+f5+ ···+fN−2+13 12fN−1+5 12fN/bracketrightbigg +O/parenleftbigg1 N3/parenrightbigg (4.1.20 ) CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), §25.4. [1] Isaacson, E., and Keller, H.B. 1966, Analysis of Numerical Methods (New York: Wiley), §7.1. 4.2 Elementary Algorithms Ourstartingpointis equation(4.1.11),theextendedtrapezoidalrule. Thereare two facts about the trapezoidal rule which make it the starting point for a variety of algorithms. One fact is rather obvious, while the second is rather “deep.” Theobviousfactisthat,forafixedfunction f(x)tobeintegratedbetweenfixed limits aand b, one can double the number of intervals in the extended trapezoidal rule without losing the benefit of previous work. The coarsest implementation of the trapezoidal rule is to average the function at its endpoints aand b. The first stage of refinementis to add to this averagethe value of the functionat the halfway point. The second stage of refinement is to add the values at the 1/4 and 3/4 points. And so on (see Figure 4.2.1). Without further ado we can write a routine with this kind of logic to it: 4.2 ElementaryAlgorithms 131Sample 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 trapzd(func,a,b,s,n) INTEGER n REAL a,b,s,func EXTERNAL func This routine computes the nth stage of refinement of an extended trapezoidal rule. func is input as the name of the function to be integrated between limits aandb, also input. When called with n=1, the routine returns as sthe crudest estimate of/integraltextb af(x)dx. Subsequent calls with n=2,3,... (in that sequential order) will improve the accuracy of sby adding 2n-2 additional interior points. sshould not be modified between sequential calls. INTEGER it,j REAL del,sum,tnm,x if (n.eq.1) then s=0.5*(b-a)*(func(a)+func(b)) else it=2**(n-2)tnm=itdel=(b-a)/tnm This is the spacing of the points to be added. x=a+0.5*del sum=0.do 11j=1,it sum=sum+func(x) x=x+del enddo 11 s=0.5*(s+(b-a)*sum/tnm) This replaces sby its refined value. endif returnEND The above routine ( trapzd) is a workhorse that can be harnessed in several ways. Thesimplestandcrudestistointegrateafunctionbytheextendedtrapezoidal rule where you know in advance (we can’t imagine how!) the number of steps youwant. If you want 2 M+1, you can accomplish this by the fragment do11j=1,m+1 call trapzd(func,a,b,s,j) enddo 11 with the answer returned as s. Much better, of course, is to refine the trapezoidal rule until some specified degree of accuracy has been achieved: SUBROUTINE qtrap(func,a,b,s) INTEGER JMAXREAL a,b,func,s,EPSEXTERNAL func PARAMETER (EPS=1.e-6, JMAX=20) C USES trapzd Returns as sthe integral of the function func from atob. The parameters EPS can be set to the desired fractional accuracy and JMAX s ot h a t2t ot h ep o w e r JMAX-1 is the maximum allowed number of steps. Integration is performed by the trapezoidal rule. INTEGER jREAL olds olds=-1.e30 Any number that is unlikely to be the average of the function at its endpoints will do here. do 11j=1,JMAX call trapzd(func,a,b,s,j) if (j.gt.5) then Avoid spurious early convergence. if (abs(s-olds).lt.EPS*abs(olds).or. * (s.eq.0..and.olds.eq.0.)) return endif olds=s enddo 11 pause ’too many steps in qtrap’ END 132 Chapter4. IntegrationofFunctionsSample 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).Unsophisticated as it is, routine qtrapis in fact a fairly robust way of doing integralsoffunctionsthatarenotverysmooth. Increasedsophisticationwill usuallytranslate into a higher-order method whose efficiency will be greater only for sufficientlysmoothintegrands. qtrapisthemethodofchoice,e.g.,foranintegrand which is a functionof a variablethat is linearly interpolatedbetweenmeasureddatapoints. Besurethatyoudonotrequiretoostringentan EPS,however: If qtraptakes too many steps in trying to achieve your required accuracy, accumulated roundoff errors may start increasing, and the routine may never converge. A value 10 −6 is just on the edge of trouble for most 32-bit machines; it is achievable when the convergence is moderately rapid, but not otherwise. We come now to the “deep” fact about the extended trapezoidal rule, equation (4.1.11). It is this: The error of the approximation, which begins with a term oforder 1/N 2,i si nf a c t entirely even whenexpressedin powersof 1/N. This follows directly from the Euler-Maclaurin Summation Formula , /integraldisplayxN x1f(x)dx=h/bracketleftbigg1 2f1+f2+f3+···+fN−1+1 2fN/bracketrightbigg −B2h2 2!(f/prime N−f/prime 1)−···−B2kh2k (2k)!(f(2k−1) N −f(2k−1) 1 )−···(4.2.1 ) Here B2kis aBernoulli number , defined by the generating function t et−1=∞/summationdisplay n=0Bntn n!(4.2.2 ) with the first few even values (odd values vanish except for B1=−1/2) B0=1 B2=1 6B4=−1 30B6=1 42 B8=−1 30B10=5 66B12=−691 2730(4.2.3 ) Equation (4.2.1) is not a convergent expansion, but rather only an asymptotic expansion whose error when truncated at any point is always less than twice the magnitude of the first neglected term. The reason that it is not convergent is that the Bernoulli numbers become very large, e.g., B50=495057205241079648212477525 66 The key point is that only evenpowers of hoccur in the errorseries of (4.2.1). This fact is not, in general, shared by the higher-order quadrature rules in §4.1. For example, equation (4.1.12) has an error series beginning with O(1/N3),b u t continuing with all subsequent powers of N:1/N4,1/N5, etc. Supposewe evaluate(4.1.11)with Nsteps, gettinga result SN, andthenagain with 2Nsteps, getting a result S2N. (This is done by any two consecutive calls of 4.2 ElementaryAlgorithms 133Sample 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).trapzd.) The leadingerror term in the second evaluationwill be 1/4 the size of the error in the first evaluation. Therefore the combination S=4 3S2N−1 3SN (4.2.4 ) willcancelouttheleadingordererrorterm. Butthere isnoerrortermoforder 1/N3, by(4.2.1). Thesurvivingerrorisoforder 1/N4,thesameasSimpson’srule. Infact, it should not take long for you to see that (4.2.4)is exactlySimpson’s rule (4.1.13), alternating2/3’s,4/3’s,andall. Thisis thepreferredmethodforevaluatingthatrule, and we can write it as a routine exactly analogous to qtrapabove: SUBROUTINE qsimp(func,a,b,s) INTEGER JMAXREAL a,b,func,s,EPS EXTERNAL func PARAMETER (EPS=1.e-6, JMAX=20) C USES trapzd Returns as sthe integral of the function func from atob. The parameters EPS can be set to the desired fractional accuracy and JMAX s ot h a t2t ot h ep o w e r JMAX-1 is the maximum allowed number of steps. Integration is performed by Simpson’s rule. INTEGER j REAL os,ost,st ost=-1.e30os= -1.e30do 11j=1,JMAX call trapzd(func,a,b,st,j) s=(4.*st-ost)/3. Compare equation (4.2.4), above. if (j.gt.5) then Avoid spurious early convergence. if (abs(s-os).lt.EPS*abs(os).or. * (s.eq.0..and.os.eq.0.)) return endifos=s ost=st enddo 11 pause ’too many steps in qsimp’END The routine qsimpwill in general be more efficient than qtrap(i.e., require fewer function evaluations) when the function to be integrated has a finite 4th derivative (i.e., a continuous 3rd derivative). The combination of qsimpand its necessary workhorse trapzdis a good one for light-duty work. CITED REFERENCES AND FURTHER READING: Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §3.3. Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), §§7.4.1–7.4.2. Forsythe, G.E., Malcolm, M.A., and Moler, C.B. 1977, Computer Methods for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall), §5.3. 134 Chapter4. IntegrationofFunctionsSample 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).4.3 Romberg Integration We can view Romberg’s method as the natural generalization of the routine qsimpin the last section to integration schemes that are of higher order than Simpson’s rule. The basic idea is to use the results from ksuccessive refinements of the extended trapezoidal rule (implemented in trapzd) to remove all terms in the error series up to but not including O(1/N2k). The routine qsimpis the case ofk=2. This is one example of a very general idea that goes by the name of Richardson’sdeferred approachto thelimit : Performsome numericalalgorithmfor various values of a parameter h, and then extrapolate the result to the continuum limit h=0. Equation(4.2.4),whichsubtractsoffthe leadingerrorterm,is a specialcase of polynomial extrapolation. In the more general Romberg case, we can use Neville’s algorithm (see §3.1) to extrapolate the successive refinements to zero stepsize. Neville’salgorithmcaninfactbecodedveryconciselywithinaRombergintegration routine. For clarity of the program, however,it seems better to do the extrapolation by subroutine call to polint, already given in §3.1. SUBROUTINE qromb(func,a,b,ss) INTEGER JMAX,JMAXP,K,KMREAL a,b,func,ss,EPSEXTERNAL func PARAMETER (EPS=1.e-6, JMAX=20, JMAXP=JMAX+1, K=5, KM=K-1) C USES polint,trapzd Returns as ssthe integral of the function func from atob. Integration is performed by Romberg’s method of order 2 K, where, e.g., K=2 is Simpson’s rule. Parameters: EPS is the fractional accuracy desired, as determined by the extrapolation error estimate; JMAX limits the total number of steps; Kis the number of points used in the extrapolation. INTEGER j REAL dss,h(JMAXP),s(JMAXP) These store the successive trapezoidal approximations and their relative stepsizes. h(1)=1. do11j=1,JMAX call trapzd(func,a,b,s(j),j)if (j.ge.K) then call polint(h(j-KM),s(j-KM),K,0.,ss,dss) if (abs(dss).le.EPS*abs(ss)) return endifs(j+1)=s(j) h(j+1)=0.25*h(j) This is a key step: The factor is 0.25 even though the stepsize is decreased by only 0.5. This makesthe extrapolation a polynomial in h 2as allowed by equation (4.2.1), not just a polynomial in h.enddo 11 pause ’too many steps in qromb’ END The routine qromb, along with its required trapzdandpolint, is quite powerfulfor sufficientlysmooth(e.g.,analytic)integrands,integratedoverintervals whichcontainnosingularities,andwheretheendpointsarealsononsingular. qromb, in such circumstances, takes many, manyfewer function evaluations than either of the routines in §4.2. For example, the integral /integraldisplay2 0x4log( x+/radicalbig x2+1 ) dx