Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 61
Текст из файла (страница 61)
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 .248Chapter 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.141592653589793void 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 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;if (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ν 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 Order249sum=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ν .*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 FunctionsFor positive x, the Airy functions are defined by)1 xAi(x) =K1/3 (z)(6.7.41)π 3)x[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)π3so 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 expressionsxAi (x) = − √ K2/3 (z)π 3(6.7.45)21Bi (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) = −23x1Ai (−x) =J2/3 (z) + √ Y2/3 (z)231x√ 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 Ai (x), Bi (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)6.7 Bessel Functions of Fractional Order251*bi=rootx*(rk/PI+2.0*ONOVRT*ri);bessik(z,TWOTHR,&ri,&rk,&rip,&rkp);*aip = -x*ONOVRT*rk/PI;*bip=x*(rk/PI+2.0*ONOVRT*ri);} else if (x < 0.0) {bessjy(z,THIRD,&rj,&ry,&rjp,&ryp);*ai=0.5*rootx*(rj-ONOVRT*ry);*bi = -0.5*rootx*(ry+ONOVRT*rj);bessjy(z,TWOTHR,&rj,&ry,&rjp,&ryp);*aip=0.5*absx*(ONOVRT*ry+rj);*bip=0.5*absx*(ONOVRT*rj-ry);} else {Case x = 0.*ai=0.35502805;*bi=(*ai)/ONOVRT;*aip = -0.25881940;*bip = -(*aip)/ONOVRT;}}Spherical Bessel FunctionsFor integer n, spherical Bessel functions are defined by)πjn (x) =Jn+(1/2) (x)2x)π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) (x), y (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);}252Chapter 6.Special FunctionsCITED REFERENCES AND FURTHER READING:Barnett, A.R., Feng, D.H., Steed, J.W., and Goldfarb, L.J.B.
1974, Computer Physics Communications, vol. 8, pp. 377–395. [1]Temme, N.M. 1976, Journal of Computational Physics, vol. 21, pp. 343–350 [2]; 1975, op. cit.,vol. 19, pp. 324–337. [3]Thompson, I.J., and Barnett, A.R. 1987, Computer Physics Communications, vol. 47, pp. 245–257. [4]Barnett, A.R. 1981, Computer Physics Communications, vol.
21, pp. 297–314.Thompson, I.J., and Barnett, A.R. 1986, Journal of Computational Physics, vol. 64, pp. 490–509.Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), Chapter 10.6.8 Spherical HarmonicsSpherical harmonics occur in a large variety of physical problems, for example, whenever a wave equation, or Laplace’s equation, is solved by separation of variables in spherical coordinates.
The spherical harmonic Ylm (θ, φ),−l ≤ m ≤ l, is a function of the two coordinates θ, φ on the surface of a sphere.The spherical harmonics are orthogonal for different l and m, and they arenormalized so that their integrated square over the sphere is unity:&&2π1dφ−10d(cos θ)Yl m *(θ, φ)Ylm (θ, φ) = δl l δm m(6.8.1)Here asterisk denotes complex conjugation.Mathematically, the spherical harmonics are related to associated Legendrepolynomials by the equation*Ylm (θ, φ) =2l + 1 (l − m)! mP (cos θ)eimφ4π (l + m)! l(6.8.2)By using the relationYl,−m (θ, φ) = (−1)m Ylm *(θ, φ)(6.8.3)we can always relate a spherical harmonic to an associated Legendre polynomialwith m ≥ 0. With x ≡ cos θ, these are defined in terms of the ordinary Legendrepolynomials (cf. §4.5 and §5.5) byPlm (x) = (−1)m (1 − x2 )m/2dmPl (x)dxm(6.8.4)2536.8 Spherical HarmonicsThe first few associated Legendre polynomials, and their corresponding normalized spherical harmonics, are/1P00 (x) =1Y00 =4π/3P11 (x) = − (1 − x2 )1/2Y11 = − 8πsin θeiφ/3P10 (x) =xY10 =cos θ4π/15P22 (x) = 3 (1 − x2 )Y22 = 14 2πsin2 θe2iφ/15P21 (x) = −3 (1 − x2 )1/2 xY21 = − 8πsin θ cos θeiφ/5 312P20 (x) = 12 (3x2 − 1)Y20 =4π ( 2 cos θ − 2 )(6.8.5)There are many bad ways to evaluate associated Legendre polynomials numerically.
For example, there are explicit expressions, such as(l − m)(m + l + 1) 1 − x(−1)m (l + m)!(1 − x2 )m/2 1 −Plm (x) = m2 m!(l − m)!1!(m + 1)2!2(l − m)(l − m − 1)(m + l + 1)(m + l + 2) 1 − x−···+2!(m + 1)(m + 2)2(6.8.6)l−mwhere the polynomial continues up through the term in (1 − x). (See [1] forthis and related formulas.) This is not a satisfactory method because evaluationof the polynomial involves delicate cancellations between successive terms, whichalternate in sign. For large l, the individual terms in the polynomial become verymuch larger than their sum, and all accuracy is lost.In practice, (6.8.6) can be used only in single precision (32-bit) for l upto 6 or 8, and in double precision (64-bit) for l up to 15 or 18, depending onthe precision required for the answer.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
