c6-7 (779508), страница 3
Текст из файла (страница 3)
In both cases, therefore, we recurse Iν down to a value ν = µ in this interval, findKµ there, and recurse Kν back up to the original value of ν.The routine assumes ν ≥ 0. For negative ν use the reflection formulas2I−ν = Iν + sin(νπ) Kνπ(6.7.40)K−ν = Kνpk =Note that for large x, Iν ∼ ex , Kν ∼ e−x , and so these functions will overflow orunderflow. It is often desirable to be able to compute the scaled quantities e−x Iν and ex Kν .Simply omitting the factor e−x in equation (6.7.23) will ensure that all four quantities willhave the appropriate scaling. If you also want to scale the four quantities for small x whenthe series in equation (6.7.37) are used, you must multiply each series by ex .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).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 thecontinued fraction is being evaluated as∞Xz1=∆hn(6.7.33)z0n=0248Chapter 6.Special Functions#include <math.h>#define EPS 1.0e-10#define FPMIN 1.0e-30#define MAXIT 10000#define XMIN 2.0#define PI 3.141592653589793if (x <= 0.0 || xnu < 0.0) nrerror("bad arguments in bessik");nl=(int)(xnu+0.5);nl is the number of downward rexmu=xnu-nl;currences of the I’s and upwardxmu2=xmu*xmu;recurrences of K’s. xmu lies bexi=1.0/x;tween −1/2 and 1/2.xi2=2.0*xi;h=xnu*xi;Evaluate CF1 by modified Lentz’sif (h < FPMIN) h=FPMIN;method (§5.2).b=xi2*xnu;d=0.0;c=h;for (i=1;i<=MAXIT;i++) {b += xi2;d=1.0/(b+d);Denominators cannot be zero here,c=b+1.0/c;so no need for special precaudel=c*d;tions.h=del*h;if (fabs(del-1.0) < EPS) break;}if (i > MAXIT) nrerror("x too large in bessik; try asymptotic expansion");ril=FPMIN;Initialize Iν and Iν0 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µ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=fact*(gam1*cosh(e)+gam2*fact2*d);f0 .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).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ν0 ,rkp = Kν0 , for positive x and for xnu = ν ≥ 0. The relative accuracy is within one or twosignificant digits of EPS.
FPMIN is a number close to the machine’s smallest floating-pointnumber. All internal arithmetic is in double precision. To convert the entire routine to doubleprecision, change the float declarations above to double and decrease EPS to 10−16 . Alsoconvert 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;6.7 Bessel Functions of Fractional Order249}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).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;i++) {ff=(i*ff+p+q)/(i*i-xmu2);c *= (d/i);p /= (i-xmu);q /= (i+xmu);del=c*ff;sum += del;del1=c*(p-i*ff);sum1 += del1;if (fabs(del) < fabs(sum)*EPS) break;}if (i > MAXIT) nrerror("bessk series failed to converge");rkmu=sum;rk1=sum1*xi2;} else {Evaluate CF2 by Steed’s algorithmb=2.0*(1.0+x);(§5.2), which is OK because thered=1.0/b;can be no zero denominators.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;i++) {a -= 2*(i-1);c = -a*c/i;qnew=(q1-b*q2)/a;q1=q2;q2=qnew;q += c*qnew;b += 2.0;d=1.0/(b+a*d);delh=(b*d-1.0)*delh;h += delh;dels=q*delh;s += dels;if (fabs(dels/s) < EPS) break;Need only test convergence of sum since CF2 itself converges more quickly.}if (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 scalerk1=rkmu*(xmu+x+0.5-h)*xi;all the returned functions by exp(x)}for x ≥ XMIN.rkmup=xmu*xi*rkmu-rk1;rimu=xi/(f*rkmu-rkmup);Get Iµ from Wronskian.*ri=(rimu*ril1)/ril;Scale original Iν and Iν0 .*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;250Chapter 6.Special FunctionsAiry Functionsso that Ai and Bi can be evaluated with a single call to bessik.The derivatives should not be evaluated by simply differentiating the above expressionsbecause of possible subtraction errors near x = 0.
Instead, use the equivalent expressionsxAi0 (x) = − √ K2/3 (z)π 3(6.7.45)21Bi0 (x) = x √ I2/3 (z) + K2/3 (z)π3The corresponding formulas for negative arguments are√ x1Ai(−x) =J1/3 (z) − √ Y1/3 (z)23√ x 1√ J1/3 (z) + Y1/3 (z)Bi(−x) = −23x1Ai0 (−x) =J2/3 (z) + √ Y2/3 (z)23x 10√ J2/3 (z) − Y2/3 (z)Bi (−x) =23#include <math.h>#define PI 3.1415927#define THIRD (1.0/3.0)#define TWOTHR (2.0*THIRD)#define ONOVRT 0.57735027void airy(float x, float *ai, float *bi, float *aip, float *bip)Returns Airy functions Ai(x), Bi(x), and their derivatives Ai0 (x), Bi0 (x).{void bessik(float x, float xnu, float *ri, float *rk, float *rip,float *rkp);void bessjy(float x, float xnu, float *rj, float *ry, float *rjp,float *ryp);float absx,ri,rip,rj,rjp,rk,rkp,rootx,ry,ryp,z;absx=fabs(x);rootx=sqrt(absx);z=TWOTHR*absx*rootx;if (x > 0.0) {bessik(z,THIRD,&ri,&rk,&rip,&rkp);*ai=rootx*ONOVRT*rk/PI;(6.7.46)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).For positive x, the Airy functions are defined byr1 xAi(x) =K1/3 (z)(6.7.41)π 3rx[I1/3 (z) + I−1/3 (z)](6.7.42)Bi(x) =3where2z = x3/2(6.7.43)3By using the reflection formula (6.7.40), we can convert (6.7.42) into the computationallymore useful form√21Bi(x) = x √ I1/3 (z) + K1/3 (z)(6.7.44)π36.7 Bessel Functions of Fractional Order251}Spherical Bessel FunctionsFor integer n, spherical Bessel functions are defined byrπjn (x) =Jn+(1/2) (x)2xrπyn (x) =Yn+(1/2) (x)2x(6.7.47)They can be evaluated by a call to bessjy, and the derivatives can safely be found fromthe derivatives of equation (6.7.47).Note that in the continued fraction CF2 in (6.7.3) just the first term survives for ν = 1/2.Thus one can make a very simple algorithm for spherical Bessel functions along the lines ofbessjy by always recursing jn down to n = 0, setting p and q from the first term in CF2, andthen recursing yn up.
No special series is required near x = 0. However, bessjy is alreadyso efficient that we have not bothered to provide an independent routine for spherical Bessels.#include <math.h>#define RTPIO2 1.2533141void sphbes(int n, float x, float *sj, float *sy, float *sjp, float *syp)0 (x), y0 (x) for integer n.Returns spherical Bessel functions jn (x), yn (x), and their derivatives jnn{void bessjy(float x, float xnu, float *rj, float *ry, float *rjp,float *ryp);void nrerror(char error_text[]);float factor,order,rj,rjp,ry,ryp;if (n < 0 || x <= 0.0) nrerror("bad arguments in sphbes");order=n+0.5;bessjy(x,order,&rj,&ry,&rjp,&ryp);factor=RTPIO2/sqrt(x);*sj=factor*rj;*sy=factor*ry;*sjp=factor*rjp-(*sj)/(2.0*x);*syp=factor*ryp-(*sy)/(2.0*x);}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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
