f4-4

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∞ −∞cosxdx), 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 ) 136 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).This equation can be derived by writing out (4.2.1) with stepsize h, then writing it out again with stepsize h/2, then subtracting the first from twice the second. It is not possible to double the number of steps in the extended midpoint rule and still have the benefit of previous function evaluations (try it!). However, it is possible to triplethe number of steps and do so. Shall we do this, or double and accept the loss? On the average, tripling does a factor√ 3of unnecessary work, since the “right” number of steps for a desired accuracy criterion may in fact fall anywhere in the logarithmic interval implied by tripling. For doubling, the factor is only√ 2, but we lose an extra factor of 2 in being unable to use all the previous evaluations. Since 1.732<2×1.414, it is better to triple. Here is the resulting routine, which is directly comparable to trapzd. SUBROUTINE midpnt(func,a,b,s,n) INTEGER nREAL a,b,s,funcEXTERNAL func This routine computes the nth stage of refinement of an extended midpoint rule. funcis 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 (2/3)×3n-1additional interior points. sshould not be modified between sequential calls. INTEGER it,jREAL ddel,del,sum,tnm,xif (n.eq.1) then s=(b-a)*func(0.5*(a+b)) else it=3**(n-2)tnm=it del=(b-a)/(3.*tnm) ddel=del+del The added points alternate in spacing between delandddel. x=a+0.5*del sum=0. do 11j=1,it sum=sum+func(x)x=x+ddel sum=sum+func(x) x=x+del enddo 11 s=(s+(b-a)*sum/tnm)/3. The new sum is combined with the old integral to give a refined integral. endif returnEND The routine midpntcan exactly replace trapzdin a driver routinelike qtrap (§4.2); one simply changes call trapzd tocall midpnt , and perhaps also decreases the parameter JMAXsince 3JMAX−1(from step tripling) is a much larger number than 2JMAX−1(step doubling). TheopenformulaimplementationanalogoustoSimpson’srule( qsimpin§4.2) substitutes midpntfortrapzdand decreases JMAXas above, but now also changes the extrapolation step to be s=(9.*st-ost)/8. 4.4 ImproperIntegrals 137Sample 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).since, when the number of steps is tripled, the error decreases to 1/9th its size, not 1/4th as with step doubling. Either the modified qtrapor the modified qsimpwill fix the first problem on the list at the beginning of this section. Yet more sophisticated is to generalize Romberg integration in like manner: SUBROUTINE qromo(func,a,b,ss,choose) INTEGER JMAX,JMAXP,K,KMREAL a,b,func,ss,EPSEXTERNAL func,choose PARAMETER (EPS=1.e-6, JMAX=14, JMAXP=JMAX+1, K=5, KM=K-1) C USES polint Romberg integration on an open interval. Returns as ssthe integral of the function func from atob, using any specified integrating subroutine choose and Romberg’s method. Normally choose will be an open formula, not evaluating the function at the endpoints. It is assumed that choose triples the number of steps on each call, and that its error series contains only even powers of the number of steps. The routines midpnt ,midinf ,midsql , midsqu , are possible choices for choose . The parameters have the same meaning as in qromb. INTEGER jREAL dss,h(JMAXP),s(JMAXP) h(1)=1. do 11j=1,JMAX call choose(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 endif s(j+1)=s(j) h(j+1)=h(j)/9. This is where the assumption of step tripling and an even error series is used. enddo 11 pause ’too many steps in qromo’ END Thedifferencesbetween qromoandqromb(§4.3)aresoslightthatitisperhaps gratuitoustolist qromoinfull. It,however,is anexcellentdriverroutineforsolving all the other problems of improper integrals in our first list (except the intractablefifth), as we shall now see. The basic trick for improper integrals is to make a change of variables to eliminate the singularity, or to map an infinite range of integration to a finite one. For example, the identity /integraldisplay b af(x)dx=/integraldisplay1/a 1/b1 t2f/parenleftbigg1 t/parenrightbigg dt ab > 0( 4.4.2 ) canbeusedwith either b→∞andapositive,orwitha→− ∞andbnegative,and works for any function which decreases towards infinity faster than 1/x2. You can make the changeof variable implied by (4.4.2)either analytically and then use (e.g.) qromoandmidpntto do the numerical evaluation, oryou can let the numerical algorithm make the change of variable for you. We prefer the lattermethod as being more transparent to the user. To implement equation (4.4.2) we simply write a modified version of midpnt, called midinf, which allows bto be infinite (or, more precisely, a very large number on your particular machine, such as1×10 30), or ato be negative and infinite. 138 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).SUBROUTINE midinf(funk,aa,bb,s,n) INTEGER n REAL aa,bb,s,funk EXTERNAL funk This routine is an exact replacement for midpnt , i.e., returns as sthenth stage of refinement of the integral of funkfrom aatobb, except that the function is evaluated at evenly spaced points in 1/xrather than in x. This allows the upper limit bbto be as large and positive as the computer allows, or the lower limit aato be as large and negative, but not both. aaandbbmust have the same sign. INTEGER it,j REAL a,b,ddel,del,sum,tnm,func,xfunc(x)=funk(1./x)/x**2 This statement function effects the change of variable. b=1./aa These two statements change the limits of integration ac- cordingly. a=1./bb if (n.eq.1) then From this point on, the routine is exactly identical to midpnt . s=(b-a)*func(0.5*(a+b)) else it=3**(n-2)tnm=itdel=(b-a)/(3.*tnm) ddel=del+del x=a+0.5*delsum=0. do 11j=1,it sum=sum+func(x)x=x+ddelsum=sum+func(x) x=x+del enddo 11 s=(s+(b-a)*sum/tnm)/3. endif returnEND If you need to integrate from a negative lower limit to positive infinity, you do this bybreakingthe integralinto two pieces at some positivevalue, forexample, call qromo(funk,-5.,2.,s1,midpnt) call qromo(funk,2.,1.e30,s2,midinf)answer=s1+s2 Where should you choose the breakpoint? At a sufficiently large positive value so that the function funkis at least beginning to approach its asymptotic decrease to zero value at infinity. The polynomial extrapolation implicit in the second call to qromodeals with a polynomial in 1/x, not in x. Todealwithanintegralthathasanintegrablepower-lawsingularityatitslower limit, one also makes a change of variable. If the integrand diverges as (x−a)−γ, 0≤γ< 1, near x=a, use the identity /integraldisplayb af(x)dx=1 1−γ/integraldisplay(b−a)1−γ 0tγ 1−γf(t1 1−γ+a)dt (b>a )(4.4.3 ) If the singularity is at the upper limit, use the identity /integraldisplayb af(x)dx=1 1−γ/integraldisplay(b−a)1−γ 0tγ 1−γf(b−t1 1−γ)dt (b>a )(4.4.4 ) 4.4 ImproperIntegrals 139Sample 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).If there is a singularity at both limits, divide the integral at an interior breakpoint as in the example above. Equations (4.4.3) and (4.4.4) are particularly simple in the case of inverse square-root singularities, a case that occurs frequently in practice: /integraldisplayb af(x)dx=/integraldisplay√ b−a 02tf(a+t2)dt (b>a )( 4.4.5 ) for a singularity at a, and /integraldisplayb af(x)dx=/integraldisplay√ b−a 02tf(b−t2)dt (b>a )( 4.4.6 ) for a singularity at b. Once again, we can implement these changes of variable transparentlyto the user by defining substitute routines for midpntwhich make the change of variable automatically: SUBROUTINE midsql(funk,aa,bb,s,n) INTEGER nREAL aa,bb,s,funk EXTERNAL funk This routine is an exact replacement for midpnt , except that it allows for an inverse square- root singularity in the integrand at the lower limit aa. INTEGER it,j REAL ddel,del,sum,tnm,x,func,a,b func(x)=2.*x*funk(aa+x**2)b=sqrt(bb-aa) a=0. if (n.eq.1) then The rest of the routine is exactly like midpnt and is omitted. Similarly, SUBROUTINE midsqu(funk,aa,bb,s,n) INTEGER n REAL aa,bb,s,funkEXTERNAL funk This routine is an exact replacement for midpnt , except that it allows for an inverse square- root singularity in the integrand at the upper limit bb. INTEGER it,jREAL ddel,del,sum,tnm,x,func,a,b func(x)=2.*x*funk(bb-x**2) b=sqrt(bb-aa)a=0.if (n.eq.1) then The rest of the routine is exactly like midpnt and is omitted. One last example should suffice to show how these formulas are derived in general. Supposethe upperlimit ofintegrationis infinite,andthe integrandfalls off exponentially. Thenwewantachangeofvariablethatmaps e−xdxinto (±)dt(with the sign chosen to keep the upper limit of the new variable larger than the lower limit). Doing the integration gives by inspection t=e−xor x=−logt (4.4.7 ) 140 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).so that /integraldisplayx=∞ x=af(x)dx=/integraldisplayt=e−a t=0f(−logt)dt t(4.4.8 ) The user-transparent implementation would be SUBROUTINE midexp(funk,aa,bb,s,n) INTEGER n REAL aa,bb,s,funk EXTERNAL funk This routine is an exact replacement for midpnt , except that bbis assumed to be infinite (value passed not actually used). It is assumed that the function funkdecreases exponen- tially rapidly at infinity. INTEGER it,jREAL ddel,del,sum,tnm,x,func,a,b func(x)=funk(-log(x))/x b=exp(-aa)a=0.if (n.eq.1) then The rest of the routine is exactly like midpnt and is omitted. CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), Chapter 4. Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), §7.4.3, p. 294. Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §3.7, p. 152. 4.5 Gaussian Quadratures and Orthogonal Polynomials Intheformulasof §4.1,theintegralofafunctionwas approximatedbythesum of its functional values at a set of equally spaced points, multiplied by certain aptlychosen weighting coefficients. We saw that as we allowed ourselves more freedom in choosing the coefficients, we could achieve integration formulas of higher and higher order. The idea of Gaussian quadratures is to give ourselves the freedom to choose not only the weighting coefficients, but also the location of the abscissas at whichthe functionis to be evaluated: Theywill no longerbe equallyspaced. Thus,wewill have twicethenumberofdegreesoffreedomat ourdisposal;it will turnout that we can achieveGaussian quadratureformulaswhose orderis, essentially, twice that of the Newton-Cotesformulawith the same numberof functionevaluations. Does this sound too good to be true? Well, in a sense it is. The catch is a familiar one, which cannot be overemphasized: High order is not the same as high accuracy. High order translates to high accuracy only when the integrand is very smooth, in the sense of being “well-approximated by a polynomial.”