Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 62
Текст из файла (страница 62)
A more robust computational procedure istherefore desirable, as follows:The associated Legendre functions satisfy numerous recurrence relations, tabulated in [1-2] . These are recurrences on l alone, on m alone, and on both l andm simultaneously. Most of the recurrences involving m are unstable, and sodangerous for numerical work. The following recurrence on l is, however, stable(compare 5.5.1):mm− (l + m − 1)Pl−2(l − m)Plm = x(2l − 1)Pl−1(6.8.7)It is useful because there is a closed-form expression for the starting value,mPm= (−1)m (2m − 1)!!(1 − x2 )m/2(6.8.8)(The notation n!! denotes the product of all odd integers less than or equal to n.)mUsing (6.8.7) with l = m + 1, and setting Pm−1= 0, we findmmPm+1= x(2m + 1)Pm(6.8.9)Equations (6.8.8) and (6.8.9) provide the two starting values required for (6.8.7)for general l.The function that implements this is254Chapter 6.Special Functions#include <math.h>float plgndr(int l, int m, float x)Computes the associated Legendre polynomial Plm (x).
Here m and l are integers satisfying0 ≤ m ≤ l, while x lies in the range −1 ≤ x ≤ 1.{void nrerror(char error_text[]);float fact,pll,pmm,pmmp1,somx2;int i,ll;if (m < 0 || m > l || fabs(x) > 1.0)nrerror("Bad arguments in routine plgndr");m.pmm=1.0;Compute Pmif (m > 0) {somx2=sqrt((1.0-x)*(1.0+x));fact=1.0;for (i=1;i<=m;i++) {pmm *= -fact*somx2;fact += 2.0;}}if (l == m)return pmm;m .else {Compute Pm+1pmmp1=x*(2*m+1)*pmm;if (l == (m+1))return pmmp1;else {Compute Plm , l > m + 1.for (ll=m+2;ll<=l;ll++) {pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m);pmm=pmmp1;pmmp1=pll;}return pll;}}}CITED REFERENCES AND FURTHER READING:Magnus, W., and Oberhettinger, F. 1949, Formulas and Theorems for the Functions of Mathematical Physics (New York: Chelsea), pp.
54ff. [1]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 8. [2]6.9 Fresnel Integrals, Cosine and Sine Integrals2556.9 Fresnel Integrals, Cosine and Sine IntegralsFresnel IntegralsThe two Fresnel integrals are defined by&C(x) =0x#π $t2 dt,cos2&S(x) =xsin0#π $t2 dt2(6.9.1)The most convenient way of evaluating these functions to arbitrary precision isto use power series for small x and a continued fraction for large x. The series are# π $4 x9# π $2 x5+−···C(x) = x −2 5 · 2!2 9 · 4!# π $3 x7# π $5 x11# π $ x3−+−···S(x) =2 3 · 1!2 7 · 3!2 11 · 5!(6.9.2)There is a complex continued fraction that yields both S(x) and C(x) simultaneously:√π(1 − i)x2(6.9.3)11 1/2 1 3/2 2···erfc z = √π z+ z+ z+ z+ z+2z1·23·41···= √π 2z 2 + 1 − 2z 2 + 5 − 2z 2 + 9 −(6.9.4)C(x) + iS(x) =1+ierf z,2z=wherez2eIn the last line we have converted the “standard” form of the continued fraction toits “even” form (see §5.2), which converges twice as fast.
We must be careful notto evaluate the alternating series (6.9.2) at too large a value of x; inspection of theterms shows that x = 1.5 is a good point to switch over to the continued fraction.Note that for large xC(x) ∼#π $11+sinx2 ,2 πx2S(x) ∼#π $11−cosx22 πx2(6.9.5)Thus the precision of the routine frenel may be limited by the precision of thelibrary routines for sine and cosine for large x.256Chapter 6.Special Functions#include <math.h>#include "complex.h"#define EPS 6.0e-8#define MAXIT 100#define FPMIN 1.0e-30#define XMIN 1.5#define PI 3.1415927#define PIBY2 (PI/2.0)Here EPS is the relative error; MAXIT is the maximum number of iterations allowed; FPMINis a number near the smallest representable floating-point number; XMIN is the dividing linebetween using the series and continued fraction.#define TRUE 1#define ONE Complex(1.0,0.0)void frenel(float x, float *s, float *c)Computes the Fresnel integrals S(x) and C(x) for all real x.{void nrerror(char error_text[]);int k,n,odd;float a,ax,fact,pix2,sign,sum,sumc,sums,term,test;fcomplex b,cc,d,h,del,cs;ax=fabs(x);if (ax < sqrt(FPMIN)) {*s=0.0;*c=ax;} else if (ax <= XMIN) {sum=sums=0.0;sumc=ax;sign=1.0;fact=PIBY2*ax*ax;odd=TRUE;term=ax;n=3;for (k=1;k<=MAXIT;k++) {term *= fact/k;sum += sign*term/n;test=fabs(sum)*EPS;if (odd) {sign = -sign;sums=sum;sum=sumc;} else {sumc=sum;sum=sums;}if (term < test) break;odd=!odd;n += 2;}if (k > MAXIT) nrerror("series failed*s=sums;*c=sumc;} else {pix2=PI*ax*ax;b=Complex(1.0,-pix2);cc=Complex(1.0/FPMIN,0.0);d=h=Cdiv(ONE,b);n = -1;for (k=2;k<=MAXIT;k++) {n += 2;a = -n*(n+1);b=Cadd(b,Complex(4.0,0.0));d=Cdiv(ONE,Cadd(RCmul(a,d),b));Special case: avoid failure of convergencetest because of underflow.Evaluate both series simultaneously.in frenel");Evaluate continued fraction by modifiedLentz’s method (§5.2).Denominators cannot be zero.6.9 Fresnel Integrals, Cosine and Sine Integrals257cc=Cadd(b,Cdiv(Complex(a,0.0),cc));del=Cmul(cc,d);h=Cmul(h,del);if (fabs(del.r-1.0)+fabs(del.i) < EPS) break;}if (k > MAXIT) nrerror("cf failed in frenel");h=Cmul(Complex(ax,-ax),h);cs=Cmul(Complex(0.5,0.5),Csub(ONE,Cmul(Complex(cos(0.5*pix2),sin(0.5*pix2)),h)));*c=cs.r;*s=cs.i;}if (x < 0.0) {*c = -(*c);*s = -(*s);}Use antisymmetry.}Cosine and Sine IntegralsThe cosine and sine integrals are defined by& xcos t − 1dtCi(x) = γ + ln x +t0& xsin tSi(x) =dtt0(6.9.6)Here γ ≈ 0.5772 .
. . is Euler’s constant. We only need a way to calculate thefunctions for x > 0, becauseSi(−x) = − Si(x),Ci(−x) = Ci(x) − iπ(6.9.7)Once again we can evaluate these functions by a judicious combination ofpower series and complex continued fraction. The series arex3x5+−···3 · 3! 5 · 5!x4x2+−···Ci(x) = γ + ln x + −2 · 2! 4 · 4!Si(x) = x −(6.9.8)The continued fraction for the exponential integral E1 (ix) isE1 (ix) = − Ci(x) + i[Si(x) − π/2]11221···= e−ixix + 1 + ix + 1 + ix +12221−ix···=e1 + ix − 3 + ix − 5 + ix −(6.9.9)The “even” form of the continued fraction is given in the last line and convergestwice as fast for about the same amount of computation.
A good crossover pointfrom the alternating series to the continued fraction is x = 2 in this case. As forthe Fresnel integrals, for large x the precision may be limited by the precision ofthe sine and cosine routines.258Chapter 6.#include <math.h>#include "complex.h"#define EPS 6.0e-8#define EULER 0.57721566#define MAXIT 100#define PIBY2 1.5707963#define FPMIN 1.0e-30#define TMIN 2.0#define TRUE 1#define ONE Complex(1.0,0.0)Special FunctionsRelative error, or absolute error near a zero of Ci(x).Euler’s constant γ.Maximum number of iterations allowed.π/2.Close to smallest representable floating-point number.Dividing line between using the series and continued fraction.void cisi(float x, float *ci, float *si)Computes the cosine and sine integrals Ci(x) and Si(x).
Ci(0) is returned as a large negativenumber and no error message is generated. For x < 0 the routine returns Ci(−x) and you mustsupply the −iπ yourself.{void nrerror(char error_text[]);int i,k,odd;float a,err,fact,sign,sum,sumc,sums,t,term;fcomplex h,b,c,d,del;t=fabs(x);if (t == 0.0) {Special case.*si=0.0;*ci = -1.0/FPMIN;return;}if (t > TMIN) {Evaluate continued fraction by modifiedb=Complex(1.0,t);Lentz’s method (§5.2).c=Complex(1.0/FPMIN,0.0);d=h=Cdiv(ONE,b);for (i=2;i<=MAXIT;i++) {a = -(i-1)*(i-1);b=Cadd(b,Complex(2.0,0.0));d=Cdiv(ONE,Cadd(RCmul(a,d),b));Denominators cannot be zero.c=Cadd(b,Cdiv(Complex(a,0.0),c));del=Cmul(c,d);h=Cmul(h,del);if (fabs(del.r-1.0)+fabs(del.i) < EPS) break;}if (i > MAXIT) nrerror("cf failed in cisi");h=Cmul(Complex(cos(t),-sin(t)),h);*ci = -h.r;*si=PIBY2+h.i;} else {Evaluate both series simultaneously.if (t < sqrt(FPMIN)) {Special case: avoid failure of convergencesumc=0.0;test because of underflow.sums=t;} else {sum=sums=sumc=0.0;sign=fact=1.0;odd=TRUE;for (k=1;k<=MAXIT;k++) {fact *= t/k;term=fact/k;sum += sign*term;err=term/fabs(sum);if (odd) {sign = -sign;sums=sum;sum=sumc;} else {sumc=sum;sum=sums;2596.10 Dawson’s Integral}if (err < EPS) break;odd=!odd;}if (k > MAXIT) nrerror("maxits exceeded in cisi");}*si=sums;*ci=sumc+log(t)+EULER;}if (x < 0.0) *si = -(*si);}CITED REFERENCES AND FURTHER READING:Stegun, I.A., and Zucker, R.
1976, Journal of Research of the National Bureau of Standards,vol. 80B, pp. 291–311; 1981, op. cit., vol. 86, pp. 661–686.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), Chapters 5 and 7.6.10 Dawson’s IntegralDawson’s Integral F (x) is defined by&−x2F (x) = ex2et dt(6.10.1)0The function can also be related to the complex error function by√i π −z2e[1 − erfc(−iz)] .F (z) =2(6.10.2)A remarkable approximation for F (x), due to Rybicki [1], is21 e−(z−nh)F (z) = lim √h→0πn(6.10.3)n oddWhat makes equation (6.10.3) unusual is that its accuracy increases exponentiallyas h gets small, so that quite moderate values of h (and correspondingly quite rapidconvergence of the series) give very accurate approximations.We will discuss the theory that leads to equation (6.10.3) later, in §13.11, asan interesting application of Fourier methods.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
