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)
6.7 Bessel Functions of Fractional Order 241 You can easily derive it from the three-term recurrence relation for Bessel functions:Start with equation (6.5.6)and use equation (5.5.18).Forward evaluation of the continued fraction by one of the methods of $5.2 is essentially equivalent to backward recurrence of the recurrence relation.The rate of convergence of CFI is determined by the position of the turning point p=v(+1),beyond which the Bessel functions become oscillatory.Ifp, convergence is very rapid.If,then each iteration of the continued fraction effectively increases by one until p;thereafter rapid convergence sets in.Thus the number of iterations of CFI is of order x for large x.In the routine bessjy we set the maximum allowed number of iterations to 10,000.For larger you can use the usual asymptotic expressions for Bessel functions. One can show that the sign of is the same as the sign of the denominator of CFI once it has converged. 8 The complex continued fraction CF2 is defined by J+iy +i+/22-23/22-w2 1 学影 nted for p十gq三 一2 (6.7.3) Jvtirv x2(x+)+2(x+2i)+ (We sketch the derivation of CF2 in the analogous case of modified Bessel functions in the next subsection.)This continued fraction converges rapidly forrp,while convergence fails as z0.We have to adopt a special method for small r,which we describe below.For x not too small,we can ensure that p by a stable recurrence of andJ downwards to a value=,thus yielding the ratio fu at this lower value of v.This is the stable (North direction for the recurrence relation.The initial values for the recurrence are J=arbitrary, J=fuJv, (6.7.4) 。o with the sign of the arbitrary initial value of chosen to be the sign of the denominator of CF1.Choosing the initial value of/,very small minimizes the possibility of overflow during the recurrence.The recurrence relations are 9 Programs -1=二+ 元 (6.7.5) SCIENTIFIC U-1w-1- -1= 6 Once CF2 has been evaluated at=,then with the Wronskian(6.7.1)we have enough relations to solve for all four quantities.The formulas are simplified by introducing the quantity y≡Pfu (6.7.6) Then 10621 1/2 Numerica =±(g+p-Jm (6.7.7) 务是2 uction 431086 Ju=fuJu (6.7.8) Ya=YJμ (6.7.9) =(+)》 (6.7.10) North Software. The sign of in (6.7.7)is chosen to be the same as the sign of the initial in(6.7.4) Once all four functions have been determined at the value v=u,we can find them at the original value of v.For/and J,simply scale the values in (6.7.4)by the ratio of (6.7.7)to the value found after applying the recurrence(6.7.5).The quantities Y and Y can be found by starting with the values in(6.7.9)and (6.7.10)and using the stable upwards recurrence Y+1= _2yY-Y-1 (6.7.11) together with the relation Yo=Lyo-Yuti (6.7.12)
6.7 Bessel Functions of Fractional Order 241 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). You can easily derive it from the three-term recurrence relation for Bessel functions: Start with equation (6.5.6) and use equation (5.5.18). Forward evaluation of the continued fraction by one of the methods of §5.2 is essentially equivalent to backward recurrence of the recurrence relation. The rate of convergence of CF1 is determined by the position of the turning point xtp = ν(ν + 1) ≈ ν, beyond which the Bessel functions become oscillatory. If x ∼ xtp, then each iteration of the continued fraction effectively increases ν by one until x ∼ xtp, while convergence fails as x → 0. We have to adopt a special method for small x, which we describe below. For x not too small, we can ensure that x >∼ xtp by a stable recurrence of Jν and J ν downwards to a value ν = µ <∼ x, thus yielding the ratio fµ at this lower value of ν. This is the stable direction for the recurrence relation. The initial values for the recurrence are Jν = arbitrary, J ν = fν Jν , (6.7.4) with the sign of the arbitrary initial value of Jν chosen to be the sign of the denominator of CF1. Choosing the initial value of Jν very small minimizes the possibility of overflow during the recurrence. The recurrence relations are Jν−1 = ν xJν + J ν J ν−1 = ν − 1 x Jν−1 − Jν (6.7.5) Once CF2 has been evaluated at ν = µ, then with the Wronskian (6.7.1) we have enough relations to solve for all four quantities. The formulas are simplified by introducing the quantity γ ≡ p − fµ q (6.7.6) Then Jµ = ± W q + γ(p − fµ) 1/2 (6.7.7) J µ = fµJµ (6.7.8) Yµ = γJµ (6.7.9) Y µ = Yµ p + q γ (6.7.10) The sign of Jµ in (6.7.7) is chosen to be the same as the sign of the initial Jν in (6.7.4). Once all four functions have been determined at the value ν = µ, we can find them at the original value of ν. For Jν and J ν , simply scale the values in (6.7.4) by the ratio of (6.7.7) to the value found after applying the recurrence (6.7.5). The quantities Yν and Y ν can be found by starting with the values in (6.7.9) and (6.7.10) and using the stable upwards recurrence Yν+1 = 2ν x Yν − Yν−1 (6.7.11) together with the relation Y ν = ν xYν − Yν+1 (6.7.12)
242 Chapter 6.Special Functions Now turn to the case of small when CF2 is not suitable.Temme[2]has given a good method of evaluating Yi and Y,and hence Y from(6.7.12),by series expansions that accurately handle the singularity as x-0.The expansions work only for1/2, and so now the recurrence (6.7.5)is used to evaluate f at a value=u in this interval. Then one calculates from W Jn=y-Yfu (6.7.13) andJ from (6.7.8).The values at the original value of v are determined by scaling as before, and the Y's are recurred up as before. Temme's series are 8 =-月 ckgk 1=-2cht (6.7.14) 元4 k三0 Here 94=-z/4* (6.7.15) while the coefficients g and hk are defined in terms of quantities pk,q&,and fk that can be found by recursion: RECIPES 9=+2m2(经)q ≥ hk =-kgk +Pk %-1 Press. k-v (6.7.16) 9张二张1 k+v Programs fk= kfk-1+pk-1+9k-1 k2-w2 IENTIFIC The initial values for the recurrences are 6 ) - P0= T(1+) =(}r1- MPUTING (6.7.17) 1992 6=2vm sin ohi)+()ro间 (ISBN with 喜季E Numerica 10621 Recipes 43108 1 1 1 (w)=2元1- (6.7.18) T(1+w) (outside 1 1 1 r2(0）=2a-m+T1+o] North Software. The whole point of writing the formulas in this way is that the potential problems as0 can be controlled by evaluating v/sin v,sinh o/o,and Ti carefully.In particular,Temme gives Chebyshev expansions for Ti(v)and T2().We have rearranged his expansion for Ii to be explicitly an even series in so that we can use our routine chebev as explained in 85.8. The routine assumes>0.For negative v you can use the reflection formulas J-v=coswπJw-sin vr Yv (6.7.19) Y_v=sin J+Cosv Yu The routine also assumes>0.For 0.For =0,Y is singular
242 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). Now turn to the case of small x, when CF2 is not suitable. Temme [2] has given a good method of evaluating Yν and Yν+1, and hence Y ν from (6.7.12), by series expansions that accurately handle the singularity as x → 0. The expansions work only for |ν| ≤ 1/2, and so now the recurrence (6.7.5) is used to evaluate fν at a value ν = µ in this interval. Then one calculates Jµ from Jµ = W Y µ − Yµfµ (6.7.13) and J µ from (6.7.8). The values at the original value of ν are determined by scaling as before, and the Y ’s are recurred up as before. Temme’s series are Yν = −∞ k=0 ckgk Yν+1 = − 2 x ∞ k=0 ckhk (6.7.14) Here ck = (−x2/4)k k! (6.7.15) while the coefficients gk and hk are defined in terms of quantities pk, qk, and fk that can be found by recursion: gk = fk + 2 ν sin2 νπ 2 qk hk = −kgk + pk pk = pk−1 k − ν qk = qk−1 k + ν fk = kfk−1 + pk−1 + qk−1 k2 − ν2 (6.7.16) The initial values for the recurrences are p0 = 1 π x 2 −ν Γ(1 + ν) q0 = 1 π x 2 ν Γ(1 − ν) f0 = 2 π νπ sin νπ cosh σΓ1(ν) + sinh σ σ ln 2 x Γ2(ν) (6.7.17) with σ = ν ln 2 x Γ1(ν) = 1 2ν 1 Γ(1 − ν) − 1 Γ(1 + ν) Γ2(ν) = 1 2 1 Γ(1 − ν) + 1 Γ(1 + ν) (6.7.18) The whole point of writing the formulas in this way is that the potential problems as ν → 0 can be controlled by evaluating νπ/ sin νπ, sinh σ/σ, and Γ1 carefully. In particular, Temme gives Chebyshev expansions for Γ1(ν) and Γ2(ν). We have rearranged his expansion for Γ1 to be explicitly an even series in ν so that we can use our routine chebev as explained in §5.8. The routine assumes ν ≥ 0. For negative ν you can use the reflection formulas J−ν = cos νπ Jν − sin νπ Yν Y−ν = sin νπ Jν + cos νπ Yν (6.7.19) The routine also assumes x > 0. For x 0. For x = 0, Yν is singular.
6.7 Bessel Functions of Fractional Order 243 Internal arithmetic in the routine is carried out in double precision.The complex arithmetic is carried out explicitly with real variables. #include #include "nrutil.h" #define EPS 1.0e-10 #define FPMIN 1.0e-30 #define MAXIT 10000 #define XMIN 2.0 #def1nePI3.141592653589793 http://www. void bessjy(float x,float xnu,float *rj,float *ry,float *rjp,float *ryp) Returns the Bessel functions rj =Jv.ry =Yo and their derivatives rjp=J ryp =Y,for positive x and for xnu=>0.The relative accuracy is within one or two significant digits of EPS,except near a zero of one of the functions,where EPS controls its absolute accuracy. 8 FPMIN is a number close to the machine's smallest floating-point number.All internal arithmetic 虽 granted for 18881892 is in double precision.To convert the entire routine to double precision,change the float declarations above to double and decrease EPS to 10-16.Also convert the function beschb 1800 void beschb(double x,double *gam1,double *gam2,double *gampl, double *gammi); from NUMERICAL RECIPESI int i,isign,1,nl; 可 double a,b,br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli,dr,e,f,fact,fact2, fact3,ff,gam,gami,gam2,gammi,gampl,h,p,pimu,pimu2,q,r,rjl, (Nort rjl1,rjmu,rjpi,rjpl,rjtemp,ry1,rymu,rymup,rytemp,sum,sumi, temp,日，x2,x1,x12,Xmu,Xmu2; America server computer, one paper University Press. THE if (x =1;1--){ rjtemp=fact*rjl+rjpl; fact -xi;
6.7 Bessel Functions of Fractional Order 243 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). Internal arithmetic in the routine is carried out in double precision. The complex arithmetic is carried out explicitly with real variables. #include #include "nrutil.h" #define EPS 1.0e-10 #define FPMIN 1.0e-30 #define MAXIT 10000 #define XMIN 2.0 #define PI 3.141592653589793 void bessjy(float x, float xnu, float *rj, float *ry, float *rjp, float *ryp) Returns the Bessel functions rj = Jν, ry = Yν and their derivatives rjp = J ν, ryp = Y ν , for positive x and for xnu = ν ≥ 0. The relative accuracy is within one or two significant digits of EPS, except near a zero of one of the functions, where EPS controls its absolute accuracy. FPMIN is a number close to the machine’s smallest floating-point number. All internal arithmetic is in double precision. To convert the entire routine to double precision, change the float declarations above to double and decrease EPS to 10−16. Also convert the function beschb. { void beschb(double x, double *gam1, double *gam2, double *gampl, double *gammi); int i,isign,l,nl; double a,b,br,bi,c,cr,ci,d,del,del1,den,di,dlr,dli,dr,e,f,fact,fact2, fact3,ff,gam,gam1,gam2,gammi,gampl,h,p,pimu,pimu2,q,r,rjl, rjl1,rjmu,rjp1,rjpl,rjtemp,ry1,rymu,rymup,rytemp,sum,sum1, temp,w,x2,xi,xi2,xmu,xmu2; if (x MAXIT) nrerror("x too large in bessjy; try asymptotic expansion"); rjl=isign*FPMIN; Initialize Jν and J ν for downward recurrence. rjpl=h*rjl; rjl1=rjl; Store values for later rescaling. rjp1=rjpl; fact=xnu*xi; for (l=nl;l>=1;l--) { rjtemp=fact*rjl+rjpl; fact -= xi;
244 Chapter 6. Special Functions rjpl=fact*rjtemp-rjl; rjl=rjtemp; if (ril==0.0)rjl=EPS; f=rjpl/ril; Now have unnormalized J andJ if (x XMIN) Use series. x2=0.5*x; pimu=PI*xmu; fact =(fabs(pimu)<EPS 1.0 pimu/sin(pimu)); d=-1og(x2); e=xmu*d; fact2 =(fabs(e)<EPS 1.0 sinh(e)/e); http://www.nr. readable files Permission is beschb(xmu,&gam1,&gam2,&gampl,&gammi);Chebyshev evaluation of T1 and T2. ff=2.0/PI*fact*(gam1*cosh(e)+gam2*fact2*d); f0- e=exp(e); 83 p=e/(gampl*PI); po (including this one) granted fori g=1.0/(e*PI*gammi); 90- pimu2=0.5*pimu; fact3 =(fabs(pimu2)<EPS 1.0:sin(pimu2)/pimu2); 11-800-872 internet r=PI*pimu2*fact3*fact3; c=1.0; d=-x2*x2; from NUMERICAL RECIPES IN sum-ff+r*q; sumi=p; for(1=1;i<=MxIT;1++) ff=(i*ff+p+q)/(i*i-xmu2) -7423(North America c*=(d/i); t users to make one paper copy for their to any server computer,is strictly prohibited. 1988-1992 by Cambridge University Press.Programs Copyright (C) THE p /(i-xmu); q/=(1+xmu); del=c*(ff+r*g); only), sum +del; del1=c*p-i*del; sum1 +del1; if (fabs(del)<(1.0+fabs(sum))*EPS)break; email if (i MAXIT) nrerror("bessy series failed to converge"); rymu =-sum; ry1=“sum1米X12: rymup=xmu*xi*rymu-ry1; rjmu=w/(rymup-f*rymu); Equation (6.7.13) ART OF SCIENTIFIC COMPUTING(ISBN 0-521 else Evaluate CF2 by modified Lentz's method($5.2). a=0.25-Xmu2; p=-0.5*x1 9=1.0; br=2.0*x; b1=2.0: v@cambridge.org 1988-1992 by Numerical Recipes fact=a*xi/(p*p+q*q); -431085 cr=br+q*fact; ci=bi+p*fact; den=br*br+bi*bi: dr=br/den; Software. di =-bi/den; dlr=cr*dr-ci*di; (outside North America) dli=cr*di+ci*dr: temp=p*dlr-q*dli; ying of machine q=p*dli+q*dlr; p=temp; for(1=2;1<=MAXIT;1++)[ a+=2*(1-1); bi+=2.0; dr=a*dr+br: disa*di+bi; if (fabs(dr)+fabs(di)<FPMIN)dr=FPMIN; fact=a/(cr*cr+ci*ci);
244 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). rjpl=fact*rjtemp-rjl; rjl=rjtemp; } if (rjl == 0.0) rjl=EPS; f=rjpl/rjl; Now have unnormalized Jµ and J µ. if (x MAXIT) nrerror("bessy series failed to converge"); rymu = -sum; ry1 = -sum1*xi2; rymup=xmu*xi*rymu-ry1; rjmu=w/(rymup-f*rymu); Equation (6.7.13). } else { Evaluate CF2 by modified Lentz’s method (§5.2). a=0.25-xmu2; p = -0.5*xi; q=1.0; br=2.0*x; bi=2.0; fact=a*xi/(p*p+q*q); cr=br+q*fact; ci=bi+p*fact; den=br*br+bi*bi; dr=br/den; di = -bi/den; dlr=cr*dr-ci*di; dli=cr*di+ci*dr; temp=p*dlr-q*dli; q=p*dli+q*dlr; p=temp; for (i=2;i<=MAXIT;i++) { a += 2*(i-1); bi += 2.0; dr=a*dr+br; di=a*di+bi; if (fabs(dr)+fabs(di) < FPMIN) dr=FPMIN; fact=a/(cr*cr+ci*ci);
6.7 Bessel Functions of Fractional Order 245 cr=br+cr*fact; ci=bi-ci*fact: if (fabs(cr)+fabs(ci)MAXIT)nrerror("cf2 failed in bessjy"); gam=(p-f)/q; Equations(6.7.6)-(6.7.10) rjmu=sqrt(w/((p-f)*gam+q)); .com or call granted fori rjmu=SIGN(rjmu,rjl); rymu=rjmu*gam; rymup=rymu*(p+q/gam); ry1=xmu*xi*rymu-rymup; 11-800-872 (including this one) fact=rjmu/rjl; *ri=ril1*fact; Scale original J and J *rjp=rjpi*fact; for (i=1;i<=nl;i++) Upward recurrence of yo rytemp=(xmu+i)*xi2*ry1-rymu 7423 (North America to any server computer,is strictly prohibited. tusers to make one paper 1988-1992 by Cambridge University Press.Programs from NUMERICAL RECIPES IN C: THE rymu=ry1; ry1=rytemp; only), *ry=rymu; *ryp=xnu*xi*rymu-ry1; copy for their #define NUSE1 5 mail to directcustsen Copyright(C) #define NUSE2 5 ART OF SCIENTIFIC COMPUTING(ISBN 0-521 void beschb(double x,double *gam1,double *gam2,double *gampl,double *gammi) Evaluates I1 and T2 by Chebyshev expansion for x<1/2.Also returns 1/T(1+x)and v@cambr 1/T(1-x).If converting to double precision,set NUSE1 7,NUSE2 =8. float chebev(float a,float b,float c[],int m,float x); To order Numerical Recipes books or 1988-1992 by Numerical Recipes float xx; -431085 static f1oatc1[▣={ -1.142022680371168e0,6.5165112670737e-3, 3.087090173086e-4,-3.4706269649e-6,6.9437664e-9 3.67795e-11,-1.356e-13}: Software. static float c2[]= 1.843740587300905e0,-7.68528408447867e-2, (outside North America) 1.2719271366546e-3,-4.9717367042e-6,-3.31261198e-8, 2.423096e-10,-1.702e-13,-1.49e-15]; ying of machine Xx=8.0*x*x-1.0; Multiply x by 2 to make range be-1 to 1, *gam1=chebev(-1.0,1.0,c1,NUSE1,xx); and then apply transformation for eval- *gam2=chebev(-1.0,1.0,c2,NUSE2,xx) uating even Chebyshev series. *gampl=*gam2-x*(*gam1) *gammi=*gam2+x*(*gam1);
6.7 Bessel Functions of Fractional Order 245 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). cr=br+cr*fact; ci=bi-ci*fact; if (fabs(cr)+fabs(ci) MAXIT) nrerror("cf2 failed in bessjy"); gam=(p-f)/q; Equations (6.7.6) – (6.7.10). rjmu=sqrt(w/((p-f)*gam+q)); rjmu=SIGN(rjmu,rjl); rymu=rjmu*gam; rymup=rymu*(p+q/gam); ry1=xmu*xi*rymu-rymup; } fact=rjmu/rjl; *rj=rjl1*fact; Scale original Jν and J ν. *rjp=rjp1*fact; for (i=1;i<=nl;i++) { Upward recurrence of Yν. rytemp=(xmu+i)*xi2*ry1-rymu; rymu=ry1; ry1=rytemp; } *ry=rymu; *ryp=xnu*xi*rymu-ry1; } #define NUSE1 5 #define NUSE2 5 void beschb(double x, double *gam1, double *gam2, double *gampl, double *gammi) Evaluates Γ1 and Γ2 by Chebyshev expansion for |x| ≤ 1/2. Also returns 1/Γ(1 + x) and 1/Γ(1 − x). If converting to double precision, set NUSE1 = 7, NUSE2 = 8. { float chebev(float a, float b, float c[], int m, float x); float xx; static float c1[] = { -1.142022680371168e0,6.5165112670737e-3, 3.087090173086e-4,-3.4706269649e-6,6.9437664e-9, 3.67795e-11,-1.356e-13}; static float c2[] = { 1.843740587300905e0,-7.68528408447867e-2, 1.2719271366546e-3,-4.9717367042e-6,-3.31261198e-8, 2.423096e-10,-1.702e-13,-1.49e-15}; xx=8.0*x*x-1.0; Multiply x by 2 to make range be −1 to 1, and then apply transformation for evaluating even Chebyshev series. *gam1=chebev(-1.0,1.0,c1,NUSE1,xx); *gam2=chebev(-1.0,1.0,c2,NUSE2,xx); *gampl= *gam2-x*(*gam1); *gammi= *gam2+x*(*gam1); }
246 Chapter 6. Special Functions Modified Bessel Functions Steed's method does not work for modified Bessel functions because in this case CF2 is purely imaginary and we have only three relations among the four functions.Temme[3]has given a normalization condition that provides the fourth relation. The Wronskian relation is W=1K-K=- (6.7.20) The continued fraction CFl becomes 1 1 fv F元=五+2u+1/x+2u+2/z+ = (6.7.21) 81 To get CF2 and the normalization condition in a convenient form,consider the sequence of confluent hypergeometric functions zn(x)=U(w+1/2+n,2v+1,2x) (6.7.22) for fixed v.Then K(z)=x/2(2x)”er20(x) (6.7.23) RECIPES -+++(e-)周 K+1(x）_1[,1 (6.7.24) 2 K(r) Equation (6.7.23)is the standard expression for K in terms of a confluent hypergeometric function,while equation (6.7.24)follows from relations between contiguous confluent hy- Press. pergeometric functions (equations 13.4.16 and 13.4.18 in Abramowitz and Stegun).Now the functions n satisfy the three-term recurrence relation (equation 13.4.15 in Abramowitz and Stegun) Programs 2n-1(E)=bnzn()+an+12n+1 (6.7.25) with SCIENTIFIC( bn =2(n+T) am+1=-[(n+1/2)2-v2] (6.7.26) 6 Following the steps leading to equation(5.5.18),we get the continued fraction CF2 =本于 1a2 COMPUTING (ISBN 1392 20 (6.7.27) from which (6.7.24)gives K+1/K and thus K/K. 10521 Temme's normalization condition is that Numerical 2-() +1/2 431 (6.7.28) 容 uction Recipes where C.=-)”u+12+ (6.7.29) n!T(w+1/2-n) Note that the Cn's can be determined by recursion: C0=1, Cn+1 =-antiCn (6.7.30) n+1 We use the condition(6.7.28)by finding S=Cn in (6.7.31) n=1 20 Then +1/21 1+ (6.7.32)
246 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). Modified Bessel Functions Steed’s method does not work for modified Bessel functions because in this case CF2 is purely imaginary and we have only three relations among the four functions. Temme [3] has given a normalization condition that provides the fourth relation. The Wronskian relation is W ≡ IνK ν − KνI ν = − 1 x (6.7.20) The continued fraction CF1 becomes fν ≡ I ν Iν = ν x + 1 2(ν + 1)/x + 1 2(ν + 2)/x + ··· (6.7.21) To get CF2 and the normalization condition in a convenient form, consider the sequence of confluent hypergeometric functions zn(x) = U(ν + 1/2 + n, 2ν + 1, 2x) (6.7.22) for fixed ν. Then Kν(x) = π1/2 (2x) ν e −xz0(x) (6.7.23) Kν+1(x) Kν (x) = 1 x ν + 1 2 + x + ν2 − 1 4 z1 z0 (6.7.24) Equation (6.7.23) is the standard expression for Kν in terms of a confluent hypergeometric function, while equation (6.7.24) follows from relations between contiguous confluent hypergeometric functions (equations 13.4.16 and 13.4.18 in Abramowitz and Stegun). Now the functions zn satisfy the three-term recurrence relation (equation 13.4.15 in Abramowitz and Stegun) zn−1(x) = bnzn(x) + an+1zn+1 (6.7.25) with bn = 2(n + x) an+1 = −[(n + 1/2)2 − ν2 ] (6.7.26) Following the steps leading to equation (5.5.18), we get the continued fraction CF2 z1 z0 = 1 b1 + a2 b2 + ··· (6.7.27) from which (6.7.24) gives Kν+1/Kν and thus K ν/Kν . Temme’s normalization condition is that ∞ n=0 Cnzn = 1 2x ν+1/2 (6.7.28) where Cn = (−1)n n! Γ(ν + 1/2 + n) Γ(ν + 1/2 − n) (6.7.29) Note that the Cn’s can be determined by recursion: C0 = 1, Cn+1 = − an+1 n + 1Cn (6.7.30) We use the condition (6.7.28) by finding S = ∞ n=1 Cn zn z0 (6.7.31) Then z0 = 1 2x ν+1/2 1 1 + S (6.7.32)
6.7 Bessel Functions of Fractional Order 247 and (6.7.23)gives K. Thompson and Barnett[4]have given a clever method of doing the sum (6.7.31) simultaneously with the forward evaluation of the continued fraction CF2.Suppose the continued fraction is being evaluated as =△hn (6.7.33) n=0 where the increments Ahn are being found by,e.g.,Steed's algorithm or the modified Lentz's algorithm of $5.2.Then the approximation to S keeping the first N terms can be found as Sw=Qn△hn (6.7.34) 8 n=1 Here Qm=∑Ckg (6.7.35) and q is found by recursion from qk+1=(qk-1-bkqk)/ak+1 (6.7.36) starting with o =0,g=1.For the case at hand,approximately three times as many terms RECIPES are needed to get S to converge as are needed simply for CF2 to converge. To find K and K+for small x we use series analogous to (6.7.14): K=∑kfik K+1=2∑ (6.7.37) k=0 正k-0 Press. Here 94=2/4 Programs k！ hk=-kfk+p张 IENTIFIC P&Pk-1 k-v (6.7.38) 9张-张-1 k+v MPUTING fk= kfk-1+Pk-1+9k-1 k2-v2 1920 (ISBN The initial values for the recurrences are P0= 1(）"r(1+0 Numerical 10621 2(2 90= T(1-w) (6.7.39) uction Recipes 43108 f6= UT sin cosh ar1(v)+ (outside Both the series for small and CF2 and the normalization relation(6.7.28)require North Software. 1/2.In both cases,therefore,we recurse I down to a value=u in this interval,find K there,and recurse K back up to the original value of v. The routine assumes>0.For negative v use the reflection formulas 显 Iv=Iv+-sin(UT)Kv (6.7.40) K-v=Kv Note that for large Ie,Ke-,and so these functions will overflow or underflow.It is often desirable to be able to compute the scaled quantitiese and eK Simply omitting the factor e*in equation(6.7.23)will ensure that all four quantities will have the appropriate scaling.If you also want to scale the four quantities for small when the series in equation (6.7.37)are used,you must multiply each series by e
6.7 Bessel Functions of Fractional Order 247 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). and (6.7.23) gives Kν . Thompson and Barnett [4] have given a clever method of doing the sum (6.7.31) simultaneously with the forward evaluation of the continued fraction CF2. Suppose the continued fraction is being evaluated as z1 z0 = ∞ n=0 ∆hn (6.7.33) where the increments ∆hn are being found by, e.g., Steed’s algorithm or the modified Lentz’s algorithm of §5.2. Then the approximation to S keeping the first N terms can be found as SN = N n=1 Qn∆hn (6.7.34) Here Qn = n k=1 Ckqk (6.7.35) and qk is found by recursion from qk+1 = (qk−1 − bkqk)/ak+1 (6.7.36) starting with q0 = 0, q1 = 1. For the case at hand, approximately three times as many terms are needed to get S to converge as are needed simply for CF2 to converge. To find Kν and Kν+1 for small x we use series analogous to (6.7.14): Kν = ∞ k=0 ckfk Kν+1 = 2 x ∞ k=0 ckhk (6.7.37) Here ck = (x2/4)k k! hk = −kfk + pk pk = pk−1 k − ν qk = qk−1 k + ν fk = kfk−1 + pk−1 + qk−1 k2 − ν2 (6.7.38) The initial values for the recurrences are p0 = 1 2 x 2 −ν Γ(1 + ν) q0 = 1 2 x 2 ν Γ(1 − ν) f0 = νπ sin νπ cosh σΓ1(ν) + sinh σ σ ln 2 x Γ2(ν) (6.7.39) Both the series for small x, and CF2 and the normalization relation (6.7.28) require |ν| ≤ 1/2. In both cases, therefore, we recurse Iν down to a value ν = µ in this interval, find Kµ there, and recurse Kν back up to the original value of ν. The routine assumes ν ≥ 0. For negative ν use the reflection formulas I−ν = Iν + 2 π sin(νπ) Kν K−ν = Kν (6.7.40) Note that for large x, Iν ∼ ex, Kν ∼ e−x, and so these functions will overflow or underflow. It is often desirable to be able to compute the scaled quantities e −xIν and exKν . Simply omitting the factor e−x in equation (6.7.23) will ensure that all four quantities will have the appropriate scaling. If you also want to scale the four quantities for small x when the series in equation (6.7.37) are used, you must multiply each series by e x.
248 Chapter 6.Special Functions #include #define EPS 1.0e-10 #define FPMIN 1.0e-30 #define MAXIT 10000 #define XMIN 2.0 #def1nePI3.141592653589793 void bessik(float x,float xnu,float *ri,float *rk,float *rip,float *rkp) Returns the modified Bessel functions ri=I,rk =K and their derivatives rip rkp =K,for positive x and for xnu =>0.The relative accuracy is within one or two significant digits of EPS.FPMIN is a number close to the machine's smallest floating-point number.All internal arithmetic is in double precision.To convert the entire routine to double precision,change the float declarations above to double and decrease EPS to 10-16.Also 5常 convert the function beschb. void beschb(double x,double *gam1,double *gam2,double *gampl, 18881892 double *gammi); void nrerror(char error_text☐)； 1800 int i,l,nl; double a,a1,b,c,d,del,del1,delh,dels,e,f,fact,fact2,ff,gam1,gam2, gammi,gampl,h,p,pimu,q,q1,q2,qnew,ril,ril1,rimu,ripi,ripl, from NUMERICAL RECIPES I ritemp,rk1,rkmu,rkmup,rktemp,s,sum,sum1,x2,xi,xi2,xmu,xmu2; if (x =1;1--){ (outside ritemp=fact*ril+ripl; fact -xi; North Software. ripl=fact*ritemp+ril; ril=ritemp; Ame f=ripl/ril; Now have unnormalized I and if (x XMIN){ Use series. x2=0.5*x; pimu=PI*xmu: fact =(fabs(pimu)<EPS 1.0 pimu/sin(pimu)); d=-1og(x2); e=xmu*d; fact2 =(fabs(e)<EPS 1.0 sinh(e)/e); beschb(xmu,&gam1,&gam2,&gampl,&gammi); Chebyshev evaluation of T1 and T2 ff=fact*(gami*cosh(e)+gam2*fact2*d); f0
248 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). #include #define EPS 1.0e-10 #define FPMIN 1.0e-30 #define MAXIT 10000 #define XMIN 2.0 #define PI 3.141592653589793 void bessik(float x, float xnu, float *ri, float *rk, float *rip, float *rkp) Returns the modified Bessel functions ri = Iν, rk = Kν and their derivatives rip = I ν, rkp = K ν, for positive x and for xnu = ν ≥ 0. The relative accuracy is within one or two significant digits of EPS. FPMIN is a number close to the machine’s smallest floating-point number. All internal arithmetic is in double precision. To convert the entire routine to double precision, change the float declarations above to double and decrease EPS to 10−16. Also convert the function beschb. { void beschb(double x, double *gam1, double *gam2, double *gampl, double *gammi); void nrerror(char error_text[]); int i,l,nl; double a,a1,b,c,d,del,del1,delh,dels,e,f,fact,fact2,ff,gam1,gam2, gammi,gampl,h,p,pimu,q,q1,q2,qnew,ril,ril1,rimu,rip1,ripl, ritemp,rk1,rkmu,rkmup,rktemp,s,sum,sum1,x2,xi,xi2,xmu,xmu2; if (x MAXIT) nrerror("x too large in bessik; try asymptotic expansion"); ril=FPMIN; Initialize Iν and I ν for downward reripl=h*ril; currence. ril1=ril; Store values for later rescaling. rip1=ripl; fact=xnu*xi; for (l=nl;l>=1;l--) { ritemp=fact*ril+ripl; fact -= xi; ripl=fact*ritemp+ril; ril=ritemp; } f=ripl/ril; Now have unnormalized Iµ and I µ. if (x < XMIN) { Use series. x2=0.5*x; pimu=PI*xmu; fact = (fabs(pimu) < EPS ? 1.0 : pimu/sin(pimu)); d = -log(x2); e=xmu*d; fact2 = (fabs(e) < EPS ? 1.0 : sinh(e)/e); beschb(xmu,&gam1,&gam2,&gampl,&gammi); Chebyshev evaluation of Γ1 and Γ2. ff=fact*(gam1*cosh(e)+gam2*fact2*d); f0
6.7 Bessel Functions of Fractional Order 249 sum=ff; e=exp(e); p=0.5*e/gampl; P0- q=0.5/(e*gamm1); qo. c=1.0; d=x2*x2; sumi=p; for (i=1;iMAXIT)nrerror("bessk series failed to converge"); rkmu=sum; 1872 rki=sum1*xi2; else Evaluate CF2 by Steed's algorithm to any server computer, t users to make one paper 1988-1992 by Cambridge University Press.Programs from NUMERICAL RECIPES IN b=2.0*(1.0+x): (S5.2),which is OK because there d=1.0/b; can be no zero denominators. (North h=delh=d; q1=0.0; Initializations for recurrence (6.7.35). THE q2=1.0; America a1=0.25-xmu2; 是 ART q=c=al; First term in equation (6.7.34) a=-a1; s=1.0+q*de1h; for (i=2;i<=MAXIT;i++){ send strictly prohibited. copy for their a-=2*(1-1): c -a*c/i; gnew=(q1-b*q2)/a; email Copyright(C) 91=g2; q2=qnew; q +c*qnew; b+=2.0; d=1.0/(b+a*d) delh=(b*d-1.0)*delh; OF SCIENTIFIC COMPUTING(ISBN 0-521 h +delh; dels=q*delh; v@cambri s +dels; if (fabs(dels/s)<EPS)break; Need only test convergence of sum since CF2 itself converges more quickly. To order Numerical Recipes books or 1988-1992 by Numerical Recipes 1-431085 if (i MAXIT)nrerror("bessik:failure to converge in cf2"); h=a1*h; (outside rkmu=sqrt(PI/(2.0*x))*exp(-x)/s; Omit the factor exp(-)to scale rk1=rkmu*(xmu+x+0.5-h)*xi; all the returned functions by exp(x) Software. forx≥XMIN rkmup=xmu*xi*rkmu-rk1; rimu=xi/(f*rkmu-rkmup); Get Iu from Wronskian *ri=(r1mu*r111)/r11; Scale original I and visit website *rip=(rimu*rip1)/ril; f machine for(1=1;1<=nl;i++) Upward recurrence of Kv rktemp=(xmu+i)*xi2*rk1+rkmu; rkmu=rk1; rk1-rktemp; *rk=rkmu; *rkp=xnu*xi*rkmu-rk1;
6.7 Bessel Functions of Fractional Order 249 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). sum=ff; e=exp(e); p=0.5*e/gampl; p0. q=0.5/(e*gammi); q0. c=1.0; d=x2*x2; sum1=p; for (i=1;i MAXIT) nrerror("bessk series failed to converge"); rkmu=sum; rk1=sum1*xi2; } else { Evaluate CF2 by Steed’s algorithm (§5.2), which is OK because there can be no zero denominators. b=2.0*(1.0+x); d=1.0/b; h=delh=d; q1=0.0; Initializations for recurrence (6.7.35). q2=1.0; a1=0.25-xmu2; q=c=a1; First term in equation (6.7.34). a = -a1; s=1.0+q*delh; for (i=2;i MAXIT) nrerror("bessik: failure to converge in cf2"); h=a1*h; rkmu=sqrt(PI/(2.0*x))*exp(-x)/s; Omit the factor exp(−x) to scale all the returned functions by exp(x) for x ≥ XMIN. rk1=rkmu*(xmu+x+0.5-h)*xi; } rkmup=xmu*xi*rkmu-rk1; rimu=xi/(f*rkmu-rkmup); Get Iµ from Wronskian. *ri=(rimu*ril1)/ril; Scale original Iν and I ν. *rip=(rimu*rip1)/ril; for (i=1;i<=nl;i++) { Upward recurrence of Kν. rktemp=(xmu+i)*xi2*rk1+rkmu; rkmu=rk1; rk1=rktemp; } *rk=rkmu; *rkp=xnu*xi*rkmu-rk1; }