c6-7 (779508), страница 2
Текст из файла (страница 2)
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).WYµ0 − Yµ fµ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.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ν0 , ryp = Yν0 , 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ν0 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;Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).#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.141592653589793244Chapter 6.Special Functionsrjpl=fact*rjtemp-rjl;rjl=rjtemp;Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).}if (rjl == 0.0) rjl=EPS;f=rjpl/rjl;Now have unnormalized Jµ and Jµ0 .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 Order245}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ν0 .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.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).cr=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;246Chapter 6.Special FunctionsModified Bessel Functionsfν ≡1Iν0ν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 z111ν + + x + ν2 −=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+1(6.7.25)withbn = 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ν0 /Kν .Temme’s normalization condition is that ν+1/2∞X1Cn 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:Cn+1 = −C0 = 1,an+1Cnn+1(6.7.30)We use the condition (6.7.28) by findingS=∞Xn=1Thenz0 =12xCnznz0ν+1/211+S(6.7.31)(6.7.32)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Steed’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ν0 − Kν Iν0 = −(6.7.20)xThe continued fraction CF1 becomes6.7 Bessel Functions of Fractional Order247where 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 =NXQn ∆hn(6.7.34)n=1HereQn =nXCk qk(6.7.35)k=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):∞∞X2XKν =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 xq0 =Γ(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}
