f6-10

252 Chapter6. SpecialFunctionsSample 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).sum=sum+sign*term err=term/abs(sum) if(odd)then sign=-signsums=sum sum=sumc else sumc=sumsum=sums endif if(err.lt.EPS)goto 2odd=.not.odd enddo 12 pause ’maxits exceeded in cisi’ endif 2 si=sums ci=sumc+log(t)+EULER endifif(x.lt.0.)si=-sireturn END CITED REFERENCES AND FURTHER READING: Stegun, I.A., and Zucker, R. 1976, Journal of Research of the National Bureau of Standards , vol. 80B, pp. 291–311; 1981, op. cit., vol. 86, pp. 661–686. 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), Chapters 5 and 7. 6.10 Dawson’s Integral Dawson’s Integral F(x)is defined by F(x)=e−x2/integraldisplayx 0et2dt (6.10.1 ) The function can also be related to the complex error function by F(z)=i√π 2e−z2[1−erfc (−iz)]. (6.10.2 ) A remarkable approximation for F(z), due to Rybicki [1],i s F(z) = lim h→01√π/summationdisplay nodde−(z−nh)2 n(6.10.3 ) What makes equation (6.10.3) unusual is that its accuracy increases exponentially ashgets small, so that quite moderatevalues of h(and correspondinglyquite rapid convergence of the series) give very accurate approximations. 6.10Dawson’sIntegral 253Sample 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).We will discuss the theory that leads to equation (6.10.3) later, in §13.11, as an interesting application of Fourier methods. Here we simply implement a routinefor real values of xbased on the formula. It is first convenientto shift the summationindexto center it approximatelyon the maximum of the exponential term. Define n 0to be the even integer nearest to x/h, and x0≡n0h,x/prime≡x−x0, and n/prime≡n−n0, so that F(x)≈1√πN/summationdisplay n/prime=−N n/primeodde−(x/prime−n/primeh)2 n/prime+n0, (6.10.4 ) where the approximate equality is accurate when his sufficiently small and Nis sufficiently large. The computation of this formula can be greatly speeded up if we note that e−(x/prime−n/primeh)2=e−x/prime2e−(n/primeh)2/parenleftBig e2x/primeh/parenrightBign/prime . (6.10.5 ) The first factor is computed once, the second is an array of constants to be stored, and the third can be computed recursively, so that only two exponentials need be evaluated. Advantage is also taken of the symmetry of the coefficients e−(n/primeh)2by breakingthe summation up into positive and negativevalues of n/primeseparately. In the following routine, the choices h=0.4andN=1 1are made. Because of the symmetryof the summationsand the restriction to odd valuesof n, the limits on the doloops are 1 to 6. The accuracy of the result in this REALversion is about 2×10−7. In orderto maintain relative accuracynear x=0, where F(x)vanishes, theprogrambranchestotheevaluationofthepowerseries [2]forF(x),for|x|<0.2. FUNCTION dawson(x) INTEGER NMAX REAL dawson,x,H,A1,A2,A3 PARAMETER (NMAX=6,H=0.4,A1=2./3.,A2=0.4,A3=2./7.) Returns Dawson’s integral F(x)=e x p ( −x2)/integraltextx 0exp(t2)dtfor any real x. INTEGER i,init,n0 REAL d1,d2,e1,e2,sum,x2,xp,xx,c(NMAX) SAVE init,cDATA init/0/ Flag is 0if we need to initialize, else 1. if(init.eq.0)then init=1 do 11i=1,NMAX c(i)=exp(-((2.*float(i)-1.)*H)**2) enddo 11 endif if(abs(x).lt.0.2)then Use series expansion. x2=x**2 dawson=x*(1.-A1*x2*(1.-A2*x2*(1.-A3*x2))) else Use sampling theorem representation. xx=abs(x) n0=2*nint(0.5*xx/H) xp=xx-float(n0)*He1=exp(2.*xp*H)e2=e1**2 d1=float(n0+1) d2=d1-2.sum=0. do 12i=1,NMAX 254 Chapter6. SpecialFunctionsSample 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).sum=sum+c(i)*(e1/d1+1./(d2*e1)) d1=d1+2. d2=d2-2. e1=e2*e1 enddo 12 dawson=0.5641895835*sign(exp(-xp**2),x)*sum Constant is 1/√π. endifreturnEND Other methods for computingDawson’s integral are also known [2,3]. CITED REFERENCES AND FURTHER READING: Rybicki, G.B. 1989, Computers in Physics , vol. 3, no. 2, pp. 85–87. [1] Cody, W.J., Pociorek, K.A., and Thatcher, H.C. 1970, Mathematics of Computation , vol. 24, pp. 171–178. [2] McCabe, J.H. 1974, Mathematics of Computation , vol. 28, pp. 811–816. [3] 6.11 Elliptic Integrals and Jacobian Elliptic Functions Elliptic integralsoccurin manyapplications,becauseanyintegralof theform /integraldisplay R(t, s)dt (6.11.1 ) where Ris a rational function of tands, and sis the square root of a cubic or quartic polynomial in t, can be evaluated in terms of elliptic integrals. Standard references [1]describe how to carry out the reduction, which was originally done by Legendre. Legendre showed that only three basic elliptic integrals are required. The simplest of these is I1=/integraldisplayx ydt/radicalbig (a1+b1t)(a2+b2t)(a3+b3t)(a4+b4t)(6.11.2 ) wherewehavewrittenthequartic s2infactoredform. Instandardintegraltables [2], one of the limits of integration is always a zero of the quartic, while the other limit lies closer than the next zero, so that there is no singularity within the interval. To evaluate I1,wesimplybreaktheinterval [y,x ]intosubintervals,eachofwhicheither beginsorendsonasingularity. Thetables,therefore,needonlydistinguishtheeightcasesinwhicheachofthefourzeros(orderedaccordingtosize)appearsastheupper or lower limit of integration. In addition, when one of the b’s in (6.11.2) tends to zero, the quartic reduces to a cubic, with the largest or smallest singularity movingto±∞; this leads to eight more cases (actually just special cases of the first eight). The sixteen cases in total are then usually tabulatedin terms of Legendre’sstandard elliptic integral of the 1st kind, which we will define below. By a change of the variable of integration t, the zeros of the quartic are mapped to standard locations