Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 60
Текст из файла (страница 60)
We have to adopt a special method for small x, which we describe below. Forx not too small, we can ensure that x >∼ xtp by a stable recurrence of Jν and Jν downwardsto a value ν = µ <x,thusyieldingtheratiofatthislowervalueofν.Thisis the stableµ∼direction for the recurrence relation. The initial values for the recurrence areJν = 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 ofCF1.
Choosing the initial value of Jν very small minimizes the possibility of overflow duringthe recurrence. The recurrence relations areνJν−1 = Jν + Jνx(6.7.5)ν−1Jν−1 =Jν−1 − JνxOnce CF2 has been evaluated at ν = µ, then with the Wronskian (6.7.1) we have enoughrelations to solve for all four quantities. The formulas are simplified by introducing the quantityγ≡ThenJµ = ±p − fµqWq + γ(p − fµ )(6.7.6)1/2Jµ = fµ JµYµ = γJµqYµ = Yµ p +γ(6.7.7)(6.7.8)(6.7.9)(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 theoriginal value of ν.
For Jν and Jν , simply scale the values in (6.7.4) by the ratio of (6.7.7) tothe value found after applying the recurrence (6.7.5). The quantities Yν and Yν can be foundby starting with the values in (6.7.9) and (6.7.10) and using the stable upwards recurrenceYν+1 =2νYν − Yν−1x(6.7.11)νYν − Yν+1x(6.7.12)together with the relationYν =242Chapter 6.Special FunctionsNow turn to the case of small x, when CF2 is not suitable. Temme [2] has given agood method of evaluating Yν and Yν+1 , and hence Yν from (6.7.12), by series expansionsthat 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µ fromJµ =WYµ − 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∞∞2Yν = −ck gkYν+1 = −ck hk(6.7.14)xk=0k=0Here(−x2 /4)kck =(6.7.15)k!while the coefficients gk and hk are defined in terms of quantities pk , qk , and fk that canbe found by recursion:# νπ $2gk = fk + sin2qkν2hk = −kgk + pkpk−1k−νqk−1qk =k+νkfk−1 + pk−1 + qk−1fk =k2 − ν 2The initial values for the recurrences are1 # x $−νp0 =Γ(1 + ν)π 2#$1 x νq0 =Γ(1 − ν)π 2 22 νπsinh σcosh σΓ1 (ν) +Γ2 (ν)lnf0 =π sin νπσxwith 2σ = ν lnx111Γ1 (ν) =−2ν Γ(1 − ν)Γ(1 + ν)111Γ2 (ν) =+2 Γ(1 − ν)Γ(1 + ν)pk =(6.7.16)(6.7.17)(6.7.18)The whole point of writing the formulas in this way is that the potential problems as ν → 0can be controlled by evaluating νπ/ sin νπ, sinh σ/σ, and Γ1 carefully.
In particular, Temmegives Chebyshev expansions for Γ1 (ν) and Γ2 (ν). We have rearranged his expansion for Γ1to 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 formulasJ−ν = cos νπ Jν − sin νπ YνY−ν = sin νπ Jν + cos νπ Yν(6.7.19)The routine also assumes x > 0. For x < 0 the functions are in general complex, butexpressible in terms of functions with x > 0. For x = 0, Yν is singular.6.7 Bessel Functions of Fractional Order243Internal arithmetic in the routine is carried out in double precision. The complexarithmetic is carried out explicitly with real variables.#include <math.h>#include "nrutil.h"#define EPS 1.0e-10#define FPMIN 1.0e-30#define MAXIT 10000#define XMIN 2.0#define PI 3.141592653589793void 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ν , forpositive x and for xnu = ν ≥ 0.
The relative accuracy is within one or two significant digitsof 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 arithmeticis in double precision. To convert the entire routine to double precision, change the floatdeclarations 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 <= 0.0 || xnu < 0.0) nrerror("bad arguments in bessjy");nl=(x < XMIN ? (int)(xnu+0.5) : IMAX(0,(int)(xnu-x+1.5)));nl is the number of downward recurrences of the J’s and upward recurrences of Y ’s.
xmulies between −1/2 and 1/2 for x < XMIN, while it is chosen so that x is greater than theturning point for x ≥ XMIN.xmu=xnu-nl;xmu2=xmu*xmu;xi=1.0/x;xi2=2.0*xi;w=xi2/PI;The Wronskian.isign=1;Evaluate CF1 by modified Lentz’s method (§5.2).h=xnu*xi;isign keeps track of sign changes in the deif (h < FPMIN) h=FPMIN;nominator.b=xi2*xnu;d=0.0;c=h;for (i=1;i<=MAXIT;i++) {b += xi2;d=b-d;if (fabs(d) < FPMIN) d=FPMIN;c=b-1.0/c;if (fabs(c) < FPMIN) c=FPMIN;d=1.0/d;del=c*d;h=del*h;if (d < 0.0) isign = -isign;if (fabs(del-1.0) < EPS) break;}if (i > 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;244Chapter 6.Special Functionsrjpl=fact*rjtemp-rjl;rjl=rjtemp;}if (rjl == 0.0) rjl=EPS;f=rjpl/rjl;Now have unnormalized Jµ and Jµ .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=2.0/PI*fact*(gam1*cosh(e)+gam2*fact2*d);f0 .e=exp(e);p=e/(gampl*PI);p0 .q=1.0/(e*PI*gammi);q0 .pimu2=0.5*pimu;fact3 = (fabs(pimu2) < EPS ? 1.0 : sin(pimu2)/pimu2);r=PI*pimu2*fact3*fact3;c=1.0;d = -x2*x2;sum=ff+r*q;sum1=p;for (i=1;i<=MAXIT;i++) {ff=(i*ff+p+q)/(i*i-xmu2);c *= (d/i);p /= (i-xmu);q /= (i+xmu);del=c*(ff+r*q);sum += del;del1=c*p-i*del;sum1 += del1;if (fabs(del) < (1.0+fabs(sum))*EPS) break;}if (i > 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 Order245cr=br+cr*fact;ci=bi-ci*fact;if (fabs(cr)+fabs(ci) < FPMIN) cr=FPMIN;den=dr*dr+di*di;dr /= den;di /= -den;dlr=cr*dr-ci*di;dli=cr*di+ci*dr;temp=p*dlr-q*dli;q=p*dli+q*dlr;p=temp;if (fabs(dlr-1.0)+fabs(dli) < EPS) break;}if (i > 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;*rjp=rjp1*fact;for (i=1;i<=nl;i++) {rytemp=(xmu+i)*xi2*ry1-rymu;rymu=ry1;ry1=rytemp;}*ry=rymu;*ryp=xnu*xi*rymu-ry1;Scale original Jν and Jν .Upward recurrence of Yν .}#define NUSE1 5#define NUSE2 5void 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) and1/Γ(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;*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);}Multiply x by 2 to make range be −1 to 1,and then apply transformation for evaluating even Chebyshev series.246Chapter 6.Special FunctionsModified Bessel FunctionsSteed’s method does not work for modified Bessel functions because in this case CF2 ispurely imaginary and we have only three relations among the four functions.
Temme [3] hasgiven a normalization condition that provides the fourth relation.The Wronskian relation is1W ≡ Iν Kν − Kν Iν = −(6.7.20)xThe continued fraction CF1 becomesfν ≡1Iνν1= +···Iνx2(ν + 1)/x + 2(ν + 2)/x +(6.7.21)To get CF2 and the normalization condition in a convenient form, consider the sequenceof confluent hypergeometric functionsfor fixed ν.zn (x) = U (ν + 1/2 + n, 2ν + 1, 2x)(6.7.22)Kν (x) = π1/2 (2x)ν e−x z0 (x) Kν+1 (x)1 z1112ν + +x+ ν −=Kν (x)x24 z0(6.7.23)Then(6.7.24)Equation (6.7.23) is the standard expression for Kν in terms of a confluent hypergeometricfunction, 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). Nowthe functions zn satisfy the three-term recurrence relation (equation 13.4.15 in Abramowitzand Stegun)zn−1 (x) = bn zn (x) + an+1 zn+1with(6.7.25)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 CF2z1a21=···z0b1 + b2 +from which (6.7.24) gives Kν+1 /Kν and thus Kν /Kν .Temme’s normalization condition is that ν+1/2∞1Cn zn =2xn=0where(−1)n Γ(ν + 1/2 + n)Cn =n! Γ(ν + 1/2 − n)(6.7.27)(6.7.28)(6.7.29)Note that the Cn ’s can be determined by recursion:C0 = 1,Cn+1 = −an+1Cnn+1(6.7.30)We use the condition (6.7.28) by findingS=∞n=1Thenz0 =12xCnznz0ν+1/211+S(6.7.31)(6.7.32)6.7 Bessel Functions of Fractional Order247and (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 thecontinued fraction is being evaluated as∞z1=∆hn(6.7.33)z0n=0where the increments ∆hn are being found by, e.g., Steed’s algorithm or the modified Lentz’salgorithm of §5.2. Then the approximation to S keeping the first N terms can be found asSN =NQn ∆hn(6.7.34)n=1HereQn =n(6.7.35)Ck qkk=1and qk is found by recursion fromqk+1 = (qk−1 − bk qk )/ak+1(6.7.36)starting with q0 = 0, q1 = 1. For the case at hand, approximately three times as many termsare 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):∞∞2Kν =ck fkKν+1 =ck hk(6.7.37)xk=0k=0Here(x2 /4)kk!hk = −kfk + pkck =pk−1(6.7.38)k−νqk−1qk =k+νkfk−1 + pk−1 + qk−1fk =k2 − ν 2The initial values for the recurrences are1 # x $−νp0 =Γ(1 + ν)2 2#$1 x νq0 =Γ(1 − ν)(6.7.39)2 2 2νπsinh σcosh σΓ1 (ν) +lnΓ2 (ν)f0 =sin νπσxBoth the series for small x, and CF2 and the normalization relation (6.7.28) require|ν| ≤ 1/2.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
