f4-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 4.4 ImproperIntegrals 135Sample 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).converges (with parameters as shown above) on the very first extrapolation, after just 5 calls to trapzd, while qsimprequires8calls (8 times as manyevaluationsof the integrand) and qtraprequires 13 calls (making 256 times as many evaluations of the integrand). CITED REFERENCES AND FURTHER READING: Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §§3.4–3.5. Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), §§7.4.1–7.4.2. Ralston, A., and Rabinowitz, P. 1978, A First Course in Numerical Analysis , 2nd ed. (New York: McGraw-Hill), §4.10–2. 4.4 Improper Integrals For our present purposes, an integral will be “improper” if it has any of the following problems: •its integrandgoestoafinitelimitingvalueatfiniteupperandlowerlimits, butcannotbeevaluated rightononeofthoselimits(e.g., sinx/xatx=0) •its upper limit is ∞, or its lower limit is −∞ •it has an integrablesingularity at either limit (e.g., x−1/2atx=0) •it has an integrable singularity at a known place between its upper and lower limits •it has an integrable singularity at an unknown place between its upper and lower limits If an integral is infinite (e.g.,/integraltext∞ 1x−1dx), or does not exist in a limiting sense (e.g.,/integraltext∞ −∞cos xdx), wedonotcallitimproper;wecallitimpossible. Noamountof clever algorithmics will return a meaningfulanswer to an ill-posed problem. In this section we will generalize the techniques of the preceding two sections to cover the first four problems on the above list. A more advanced discussion of quadrature with integrable singularities occurs in Chapter 18, notably §18.3. The fifth problem, singularity at unknown location, can really only be handled by theuse of a variable stepsize differential equation integration routine, as will be given in Chapter 16. We need a workhorse like the extended trapezoidal rule (equation 4.1.11), but onewhichisan openformulainthesenseof §4.1,i.e.,doesnotrequiretheintegrand tobeevaluatedattheendpoints. Equation(4.1.19),theextendedmidpointrule,isthe best choice. The reason is that (4.1.19) shares with (4.1.11) the “deep” property of havinganerrorseriesthatisentirelyevenin h. Indeedthereisaformula,notaswell knownas it oughtto be, called the SecondEuler-Maclaurinsummationformula , /integraldisplay xN x1f(x)dx =h[f3/2+f5/2+f7/2+··· +fN−3/2+fN−1/2] +B2h2 4(f/prime N−f/prime 1)+··· +B2kh2k (2k)!(1−2−2k+1)(f(2k−1) N −f(2k−1) 1 )+···(4.4.1 )