230 Chapter 6.Special Functions 6.5 Bessel Functions of Integer Order This section and the next one present practical algorithms for computing various kinds of Bessel functions of integer order.In $6.7 we deal with fractional order.In fact,the more complicated routines for fractional order work fine for integer order too.For integer order,however,the routines in this section (and 86.6)are simpler and faster.Their only drawback is that they are limited by the precision of the underlying rational approximations.For full double precision,it is best to work with the routines for fractional order in 86.7. For any real v,the Bessel function J(z)can be defined by the series representation 菲 J(x) (-7x2) (6.5.1) 3 The series converges for all z,but it is not computationally very useful for >1. For v not an integer the Bessel function Y()is given by Y(z)= J(z)cos(vT)-Jv() (6.5.2) sin(wT)】 The right-hand side goes to the correct limiting value Yn()as v goes to some integer n,but this is also not computationally useful. For arguments x0 Numerica 10621 For >v,both Bessel functions look qualitatively like sine or cosine waves whose amplitude decays as x-1/2.The asymptotic forms forvare 43126 2 1 (outside (6.5.4) 2. 71 1 In the transition region wherev,the typical amplitudes of the Bessel functions are on the order 21/31 0.4473 四）~3亦~ v1/3 21/3 (6.5.5) 3VT方~-0.7748 1 Y(w) v1/3
230 Chapter 6. Special Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 6.5 Bessel Functions of Integer Order This section and the next one present practical algorithms for computing various kinds of Bessel functions of integer order. In §6.7 we deal with fractional order. In fact, the more complicated routines for fractional order work fine for integer order too. For integer order, however, the routines in this section (and §6.6) are simpler and faster. Their only drawback is that they are limited by the precision of the underlying rational approximations. For full double precision, it is best to work with the routines for fractional order in §6.7. For any real ν, the Bessel function Jν(x) can be defined by the series representation Jν (x) = 1 2 x ν ∞ k=0 (−1 4x2)k k!Γ(ν + k + 1) (6.5.1) The series converges for all x, but it is not computationally very useful for x 1. For ν not an integer the Bessel function Yν (x) is given by Yν(x) = Jν(x) cos(νπ) − J−ν(x) sin(νπ) (6.5.2) The right-hand side goes to the correct limiting value Y n(x) as ν goes to some integer n, but this is also not computationally useful. For arguments x 0 (6.5.3) For x>ν, both Bessel functions look qualitatively like sine or cosine waves whose amplitude decays as x−1/2. The asymptotic forms for x ν are Jν(x) ∼ 2 πx cos x − 1 2 νπ − 1 4 π Yν(x) ∼ 2 πx sin x − 1 2 νπ − 1 4 π (6.5.4) In the transition region where x ∼ ν, the typical amplitudes of the Bessel functions are on the order Jν (ν) ∼ 21/3 32/3Γ( 2 3 ) 1 ν1/3 ∼ 0.4473 ν1/3 Yν (ν) ∼ − 21/3 31/6Γ( 2 3 ) 1 ν1/3 ∼ −0.7748 ν1/3 (6.5.5)
6.5 Bessel Functions of Integer Order 231 Permission is .com or call 1-800-872- (including this one) internet -1.5 from NUMERICAL RECIPES IN C: -20 (North 0 2 4 6 10 to any server computer, 1988-1992 by Cambridge University Press. America t users to make one paper THE Figure 6.5.1.Bessel functionsJo()through J3()and Yo(r)through Y2(). ART which holds asymptotically for large v.Figure 6.5.1 plots the first few Bessel Programs functions of each kind. The Bessel functions satisfy the recurrence relations Jn+1(回= 2”Jn()-n-1( (6.5.6) and Xa回-x(国)--回 (6.5.7) OF SCIENTIFIC COMPUTING (ISBN 0-521- 1988-1992 As already mentioned in $5.5,only the second of these (6.5.7)is stable in the direction of increasing n for x<n.The reason that(6.5.6)is unstable in the direction of increasing n is simply that it is the same recurrence as (6.5.7):A small amount of"polluting"Yn introduced by roundofferror will quickly come to swamp Fuurggoglrion the desired Jn,according to equation (6.5.3). Numerical Recipes :43198.5 A practical strategy for computing the Bessel functions of integer order divides into two tasks:first,how to compute Jo,J1,Yo,and Yi,and second,how to use the (outside recurrence relations stably to find otherJ's and Y's.We treat the first task first: North Software. For x between zero and some arbitrary value (we will use the value 8), approximate Jo()andJ1()by rational functions in z.Likewise approximate by rational functions the"regular part"of Yo(x)and Yi(r),defined as visit website 6a)-2ha回ande-回he- (6.5.8) For 8<x<o,use the approximating forms (n=0,1) a=Vr()）=x)-o.()mx】 (6.5.9)
6.5 Bessel Functions of Integer Order 231 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Bessel functions 1 .5 0 −.5 −1 −1.5 −2 0 2 4 6 8 10 Y0 Y1 Y2 J0 J1 J2 J3 x Figure 6.5.1. Bessel functions J0(x) through J3(x) and Y0(x) through Y2(x). which holds asymptotically for large ν. Figure 6.5.1 plots the first few Bessel functions of each kind. The Bessel functions satisfy the recurrence relations Jn+1(x) = 2n x Jn(x) − Jn−1(x) (6.5.6) and Yn+1(x) = 2n x Yn(x) − Yn−1(x) (6.5.7) As already mentioned in §5.5, only the second of these (6.5.7) is stable in the direction of increasing n for x<n. The reason that (6.5.6) is unstable in the direction of increasing n is simply that it is the same recurrence as (6.5.7): A small amount of “polluting” Yn introduced by roundoff error will quickly come to swamp the desired Jn, according to equation (6.5.3). A practical strategy for computing the Bessel functions of integer order divides into two tasks: first, how to compute J0, J1, Y0, and Y1, and second, how to use the recurrence relations stably to find other J’s and Y ’s. We treat the first task first: For x between zero and some arbitrary value (we will use the value 8), approximate J0(x) and J1(x) by rational functions in x. Likewise approximate by rational functions the “regular part” of Y0(x) and Y1(x), defined as Y0(x) − 2 π J0(x) ln(x) and Y1(x) − 2 π J1(x) ln(x) − 1 x (6.5.8) For 8 <x< ∞, use the approximating forms (n = 0, 1) Jn(x) = 2 πx Pn 8 x cos(Xn) − Qn 8 x sin(Xn) (6.5.9)
232 Chapter 6.Special Functions a=Vr()m(x)+o.()x 6.5.10) where 2m+1 Xn三x- (6.5.11) and where Po,P,Qo,and Q1 are each polynomials in their arguments,for 0 1-800 float bessj0(float x) Returns the Bessel function Jo(x)for any real x. float ax,z; double xx,y,ans,ans1,ans2; Accumulate polynomials in double precision. (Nor if ((ax=fabs(x)) float bessy0(float x) Software. Returns the Bessel function Yo(x)for positive x. float bessjo(float x); (outside North America) float zi double xx,y,ans,ans1,ans2; Accumulate polynomials in double precision. if(x<8.0) Rational function approximation of (6.5.8) y=x*x: ans1=-2957821389.0+y*(7062834065.0+y*(-512359803.6 +y*(10879881.29+y*(-86327.92757+y*228.4622733))); ans2=40076544269.0+y*(745249964.8+y*(7189466.438 +y*(47447.26470+y*(226.1030244+y*1.0)))); ans=(ans1/ans2)+0.636619772*bessj0(x)*1og(x); else Fitting function (6.5.10)
232 Chapter 6. Special Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Yn(x) = 2 πx Pn 8 x sin(Xn) + Qn 8 x cos(Xn) (6.5.10) where Xn ≡ x − 2n + 1 4 π (6.5.11) and where P0, P1, Q0, and Q1 are each polynomials in their arguments, for 0 float bessj0(float x) Returns the Bessel function J0(x) for any real x. { float ax,z; double xx,y,ans,ans1,ans2; Accumulate polynomials in double precision. if ((ax=fabs(x)) float bessy0(float x) Returns the Bessel function Y0(x) for positive x. { float bessj0(float x); float z; double xx,y,ans,ans1,ans2; Accumulate polynomials in double precision. if (x < 8.0) { Rational function approximation of (6.5.8). y=x*x; ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6 +y*(10879881.29+y*(-86327.92757+y*228.4622733)))); ans2=40076544269.0+y*(745249964.8+y*(7189466.438 +y*(47447.26470+y*(226.1030244+y*1.0)))); ans=(ans1/ans2)+0.636619772*bessj0(x)*log(x); } else { Fitting function (6.5.10).
6.5 Bessel Functions of Integer Order 233 z=8.0/x; y=2*2; xx=x-0.785398164; ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4 +y*(-0.2073370639e-5+y*0.2093887211a-6)); ans2=-0.1562499995e-1+y*(0.1430488765e-3 +y*(-0.6911147651e-5+y*(0.7621095161e-6 +y*(-0.934945152e-7)))): ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2); return ans; 2 http://www.n Permission is read able files #include .com or call (including this one) granted fori 19881992 float bessj1(float x) 11-800 Returns the Bessel function J(x)for any real x. float ax,z; to any double xx,y,ans,ans1,ans2; Accumulate polynomials in double precision. from NUMERICAL RECIPES IN if ((ax=fabs(x)) float bessy1(float x) Returns the Bessel function Y1(x)for positive x. float bessj1(float x); (outside North America) Software. float zi double xx,y,ans,ans1,ans2; Accumulate polynomials in double precision. if(x<8.0) Rational function approximation of(6.5.8) y=*x: ans1=x*(-0.4900604943e13+y*(0.1275274390e13 +y*(-0.5153438139e11+y*(0.7349264551e9 +y*(-0.4237922726e7+y*0.8511937935e4))); ans2=0.2499580570e14+y*(0.4244419664e12 +y*(0.3733650367e10+y*(0.2245904002e8 +y*(0.1020426050e6+y*(0.3549632885e3+y)))):
6.5 Bessel Functions of Integer Order 233 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). z=8.0/x; y=z*z; xx=x-0.785398164; ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4 +y*(-0.2073370639e-5+y*0.2093887211e-6))); ans2 = -0.1562499995e-1+y*(0.1430488765e-3 +y*(-0.6911147651e-5+y*(0.7621095161e-6 +y*(-0.934945152e-7)))); ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2); } return ans; } #include float bessj1(float x) Returns the Bessel function J1(x) for any real x. { float ax,z; double xx,y,ans,ans1,ans2; Accumulate polynomials in double precision. if ((ax=fabs(x)) float bessy1(float x) Returns the Bessel function Y1(x) for positive x. { float bessj1(float x); float z; double xx,y,ans,ans1,ans2; Accumulate polynomials in double precision. if (x < 8.0) { Rational function approximation of (6.5.8). y=x*x; ans1=x*(-0.4900604943e13+y*(0.1275274390e13 +y*(-0.5153438139e11+y*(0.7349264551e9 +y*(-0.4237922726e7+y*0.8511937935e4))))); ans2=0.2499580570e14+y*(0.4244419664e12 +y*(0.3733650367e10+y*(0.2245904002e8 +y*(0.1020426050e6+y*(0.3549632885e3+y)))));
234 Chapter 6.Special Functions ans=(ans1/ans2)+0.636619772*(bessj1(x)*1og(x)-1.0/x); else Fitting function (6.5.10). z=8.0/x; y=z*z; xx=x-2.356194491: ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4 +y*(0.2457520174e-5+y*(-0.240337019a-6))); ans2=0.04687499995+y*(-0.2002690873e-3 +y*(0.8449199096e-5+y*(-0.88228987e-6 +y*0.105787412e-6))); ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2); return ans; 83 鱼 19881992 We now turn to the second task,namely how to use the recurrence formulas (6.5.6)and(6.5.7)to get the Bessel functions n(x)and Yn(z)for n 2.The latter of these is straightforward,since its upward recurrence is always stable: from NUMERICAL RECIPESI float bessy(int n,float x) Returns the Bessel function Yn(x)for positive x and n 2. float bessy0(float x); float bessy1(float x); (North America 州bMe se to make one paper University Press. THE void nrerror(char error_text []) int ji ART float by,bym,byp,tox; if (n 2)nrerror("Index n less than 2 in bessy"); Programs tox=2.0/x; by=bessy1(x); Starting values for the recurrence. bym=bessy0(x); for (j=1;j<n;j++){ Recurrence (6.5.7). byp=j*tox*by-bym; to dir bym=by; by=byp; 1988199200 OF SCIENTIFIC COMPUTING(ISBN return by; v@cam 10-:6211 The cost of this algorithm is the call to bessy1 and bessyo(which generate a call to each of bessj1 and bessjo),plus O(n)operations in the recurrence. Numerical Recipes 43108 As for.().things are a bit more complicated.We can start the recurrence upward on n from.Jo and/1.but it will remain stable only while n does not exceed z.This is,however,just fine for calls with large z and small n,a case which (outside occurs frequently in practice. North Software. The harder case to provide for is that with x<n.The best thing to do here is to use Miller's algorithm(see discussion preceding equation 5.5.16),applying the recurrence dowmard from some arbitrary starting value and making use of the upward-unstable nature of the recurrence to put us onto the correct solution.When machine we finally arrive atJo orJ we are able to normalize the solution with the sum (5.5.16)accumulated along the way. The only subtlety is in deciding at how large an n we need start the downward recurrence so as to obtain a desired accuracy by the time we reach the n that we really want.If you play with the asymptotic forms(6.5.3)and(6.5.5),you should be able to convince yourself that the answer is to start larger than the desired n by
234 Chapter 6. Special Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). ans=(ans1/ans2)+0.636619772*(bessj1(x)*log(x)-1.0/x); } else { Fitting function (6.5.10). z=8.0/x; y=z*z; xx=x-2.356194491; ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4 +y*(0.2457520174e-5+y*(-0.240337019e-6)))); ans2=0.04687499995+y*(-0.2002690873e-3 +y*(0.8449199096e-5+y*(-0.88228987e-6 +y*0.105787412e-6))); ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2); } return ans; } We now turn to the second task, namely how to use the recurrence formulas (6.5.6) and (6.5.7) to get the Bessel functions J n(x) and Yn(x) for n ≥ 2. The latter of these is straightforward, since its upward recurrence is always stable: float bessy(int n, float x) Returns the Bessel function Yn(x) for positive x and n ≥ 2. { float bessy0(float x); float bessy1(float x); void nrerror(char error_text[]); int j; float by,bym,byp,tox; if (n < 2) nrerror("Index n less than 2 in bessy"); tox=2.0/x; by=bessy1(x); Starting values for the recurrence. bym=bessy0(x); for (j=1;j<n;j++) { Recurrence (6.5.7). byp=j*tox*by-bym; bym=by; by=byp; } return by; } The cost of this algorithm is the call to bessy1 and bessy0 (which generate a call to each of bessj1 and bessj0), plus O(n) operations in the recurrence. As for Jn(x), things are a bit more complicated. We can start the recurrence upward on n from J0 and J1, but it will remain stable only while n does not exceed x. This is, however, just fine for calls with large x and small n, a case which occurs frequently in practice. The harder case to provide for is that with x<n. The best thing to do here is to use Miller’s algorithm (see discussion preceding equation 5.5.16), applying the recurrence downward from some arbitrary starting value and making use of the upward-unstable nature of the recurrence to put us onto the correct solution. When we finally arrive at J0 or J1 we are able to normalize the solution with the sum (5.5.16) accumulated along the way. The only subtlety is in deciding at how large an n we need start the downward recurrence so as to obtain a desired accuracy by the time we reach the n that we really want. If you play with the asymptotic forms (6.5.3) and (6.5.5), you should be able to convince yourself that the answer is to start larger than the desired n by
6.5 Bessel Functions of Integer Order 235 an additive amount oforder [constant xn1/2,where the square root of the constant is,very roughly,the number of significant figures of accuracy The above considerations lead to the following function. #include #define ACC 40.0 Make larger to increase accuracy. #define BIGNO 1.0e10 http://www.nr. readable files Permission is #define BIGNI 1.0e-10 float bessj(int n,float x) Returns the Bessel function n(x)for any real x and n 2. .com or call { (including this one) granted for i float bessio(float x); float bessj1(float x); void nrerror(char error_text []) 11-800-872 int j,jsum,m; float ax,bj,bim,bip,sum,tox,ans; from NUMERICAL RECIPES IN C: if (n (float)n){ Upwards recurrence from Jo and J1. tox=2.0/ax; bim=bessi0(ax); only). bj=bessj1(ax); Programs for (j=1;j0;j--){ The downward recurrence. bjm=j*tox*bj-bjp; bjp-bj; bi=bim; idge.org 之 1988-1992 by Numerical Recipes -431085 if (fabs(bj)>BIGNO) Renormalize to prevent overflows. bj *BIGNI; bjp米=BIGNI: ans米=BIGNI; sum *BIGNI; (outside North Amer Software. if (jsum)sum +=bj; Accumulate the sum. jsum=!jsum; Change 0 to 1 or vice versa. if (j ==n)ans=bjp; Save the unnormalized answer sum=2.0+sum-bj; Compute (5.5.16) ans /sum; and use it to normalize the answer return x <0.0&(n&1)?-ans ans;
6.5 Bessel Functions of Integer Order 235 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). an additive amount of order [constant × n]1/2, where the square root of the constant is, very roughly, the number of significant figures of accuracy. The above considerations lead to the following function. #include #define ACC 40.0 Make larger to increase accuracy. #define BIGNO 1.0e10 #define BIGNI 1.0e-10 float bessj(int n, float x) Returns the Bessel function Jn(x) for any real x and n ≥ 2. { float bessj0(float x); float bessj1(float x); void nrerror(char error_text[]); int j,jsum,m; float ax,bj,bjm,bjp,sum,tox,ans; if (n (float) n) { Upwards recurrence from J0 and J1. tox=2.0/ax; bjm=bessj0(ax); bj=bessj1(ax); for (j=1;j0;j--) { The downward recurrence. bjm=j*tox*bj-bjp; bjp=bj; bj=bjm; if (fabs(bj) > BIGNO) { Renormalize to prevent overflows. bj *= BIGNI; bjp *= BIGNI; ans *= BIGNI; sum *= BIGNI; } if (jsum) sum += bj; Accumulate the sum. jsum=!jsum; Change 0 to 1 or vice versa. if (j == n) ans=bjp; Save the unnormalized answer. } sum=2.0*sum-bj; Compute (5.5.16) ans /= sum; and use it to normalize the answer. } return x < 0.0 && (n & 1) ? -ans : ans; }
236 Chapter 6.Special Functions 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 by Dover Publications,New York),Chapter 9. Hart,J.F.,et al.1968,Computer Approximations (New York:Wiley),86.8,p.141.[1] 6.6 Modified Bessel Functions of Integer Order The modified Bessel functions In(z)and Kn(x)are equivalent to the usual Bessel functionsIn and Yn evaluated for purely imaginary arguments.In detail, the relationship is In(x)=(-i)"In(ix) K国=+u(创)+.(m (6.6.1) 芝3、 9 The particular choice of prefactor and of the linear combination ofJn and Y to form K are simply choices that make the functions real-valued for real arguments z. For small arguments n,both In(z)and Kn(x)become,asymptotically, simple powers of their argument A三9 In()≈ ()” n>0 IENTIFIC 6 Ko(x)≈-n(x) (6.6.2) Kn(a)≈n- 2 2 n>0 (ISBN These expressions are virtually identical to those for Jn()and Yn(x)in this region, 10521 except for the factor of-2/m difference between Yn(z)and Kn(z).In the region Numerical xn,however,the modified functions have quite different behavior than the 431 Bessel functions, (outside Recipes In(d)≈√2 =exp(x) (6.6.3) North Kn(e)≈π =exp(-z) 2π The modified functions evidently have exponential rather than sinusoidal be- havior for large arguments(see Figure 6.6.1).The smoothness of the modified Bessel functions,once the exponential factor is removed,makes a simple polynomial approximation of a few terms quite suitable for the functions Io,I1,Ko,and K1 The following routines,based on polynomial coefficients given by Abramowitz and Stegun [1],evaluate these four functions,and will provide the basis for upward recursion for n 1 when n
236 Chapter 6. Special Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 by Dover Publications, New York), Chapter 9. Hart, J.F., et al. 1968, Computer Approximations (New York: Wiley), §6.8, p. 141. [1] 6.6 Modified Bessel Functions of Integer Order The modified Bessel functions In(x) and Kn(x) are equivalent to the usual Bessel functions Jn and Yn evaluated for purely imaginary arguments. In detail, the relationship is In(x)=(−i) nJn(ix) Kn(x) = π 2 i n+1[Jn(ix) + iYn(ix)] (6.6.1) The particular choice of prefactor and of the linear combination of J n and Yn to form Kn are simply choices that make the functions real-valued for real arguments x. For small arguments x n, both In(x) and Kn(x) become, asymptotically, simple powers of their argument In(x) ≈ 1 n! x 2 n n ≥ 0 K0(x) ≈ − ln(x) Kn(x) ≈ (n − 1)! 2 x 2 −n n > 0 (6.6.2) These expressions are virtually identical to those for Jn(x) and Yn(x) in this region, except for the factor of −2/π difference between Y n(x) and Kn(x). In the region x n, however, the modified functions have quite different behavior than the Bessel functions, In(x) ≈ 1 √ 2πx exp(x) Kn(x) ≈ π √2πx exp(−x) (6.6.3) The modified functions evidently have exponential rather than sinusoidal behavior for large arguments (see Figure 6.6.1). The smoothness of the modified Bessel functions, once the exponential factor is removed, makes a simple polynomial approximation of a few terms quite suitable for the functions I 0, I1, K0, and K1. The following routines, based on polynomial coefficients given by Abramowitz and Stegun [1], evaluate these four functions, and will provide the basis for upward recursion for n > 1 when x>n.
