6.9 Fresnel Integrals,Cosine and Sine Integrals 255 6.9 Fresnel Integrals,Cosine and Sine Integrals Fresnel Integrals The two Fresnel integrals are defined by ca==(后到，s回=血(g9m (6.9.1) 8= The most convenient way of evaluating these functions to arbitrary precision is to use power series for small z and a continued fraction for large z.The series are C(x)=x- /T)2x5 /T)4x9 52+(594 (6.9.2) s-(-(（)+( (T3x7 )5x11 (Nort server University 2 11.50-… America Press. THE There is a complex continued fraction that yields both S()and C()simul- taneously: Progra ca+isa=l生af,=9a-0 (6.9.3) 天复 OF SCIENTIFIC( where 61 erfc z 1 11/213/22 2+2+z+2+2十 1920 COMPUTING (ISBN (6.9.4) 22 1 1.2 3·4 222+1-222+5-222+9- In the last line we have converted the "standard"form of the continued fraction to Numerical Recipes 10-621 4310-5 its"even"form (see 85.2),which converges twice as fast.We must be careful not to evaluate the alternating series(6.9.2)at too large a value of inspection of the (outside terms shows that x =1.5 is a good point to switch over to the continued fraction. Note that for large x North Software. c~+am() (6.9.5) Thus the precision of the routine frenel may be limited by the precision of the library routines for sine and cosine for large x
6.9 Fresnel Integrals, Cosine and Sine Integrals 255 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 6.9 Fresnel Integrals, Cosine and Sine Integrals Fresnel Integrals The two Fresnel integrals are defined by C(x) = x 0 cos π 2 t 2 dt, S(x) = x 0 sin π 2 t 2 dt (6.9.1) The most convenient way of evaluating these functions to arbitrary precision is to use power series for small x and a continued fraction for large x. The series are C(x) = x − π 2 2 x5 5 · 2! + π 2 4 x9 9 · 4! −··· S(x) = π 2 x3 3 · 1! − π 2 3 x7 7 · 3! + π 2 5 x11 11 · 5! −··· (6.9.2) There is a complex continued fraction that yields both S(x) and C(x) simultaneously: C(x) + iS(x) = 1 + i 2 erf z, z = √π 2 (1 − i)x (6.9.3) where ez2 erfc z = 1 √π 1 z + 1/2 z + 1 z + 3/2 z + 2 z + ··· = 2z √π 1 2z2 + 1 − 1 · 2 2z2 + 5 − 3 · 4 2z2 + 9 − ··· (6.9.4) In the last line we have converted the “standard” form of the continued fraction to its “even” form (see §5.2), which converges twice as fast. We must be careful not to evaluate the alternating series (6.9.2) at too large a value of x; inspection of the terms shows that x = 1.5 is a good point to switch over to the continued fraction. Note that for large x C(x) ∼ 1 2 + 1 πx sin π 2 x2 , S(x) ∼ 1 2 − 1 πx cos π 2 x2 (6.9.5) Thus the precision of the routine frenel may be limited by the precision of the library routines for sine and cosine for large x.
256 Chapter 6. Special Functions #include #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;FPMIN is a number near the smallest representable floating-point number:XMIN is the dividing line between using the series and continued fraction #define TRUE 1 #define ONE Complex(1.0,0.0) void frenel(float x,float *s,float *c) 83g Computes the Fresnel integrals S(x)and C(z)for all real x. granted for 19881992 void nrerror(char error_text[]); 1600 (including this one) int k,n,odd; float a,ax,fact,pix2,sign,sum,sumc,sums,term,test; fcomplex b,cc,d,h,del,cs; 7423 from NUMERICAL RECIPES IN ax=fabs(x); if (ax sqrt(FPMIN)){ Special case:avoid failure of convergence *8=0.0; test because of underflow. 米C=aX: to any server computer, else if (ax <XMIN){ Evaluate both series simultaneously. (North America sum=sums=0.0; one paper /Cambridge University Press.Programs THE ART sumc=ax; s1gn=1.0: only). fact=PIBY2*ax*ax; odd=TRUE; send copy for their term=ax; n=3; for (k=1;k<=MAXIT;k++){ email term *fact/k; sum sign*term/n; test=fabs(sum)*EPS; if (odd){ sign =-sign; sums=sum; OF SCIENTIFIC COMPUTING (ISBN 0-521 sum=sumc; else sumc=sum: sum=sums; @cambridge.org 1988-1992 by Numerical Recipes if (term test)break; -431085 odd=!odd; n+=2; if (k MAXIT)nrerror("series failed in frenel"); (outside North Software. *s=sums; *c=sumc; else Evaluate continued fraction by modified pix2-PI*ax*ax; Lentz's method ($5.2). 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)); Denominators cannot be zero
256 Chapter 6. Special Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). #include #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; FPMIN is a number near the smallest representable floating-point number; XMIN is the dividing line between 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 MAXIT) nrerror("series failed in frenel"); *s=sums; *c=sumc; } else { Evaluate continued fraction by modified pix2=PI*ax*ax; Lentz’s method (§5.2). 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)); Denominators cannot be zero
6.9 Fresnel Integrals,Cosine and Sine Integrals 257 cc=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)0.because OF SCIENTIFIC Si(-z)=-Si(x), Ci(-x)=Ci(x)-iπ (6.9.7) 6 Once again we can evaluate these functions by a judicious combination of power series and complex continued fraction.The series are COMPUTING (ISBN 188810920 Si(x)=- 3.3! 5.51 23 (6.9.8) Ci(z)=y+Inz+ 22+ 4.4 Numerical Recipes 10-621 -43106 The continued fraction for the exponential integral E1(ix)is (outside B1(ix)=-Ci(x)+[Si(x)-π/2] Software. =e-i证 1112 2 ix+1+ix+1+ix+ (6.9.9) 1 12 22 1+ix-3+ix-5+ix- The "even"form of the continued fraction is given in the last line and converges twice as fast for about the same amount of computation.A good crossover point from the alternating series to the continued fraction is =2 in this case.As for the Fresnel integrals,for large the precision may be limited by the precision of the sine and cosine routines
6.9 Fresnel Integrals, Cosine and Sine Integrals 257 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). cc=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) 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, because Si(−x) = − Si(x), Ci(−x) = Ci(x) − iπ (6.9.7) Once again we can evaluate these functions by a judicious combination of power series and complex continued fraction. The series are Si(x) = x − x3 3 · 3! + x5 5 · 5! −··· Ci(x) = γ + ln x + − x2 2 · 2! + x4 4 · 4! −··· (6.9.8) The continued fraction for the exponential integral E1(ix) is E1(ix) = − Ci(x) + i[Si(x) − π/2] = e−ix 1 ix + 1 1 + 1 ix + 2 1 + 2 ix + ··· = e−ix 1 1 + ix − 12 3 + ix − 22 5 + ix − ··· (6.9.9) The “even” form of the continued fraction is given in the last line and converges twice as fast for about the same amount of computation. A good crossover point from the alternating series to the continued fraction is x = 2 in this case. As for the Fresnel integrals, for large x the precision may be limited by the precision of the sine and cosine routines
258 Chapter 6. Special Functions #include #include "complex.h" #define EPS 6.0e-8 Relative error,or absolute error near a zero of Ci(). #define EULER 0.57721566 Euler's constant y. #define MAXIT 100 Maximum number of iterations allowed. #define PIBY2 1.5707963 π/2. #define FPMIN 1.0e-30 Close to smallest representable floating-point number. #define TMIN 2.0 Dividing line between using the series and continued frac- #define TRUE 1 tion. #define ONE Complex(1.0,0.0) void cisi(float x,float *ci,float *si) Computes the cosine and sine integrals Ci(z)and Si(x).Ci(0)is returned as a large negative number and no error message is generated.For xTMIN){ Evaluate continued fraction by modified one paper ART b=Complex(1.0,t); Lentz's method ($5.2). c=Complex(1.0/FPMIN,0.0); 是 d=h=Cdiv(ONE,b); Programs for (i=2;iMAXIT)nrerror("cf failed in cisi"); h=Cmul(Complex(cos(t),-sin(t)),h); tci =-h.r; 10-621 米si=PIBY2+h.1: else Evaluate both series simultaneously. 、海 if (t sqrt(FPMIN)){ Special case:avoid failure of convergence 1988-1992 by Numerical Recipes -43108 sumc=0.0; test because of underflow. sumset e1se (outside sum=sums=sumc=0.0; sign=fact=1.0; North Software. odd=TRUE; for (k=1;k<=MAXIT;k++){ Amer fact *=t/k; term=fact/k; sum +sign*term; err=term/fabs(sum) if (odd){ s1g即=-s1g即： sums=sum; sum=sumc else sumc=sum; sum=sums
258 Chapter 6. Special Functions Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). #include #include "complex.h" #define EPS 6.0e-8 Relative error, or absolute error near a zero of Ci(x). #define EULER 0.57721566 Euler’s constant γ. #define MAXIT 100 Maximum number of iterations allowed. #define PIBY2 1.5707963 π/2. #define FPMIN 1.0e-30 Close to smallest representable floating-point number. #define TMIN 2.0 Dividing line between using the series and continued frac- #define TRUE 1 tion. #define ONE Complex(1.0,0.0) 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 negative number and no error message is generated. For x TMIN) { Evaluate continued fraction by modified b=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) 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 convergence sumc=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;
6.10 Dawson's Integral 259 if (err EPS)break; odd=！odd; if (k MAXIT)nrerror("maxits exceeded in cisi"); *si=sumsi *ci=sumc+log(t)+EULER; 1f(x<0.0)*s1=-(*s1); 83 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,0p.ct,ol.86,pp.661-686. Abramowitz,M.,and Stegun,I.A.1964,Handbook of Mathematical Functions,Applied Mathe- matics Series,Volume 55 (Washington:National Bureau of Standards;reprinted 1968 by Dover Publications,New York),Chapters 5 and 7. (Nort server 6.10 Dawson's Integral America computer, University Press. THE Dawson's Integral F(r)is defined by Progra F(a)-efed (6.10.1) OF SCIENTIFIC The function can also be related to the complex error function by iv示 F(z)= e[1-erfe(-iz)]. (6.10.2) A remarkable approximation for F(z),due to Rybicki [1],is 、商◇ 6 Numerica 10621 1 e-(s-nh) F()=√示naaa (6.10.3) n 431 What makes equation (6.10.3)unusual is that its accuracy increases exponentially (outside Recipes as h gets small,so that quite moderate values of h(and correspondingly quite rapid convergence of the series)give very accurate approximations. North We will discuss the theory that leads to equation(6.10.3)later,in $13.11,as an interesting application of Fourier methods.Here we simply implement a routine for real values of based on the formula. It is first convenient to shift the summation index to center it approximately on the maximum of the exponential term.Define no to be the even integer nearest to x/h,and zo≡noh,x'≡x-xo,andn'≡n-no,so that N e-(x'-n'h2 F(x)≈ (6.10.4) n'+no
6.10 Dawson’s Integral 259 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). } if (err 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 by Dover Publications, New York), Chapters 5 and 7. 6.10 Dawson’s Integral Dawson’s Integral F(x) is defined by F(x) = e−x2 x 0 et 2 dt (6.10.1) The function can also be related to the complex error function by F(z) = i √π 2 e−z2 [1 − erfc(−iz)] . (6.10.2) A remarkable approximation for F(z), due to Rybicki [1], is F(z) = limh→0 1 √π n odd e−(z−nh)2 n (6.10.3) What makes equation (6.10.3) unusual is that its accuracy increases exponentially as h gets small, so that quite moderate values of h (and correspondingly quite rapid convergence of the series) give very accurate approximations. We will discuss the theory that leads to equation (6.10.3) later, in §13.11, as an interesting application of Fourier methods. Here we simply implement a routine for real values of x based on the formula. It is first convenient to shift the summation index to center it approximately on the maximum of the exponential term. Define n0 to be the even integer nearest to x/h, and x0 ≡ n0h, x ≡ x − x0, and n ≡ n − n0, so that F(x) ≈ 1 √π N n=−N n odd e−(x −n h)2 n + n0 , (6.10.4)