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.
6.6 Modified Bessel Functions of Integer Order 237 3 Permission is 2 1K0 K 、K2 h 2 http://ww.nr.com or call 1-800-872-7423 (North America 0 0 2 3 Figure 6.6.1. Modified Bessel functions lo(x)through Is(x),Ko(x)through K2(x). only),or readable files (including this one)to any server computer,is strictly prohibited. #include float bessio(float Returns the modified Bessel function lo(x)for any real x. float ax,ans; double yi Accumulate polynomials in double precision. Copyright(C)1988-1992 by Cambridge University Press.Programs Copyright(C)1988-1992 by Numerical Recipes Sample page from NUMERICAL RECIPES IN C:THE ART OF SCIENTIFIC COMPUTING(ISBN 0-521-43108-5) if ((ax=fabs(x) float bessk0(float x) Returns the modified Bessel function Ko(x)for positive real x. float bessio(float x); double y,ans; Accumulate polynomials in double precision
6.6 Modified Bessel Functions of Integer Order 237 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). 0 1 2 3 4 01234 modified Bessel functions x K0 K1 K2 I0 I1 I2 I3 Figure 6.6.1. Modified Bessel functions I0(x) through I3(x), K0(x) through K2(x). #include float bessi0(float x) Returns the modified Bessel function I0(x) for any real x. { float ax,ans; double y; Accumulate polynomials in double precision. if ((ax=fabs(x)) float bessk0(float x) Returns the modified Bessel function K0(x) for positive real x. { float bessi0(float x); double y,ans; Accumulate polynomials in double precision
238 Chapter 6. Special Functions 1f(x float bessi1(float x) .com or call 1-800-872- internet Returns the modified Bessel function /1(x)for any real x. float ax,ans; 7423 double y; Accumulate polynomials in double precision. if ((ax=fabs(x)) To order Numerical Recipes books or personaluse.Further reproduction,or 1988-1992 by Numerical Recipes ART OF SCIENTIFIC COMPUTING(ISBN 0-521-43108-5) float bessk1(float x) Returns the modified Bessel function K1(x)for positive real x. float bessi1(float x); double y,ansi Accumulate polynomials in double precision. rsend email to directcustserv@cambridge.org(outside North America) 1f(x<=2.0)[ Polynomial fit. Software. y=x*x/4.0; ans=(1og(x/2.0)*bess11(x)+(1.0/x)*(1.0+y*(0.15443144 +y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1 +y*(-0.110404e-2+y*(-0.4686e-4))))))); else y=2.0/x: ans=(exp(-x)/sqrt(x)*(1.25331414+y*(0.23498619 +y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2 +y*(0.325614e-2+y*(-0.68245e-3))))); return ans;
238 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). if (x float bessi1(float x) Returns the modified Bessel function I1(x) for any real x. { float ax,ans; double y; Accumulate polynomials in double precision. if ((ax=fabs(x)) float bessk1(float x) Returns the modified Bessel function K1(x) for positive real x. { float bessi1(float x); double y,ans; Accumulate polynomials in double precision. if (x <= 2.0) { Polynomial fit. y=x*x/4.0; ans=(log(x/2.0)*bessi1(x))+(1.0/x)*(1.0+y*(0.15443144 +y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1 +y*(-0.110404e-2+y*(-0.4686e-4))))))); } else { y=2.0/x; ans=(exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619 +y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2 +y*(0.325614e-2+y*(-0.68245e-3))))))); } return ans; }
6.6 Modified Bessel Functions of Integer Order 239 The recurrence relation for In(x)and Kn(x)is the same as that for n(z) and Yn(z)provided that ix is substituted for x.This has the effect of changing a sign in the relation. In+1(x)= In(z)+In-1(z) (6.6.4) Kn+1()=+ 2n Kn(z)+Kn-1(x) These relations are always unstable for upward recurrence.For Kn,itself growing, this presents no problem.For In,however,the strategy of downward recursion is 83g therefore required once again,and the starting point for the recursion may be chosen 18881992 in the same manner as for the routine bessj.The only fundamental difference is 1.00 that the normalization formula for I(x)has an alternating minus sign in successive terms.which again arises from the substitution of ix for x in the formula used Z.l0e三eme previously for from NUMERICAL RECIPES I 1=Io(x)-2I2(x)+2I4(x)-2I6(x)+· (6.6.5) (Nort server to make UnN电.t 令 America computer one paper THE In fact,we prefer simply to normalize with a call to bessio. With this simple modification,the recursion routines bess j and bessy become ART the new routines bessi and bessk: Programs float bessk(int n,float x) send Returns the modified Bessel function Kn(x)for positive x and n 2. st st float bessko(float x); float bessk1(float x); void nrerror(char error_text[]); int j; float bk,bkm,bkp,tox; if (n 2)nrerror("Index n less than 2 in bessk"); tox=2.0/x; bkm=bessk0(x); Upward recurrence for all x... bk=bessk1(x); email to directcustserv@cambridge.org 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING(ISBN 0-521-43108-5) for (j=1;j machine #define ACC 40.0 Make larger to increase accuracy. #define BIGNO 1.0e10 #define BIGNI 1.0e-10 float bessi(int n,float x) Returns the modified Bessel function In(x)for any real x and n 2 float bessio(float x); void nrerror(char error-text☐)；
6.6 Modified Bessel Functions of Integer Order 239 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). The recurrence relation for In(x) and Kn(x) is the same as that for Jn(x) and Yn(x) provided that ix is substituted for x. This has the effect of changing a sign in the relation, In+1(x) = − 2n x In(x) + In−1(x) Kn+1(x)=+ 2n x Kn(x) + Kn−1(x) (6.6.4) These relations are always unstable for upward recurrence. For K n, itself growing, this presents no problem. For In, however, the strategy of downward recursion is therefore required once again, and the starting point for the recursion may be chosen in the same manner as for the routine bessj. The only fundamental difference is that the normalization formula for In(x) has an alternating minus sign in successive terms, which again arises from the substitution of ix for x in 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. With this simple modification, the recursion routines bessj and bessy become the new routines bessi and bessk: float bessk(int n, float x) Returns the modified Bessel function Kn(x) for positive x and n ≥ 2. { float bessk0(float x); float bessk1(float x); void nrerror(char error_text[]); int j; float bk,bkm,bkp,tox; if (n #define ACC 40.0 Make larger to increase accuracy. #define BIGNO 1.0e10 #define BIGNI 1.0e-10 float bessi(int n, float x) Returns the modified Bessel function In(x) for any real x and n ≥ 2. { float bessi0(float x); void nrerror(char error_text[]);
240 Chapter 6.Special Functions int ji float bi,bim,bip,tox,ans; if (n 2)nrerror("Index n less than 2 in bessi"); 1f(x==0.0) return 0.0; else tox=2.0/fabs(x); bip=ans=0.0; b1=1.0; for (j=2*(n+(int)sqrt(ACC*n));j>0;j--){Downward recurrence from even bim=bip+j*tox*bi; m. bip=bi; bi-bim; if (fabs(bi)>BIGNO){ Renormalize to prevent overflows. ans BIGNI: granted for 19881992 bi BIGNI; bip *BIGNI; 1.200 if (j =n)ans=bip; ans *bessio(x)/bi; Normalize with bessi0. to any from NUMERICAL RECIPESI return x <0.0&k(n&1)-ansans; (North America server computer, University Press. 令 one paper THE ART CITED REFERENCES AND FURTHER READING: Abramowitz,M.,and Stegun,I.A.1964,Handbook of Mathematical Functions,Applied Mathe- Programs matics Series,Volume 55(Washington:National Bureau of Standards;reprinted 1968 by Dover Publications,New York),89.8.[1] OF SCIENTIFIC Carrier,G.F.,Krook,M.and Pearson,C.E.1966,Functions of a Complex Variable (New York: McGraw-Hill),pp.220ff. 6.7 Bessel Functions of Fractional Order,Airy 1920 COMPUTING (ISBN Functions,Spherical Bessel Functions 10521 Many algorithms have been proposed for computing Bessel functions of fractional order numerically.Most of them are,in fact,not very good in practice.The routines given here are Numerical Recipes 43106 rather complicated,but they can be recommended wholeheartedly. Ordinary Bessel Functions (outside The basic idea is Steed's method,which was originally developed [1]for Coulomb wave North Software. functions.The method calculates,Y,and Y simultaneously,and so involves four relations among these functions.Three of the relations come from two continued fractions, one of which is complex.The fourth is provided by the Wronskian relation W三Y-上=名 (6.7.1) I The first continued fraction,CFl,is defined by (6.7.2) 王-20u+1)/x=2u+2/x=
240 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). int j; float bi,bim,bip,tox,ans; if (n 0;j--) { Downward recurrence from even bim=bip+j*tox*bi; m. bip=bi; bi=bim; if (fabs(bi) > BIGNO) { Renormalize to prevent overflows. ans *= BIGNI; bi *= BIGNI; bip *= BIGNI; } if (j == n) ans=bip; } ans *= bessi0(x)/bi; Normalize with bessi0. return x < 0.0 && (n & 1) ? -ans : ans; } } 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), §9.8. [1] Carrier, G.F., Krook, M. and Pearson, C.E. 1966, Functions of a Complex Variable (New York: McGraw-Hill), pp. 220ff. 6.7 Bessel Functions of Fractional Order, Airy Functions, Spherical Bessel Functions Many algorithms have been proposed for computing Bessel functions of fractional order numerically. Most of them are, in fact, not very good in practice. The routines given here are rather complicated, but they can be recommended wholeheartedly. Ordinary Bessel Functions The basic idea is Steed’s method, which was originally developed [1] for Coulomb wave functions. The method calculates Jν , J ν , Yν, and Y ν simultaneously, and so involves four relations among these functions. Three of the relations come from two continued fractions, one of which is complex. The fourth is provided by the Wronskian relation W ≡ Jν Y ν − YνJ ν = 2 πx (6.7.1) The first continued fraction, CF1, is defined by fν ≡ J ν Jν = ν x − Jν+1 Jν = ν x − 1 2(ν + 1)/x − 1 2(ν + 2)/x − ··· (6.7.2)