f6-12

6.12HypergeometricFunctions 263Sample 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).CITED REFERENCES AND FURTHER READING: Erd´elyi, A., Magnus, W., Oberhettinger, F., and Tricomi, F.G. 1953, Higher Transcendental Functions , Vol. II, (New York: McGraw-Hill). [1] Gradshteyn, I.S., and Ryzhik, I.W. 1980, Table of Integrals, Series, and Products (New York: Academic Press). [2] Carlson, B.C. 1977, SIAM Journal on Mathematical Analysis , vol. 8, pp. 231–242. [3] Carlson,B.C.1987, MathematicsofComputation ,vol.49,pp.595–606[4];1988, op.cit.,vol.51, pp.267–280[5];1989, op.cit.,vol.53,pp.327–333[6];1991, op.cit.,vol.56,pp.267–280. [7] Bulirsch,R.1965, NumerischeMathematik ,vol.7,pp.78–90;1965, op.cit.,vol.7,pp.353–354; 1969,op. cit., vol. 13, pp. 305–315. [8] Carlson, B.C. 1979, Numerische Mathematik , vol. 33, pp. 1–16. [9] Carlson, B.C., and Notis, E.M. 1981, ACM Transactions on Mathematical Software , vol. 7, pp. 398–403. [10] Carlson, B.C. 1978, SIAM Journal on Mathematical Analysis , vol. 9, p. 524–528. [11] 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), Chapter 17. [12] Mathews, J., and Walker, R.L. 1970, Mathematical Methods of Physics , 2nd ed. (Reading, MA: W.A. Benjamin/Addison-Wesley), pp. 78–79. 6.12 Hypergeometric Functions As was discussed in §5.14,a fast, generalroutineforthe the complexhyperge- ometricfunction 2F1(a, b, c ;z),is difficultorimpossible. Thefunctionis definedas the analytic continuation of the hypergeometric series, 2F1(a, b, c ;z)=1+ab cz 1!+a(a+1 ) b(b+1 ) c(c+1 )z2 2!+··· +a(a+1 ) ...(a+j−1)b(b+1 ) ...(b+j−1) c(c+1 ) ...(c+j−1)zj j!+··· (6.12.1 ) This series converges only within the unit circle |z|<1(see[1]), but one’s interest in the function is not confined to this region. Section 5.14 discussed the method of evaluating this function by direct path integrationin the complex plane. We here merely list the routines that result. Implementation of the function hypgeois straightforward, and is described by comments in the program. The machinery associated with Chapter 16’s routine forintegratingdifferentialequations, odeint, is onlyminimallyintrusive,and need not even be completely understood: use of odeintrequires a common block with one zeroed variable, one subroutine call, and a prescribed format for the derivative routine hypdrv. The subroutine hypgeowill fail, of course, for values of ztoo close to the singularity at 1. (If you need to approach this singularity, or the one at ∞, use the “lineartransformationformulas”in §15.3of[1].) Awayfrom z=1,andformoderate values of a, b, c, it is often remarkable how few steps are required to integrate the equations. A half-dozen is typical. 264 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).FUNCTION hypgeo(a,b,c,z) COMPLEX hypgeo,a,b,c,z REAL EPS PARAMETER (EPS=1.e-6) Accuracy parameter. C USES bsstep,hypdrv,hypser,odeint Complex hypergeometric function 2F1for complex a, b, c ,a n d z, by direct integration of the hypergeometric equation in the complex plane. The branch cut is taken to lie alongthe real axis, Re z> 1. INTEGER kmax,nbad,nok EXTERNAL bsstep,hypdrv COMPLEX z0,dz,aa,bb,cc,y(2)COMMON /hypg/ aa,bb,cc,z0,dzCOMMON /path/ kmax Used by odeint . kmax=0 if (real(z)**2+aimag(z)**2.le.0.25) then Use series... call hypser(a,b,c,z,hypgeo,y(2)) return else if (real(z).lt.0.) then ...or pick a starting point for the path inte- gration. z0=cmplx(-0.5,0.) else if (real(z).le.1.0) then z0=cmplx(0.5,0.) else z0=cmplx(0.,sign(0.5,aimag(z))) endif aa=a Load the common block, used to pass pa- rameters “over the head” of odeint to hypdrv .bb=b cc=c dz=z-z0 call hypser(aa,bb,cc,z0,y(1),y(2)) Get starting function and derivative. call odeint(y,4,0.,1.,EPS,.1,.0001,nok,nbad,hypdrv,bsstep) The arguments to odeint are the vector of independent variables, its length, the starting and ending values of the dependent variable, the accuracy parameter, an initial guess for stepsize,a minimum stepsize, the (returned) number of good and bad steps taken, and the names ofthe derivative routine and the (here Bulirsch-Stoer) stepping routine. hypgeo=y(1) returnEND SUBROUTINE hypser(a,b,c,z,series,deriv) INTEGER n COMPLEX a,b,c,z,series,deriv,aa,bb,cc,fac,temp Returns the hypergeometric series 2F1and its derivative, iterating to machine accuracy. Forcabs(z) ≤1/2convergence is quite rapid. deriv=cmplx(0.,0.) fac=cmplx(1.,0.) temp=facaa=a bb=b cc=cdo 11n=1,1000 fac=((aa*bb)/cc)*fac deriv=deriv+fac fac=fac*z/nseries=temp+fac if (series.eq.temp) return temp=seriesaa=aa+1.bb=bb+1. cc=cc+1. enddo 11 pause ’convergence failure in hypser’ END 6.12HypergeometricFunctions 265Sample 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 hypdrv(s,y,dyds) REAL s COMPLEX y(2),dyds(2),aa,bb,cc,z0,dz,z Derivative subroutine for the hypergeometric equation, see text equation (5.14.4). COMMON /hypg/ aa,bb,cc,z0,dz z=z0+s*dz dyds(1)=y(2)*dzdyds(2)=((aa*bb)*y(1)-(cc-((aa+bb)+1.)*z)*y(2))*dz/(z*(1.-z))return END 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). [1]