Home / Math and Physics Files / Math / Scheid and numerical / Numerical Recipes in Fortran

f6-6

PDF · 5 pages · 59.0 KB

Excerpt from the book Numerical Recipes in Fortran 77 (Cambridge University Press, 1986-1992), pages 229-233, covering the end of section 6.5 and all of section 6.6. It gives the relations of I_n and K_n to J_n and Y_n, small and large argument asymptotics, and polynomial-fit routines bessi0, bessk0, bessi1, bessk1. It also covers the recurrence relations, with upward recurrence for K_n and downward for I_n, in routines bessk and bessi.

AI-written summary; may contain errors.

Extracted text (machine-read; may contain errors)

6.6ModiﬁedBesselFunctionsofIntegerOrder 229Sample 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).bjp=bjp*BIGNI bessj=bessj*BIGNI sum=sum*BIGNI endifif(jsum.ne.0)sum=sum+bj Accumulate the sum. jsum=1-jsum C h a n g e0t o1o rv i c ev e r s a . if(j.eq.n)bessj=bjp Save the unnormalized answer. enddo 12 sum=2.*sum-bj Compute (5.5.16) bessj=bessj/sum and use it to normalize the answer. endifif(x.lt.0..and.mod(n,2).eq.1)bessj=-bessjreturn 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), Chapter 9. Hart, J.F., et al. 1968, Computer Approximations (New York: Wiley), §6.8, p. 141. [1] 6.6 ModiﬁedBessel Functions of Integer Order The modiﬁed Bessel functions In(x)andKn(x)are equivalent to the usual Bessel functions JnandYnevaluated for purely imaginary arguments. In detail, the relationship is In(x)=(−i)nJn(ix) Kn(x)=π 2in+1[Jn(ix)+iY n(ix)](6.6.1 ) Theparticularchoiceofprefactorandofthelinearcombinationof JnandYntoform Knare simply choicesthat makethe functionsreal-valuedforreal arguments x. For small arguments x/lessmuchn, both In(x)andKn(x)become, asymptotically, simple powers of their argument In(x)≈1 n!/parenleftBigx 2/parenrightBign n≥0 K0(x)≈− ln(x) Kn(x)≈(n−1)! 2/parenleftBigx 2/parenrightBig−n n> 0(6.6.2 ) Theseexpressionsarevirtuallyidenticaltothosefor Jn(x)andYn(x)inthisregion, except for the factor of −2/πdifference between Yn(x)andKn(x). In the region 230 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).01234 01234modified Bessel functions xK0K1K2I0 I1 I2 I3 Figure 6.6.1. Modi ﬁed Bessel functions I0(x)through I3(x),K0(x)through K2(x). x/greatermuchn, however, the modi ﬁed functions have quite different behavior than the Bessel functions, In(x)≈1√ 2πxexp(x) Kn(x)≈π√ 2πxexp(−x)(6.6.3 ) The modi ﬁed functions evidently have exponential rather than sinusoidal be- havior for large arguments (see Figure 6.6.1). The smoothness of the modi ﬁed Besselfunctions,oncetheexponentialfactorisremoved,makesasimplepolynomial approximation of a few terms quite suitable for the functions I0,I1,K0, and K1. The following routines,based on polynomialcoef ﬁcients givenby Abramowitz and Stegun[1], evaluate these four functions, and will provide the basis for upward recursion for n> 1when x>n. FUNCTION bessi0(x) REAL bessi0,x Returns the modiﬁed Bessel function I0(x)for any real x. REAL ax DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7, * q8,q9,y Accumulate polynomials in double precision. SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,3.5156229d0,3.0899424d0,1.2067492d0, * 0.2659732d0,0.360768d-1,0.45813d-2/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,0.1328592d-1, * 0.225319d-2,-0.157565d-2,0.916281d-2,-0.2057706d-1, * 0.2635537d-1,-0.1647633d-1,0.392377d-2/ 6.6ModiﬁedBesselFunctionsofIntegerOrder 231Sample 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 (abs(x).lt.3.75) then y=(x/3.75)**2 bessi0=p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7))))) else ax=abs(x) y=3.75/ax bessi0=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4 * +y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))) endif return END FUNCTION bessk0(x) REAL bessk0,x C USES bessi0 Returns the modiﬁed Bessel function K0(x)for positive real x. REAL bessi0DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1, * q2,q3,q4,q5,q6,q7,y Accumulate polynomials in double precision. SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7DATA p1,p2,p3,p4,p5,p6,p7/-0.57721566d0,0.42278420d0,0.23069756d0, * 0.3488590d-1,0.262698d-2,0.10750d-3,0.74d-5/ DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,-0.7832358d-1,0.2189568d-1, * -0.1062446d-1,0.587872d-2,-0.251540d-2,0.53208d-3/ if (x.le.2.0) then Polynomial ﬁt. y=x*x/4.0 bessk0=(-log(x/2.0)*bessi0(x))+(p1+y*(p2+y*(p3+ * y*(p4+y*(p5+y*(p6+y*p7)))))) else y=(2.0/x) bessk0=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+ * y*(q4+y*(q5+y*(q6+y*q7)))))) endif returnEND FUNCTION bessi1(x) REAL bessi1,x Returns the modiﬁed Bessel function I 1(x)for any real x. REAL axDOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7, * q8,q9,y Accumulate polynomials in double precision. SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9 DATA p1,p2,p3,p4,p5,p6,p7/0.5d0,0.87890594d0,0.51498869d0, * 0.15084934d0,0.2658733d-1,0.301532d-2,0.32411d-3/ DATA q1,q2,q3,q4,q5,q6,q7,q8,q9/0.39894228d0,-0.3988024d-1, * -0.362018d-2,0.163801d-2,-0.1031555d-1,0.2282967d-1,* -0.2895312d-1,0.1787654d-1,-0.420059d-2/ if (abs(x).lt.3.75) then Polynomial ﬁt. y=(x/3.75)**2 bessi1=x*(p1+y*(p2+y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else ax=abs(x) y=3.75/axbessi1=(exp(ax)/sqrt(ax))*(q1+y*(q2+y*(q3+y*(q4+ * y*(q5+y*(q6+y*(q7+y*(q8+y*q9)))))))) if(x.lt.0.)bessi1=-bessi1 endifreturn END 232 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 bessk1(x) REAL bessk1,x C USES bessi1 Returns the modiﬁed Bessel function K1(x)for positive real x. REAL bessi1 DOUBLE PRECISION p1,p2,p3,p4,p5,p6,p7,q1, * q2,q3,q4,q5,q6,q7,y Accumulate polynomials in double precision. SAVE p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7 DATA p1,p2,p3,p4,p5,p6,p7/1.0d0,0.15443144d0,-0.67278579d0, * -0.18156897d0,-0.1919402d-1,-0.110404d-2,-0.4686d-4/ DATA q1,q2,q3,q4,q5,q6,q7/1.25331414d0,0.23498619d0,-0.3655620d-1, * 0.1504268d-1,-0.780353d-2,0.325614d-2,-0.68245d-3/ if (x.le.2.0) then Polynomial ﬁt. y=x*x/4.0bessk1=(log(x/2.0)*bessi1(x))+(1.0/x)*(p1+y*(p2+ * y*(p3+y*(p4+y*(p5+y*(p6+y*p7)))))) else y=2.0/xbessk1=(exp(-x)/sqrt(x))*(q1+y*(q2+y*(q3+ * y*(q4+y*(q5+y*(q6+y*q7)))))) endifreturn END The recurrence relation for In(x)andKn(x)is the same as that for Jn(x) andYn(x)provided that ixis substituted for x. This has the effect of changing a sign in the relation, In+1(x)=−/parenleftbigg2n x/parenrightbigg In(x)+In−1(x) Kn+1(x)=+/parenleftbigg2n x/parenrightbigg Kn(x)+Kn−1(x)(6.6.4 ) These relations are always unstablefor upward recurrence. For K n, itself growing, this presents no problem. For In, however, the strategy of downward recursion is thereforerequiredonceagain,andthestartingpointfortherecursionmaybechosenin the same manner as for the routine bessj. The only fundamental difference is that thenormalizationformulafor I n(x)has analternatingminussign insuccessive terms, which again arises from the substitution of ixforxin the formula used previously for Jn 1=I0(x)−2I2(x)+2I4(x)−2I6(x)+··· (6.6.5 ) In fact, we prefer simply to normalize with a call to bessi0. Withthissimplemodi ﬁcation,therecursionroutines bessjandbessybecome the new routines bessiandbessk: FUNCTION bessk(n,x) INTEGER n REAL bessk,x C USES bessk0,bessk1 Returns the modiﬁed Bessel function Kn(x)for positive xandn≥2. INTEGER j REAL bk,bkm,bkp,tox,bessk0,bessk1if (n.lt.2) pause ’bad argument n in bessk’ tox=2.0/x 6.6ModiﬁedBesselFunctionsofIntegerOrder 233Sample 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).bkm=bessk0(x) Upward recurrence for all x... bk=bessk1(x) do11j=1,n-1 ...and here it is. bkp=bkm+j*tox*bkbkm=bk bk=bkp enddo 11 bessk=bkreturn END FUNCTION bessi(n,x) INTEGER n,IACC REAL bessi,x,BIGNO,BIGNIPARAMETER (IACC=40,BIGNO=1.0e10,BIGNI=1.0e-10) C USES bessi0 Returns the modiﬁed Bessel function In(x)for any real xandn≥2. INTEGER j,mREAL bi,bim,bip,tox,bessi0 if (n.lt.2) pause ’bad argument n in bessi’ if (x.eq.0.) then bessi=0. else tox=2.0/abs(x)bip=0.0bi=1.0 bessi=0. m=2*((n+int(sqrt(float(IACC*n))))) Downward recurrence from even m. do 11j=m,1,-1 Make IACC larger to increase accuracy. bim=bip+float(j)*tox*bi The downward recurrence. bip=bibi=bimif (abs(bi).gt.BIGNO) then Renormalize to prevent overﬂows. bessi=bessi*BIGNI bi=bi*BIGNIbip=bip*BIGNI endif if (j.eq.n) bessi=bip enddo 11 bessi=bessi*bessi0(x)/bi Normalize with bessi0 . if (x.lt.0..and.mod(n,2).eq.1) bessi=-bessi endifreturnEND 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), §9.8. [1] Carrier, G.F., Krook, M. and Pearson, C.E. 1966, Functions of a Complex Variable (New York: McGraw-Hill), pp. 220ff.