f6-3

6.3ExponentialIntegrals 215Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Chi-Square ProbabilityFunction P(χ2|ν)is definedas the probabilitythat the observedchi-squarefor a correct model should be less than a value χ2. (We will discuss the use of this function in Chapter15.) Itscomplement Q(χ2|ν)istheprobabilitythattheobservedchi-square will exceed the value χ2by chance evenfor a correct model. In both cases νis an integer,the numberof degreesof freedom. Thefunctionshavethe limitingvalues P(0|ν)=0 P(∞|ν)=1 ( 6.2.16 ) Q(0|ν)=1 Q(∞|ν)=0 ( 6.2.17 ) and the following relation to the incomplete gamma functions, P(χ2|ν)=P/parenleftbiggν 2,χ2 2/parenrightbigg =gammp/parenleftbiggν 2,χ2 2/parenrightbigg (6.2.18 ) Q(χ2|ν)=Q/parenleftbiggν 2,χ2 2/parenrightbigg =gammq/parenleftbiggν 2,χ2 2/parenrightbigg (6.2.19 ) CITED REFERENCES AND FURTHER READING: Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), Chapters 6, 7, and 26. Pearson, K. (ed.) 1951, Tables of the Incomplete Gamma Function (Cambridge: Cambridge University Press). 6.3 Exponential Integrals The standard definition of the exponential integral is En(x)=/integraldisplay∞ 1e−xt tndt, x > 0,n =0,1,... (6.3.1 ) The function defined by the principal value of the integral Ei(x)=−/integraldisplay∞ −xe−t tdt=/integraldisplayx −∞et tdt, x > 0( 6.3.2 ) is also called an exponential integral. Note that Ei(−x)is related to −E1(x)by analytic continuation. The function En(x)is a special case of the incomplete gamma function En(x)=xn−1Γ(1−n, x)( 6.3.3 ) 216 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).We can therefore use a similar strategy for evaluating it. The continued fraction — just equation (6.2.6) rewritten — converges for all x> 0: En(x)=e−x/parenleftbigg1 x+n 1+1 x+n+1 1+2 x+···/parenrightbigg (6.3.4 ) We use it in its more rapidly converging even form, En(x)=e−x/parenleftbigg1 x+n−1·n x+n+2−2(n+1 ) x+n+4−···/parenrightbigg (6.3.5 ) Thecontinuedfractiononlyreallyconvergesfastenoughtobeusefulfor x>∼1. For0<x<∼1, we can use the series representation En(x)=(−x)n−1 (n−1)![−lnx+ψ(n)]−∞/summationdisplay m=0 m/negationslash=n−1(−x)m (m−n+1 )m!(6.3.6 ) The quantity ψ(n)here is the digammafunction,givenforinteger argumentsby ψ(1) = −γ, ψ (n)=−γ+n−1/summationdisplay m=11 m(6.3.7 ) where γ=0.5772156649 ...is Euler’sconstant. We evaluatetheexpression(6.3.6) in order of ascending powers of x: En(x)=−/bracketleftbigg1 (1−n)−x (2−n)·1+x2 (3−n)(1·2)−··· +(−x)n−2 (−1)(n−2)!/bracketrightbigg +(−x)n−1 (n−1)![−lnx+ψ(n)]−/bracketleftbigg(−x)n 1·n!+(−x)n+1 2·(n+1 ) !+···/bracketrightbigg (6.3.8 ) The first square bracket is omitted when n=1. This method of evaluation has the advantage that for large nthe series converges before reaching the term containing ψ(n). Accordingly, one needs an algorithm for evaluating ψ(n)only for small n, n<∼20– 40. We use equation (6.3.7), although a table look-up would improve efficiency slightly. Amos[1]presents a careful discussion of the truncation error in evaluating equation (6.3.8), and gives a fairly elaborate termination criterion. We have found that simply stoppingwhen the last term addedis smaller than the requiredtolerance works about as well. Two special cases have to be handled separately: E0(x)=e−x x En(0) =1 n−1,n > 1(6.3.9 ) 6.3ExponentialIntegrals 217Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).The routine expintallows fast evaluation of En(x)to any accuracy EPS within the reach of your machine’s word length for floating-point numbers. Theonlymodificationrequiredforincreased accuracyis to supplyEuler’s constant with enough significant digits. Wrench [2]can provide you with the first 328 digits if necessary! FUNCTION expint(n,x) INTEGER n,MAXIT REAL expint,x,EPS,FPMIN,EULER PARAMETER (MAXIT=100,EPS=1.e-7,FPMIN=1.e-30,EULER=.5772156649) Evaluates the exponential integral En(x). Parameters: MAXIT is the maximum allowed number of iterations; EPS is the desired rel- ative error, not smaller than the machine precision; FPMIN is a number near the smallest representable floating-point number; EULER is Euler’s constant γ. INTEGER i,ii,nm1 REAL a,b,c,d,del,fact,h,psi nm1=n-1if(n.lt.0.or.x.lt.0..or.(x.eq.0..and.(n.eq.0.or.n.eq.1)))then pause ’bad arguments in expint’ else if(n.eq.0)then Special case. expint=exp(-x)/x else if(x.eq.0.)then Another special case. expint=1./nm1 else if(x.gt.1.)then Lentz’s algorithm ( §5.2). b=x+nc=1./FPMIN d=1./b h=ddo 11i=1,MAXIT a=-i*(nm1+i) b=b+2.d=1./(a*d+b) Denominators cannot be zero. c=b+a/c del=c*d h=h*delif(abs(del-1.).lt.EPS)then expint=h*exp(-x) return endif enddo 11 pause ’continued fraction failed in expint’ else Evaluate series. if(nm1.ne.0)then S e tfi r s tt e r m . expint=1./nm1 else expint=-log(x)-EULER endif fact=1. do13i=1,MAXIT fact=-fact*x/iif(i.ne.nm1)then del=-fact/(i-nm1) else psi=-EULER Compute ψ(n). do 12ii=1,nm1 psi=psi+1./ii enddo 12 del=fact*(-log(x)+psi) endif expint=expint+delif(abs(del).lt.abs(expint)*EPS) return enddo 13 218 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).pause ’series failed in expint’ endif return END A good algorithm for evaluating Eiis to use the power series for small xand the asymptotic series for large x. The power series is Ei(x)=γ+l nx+x 1·1!+x2 2·2!+··· (6.3.10 ) where γis Euler’s constant. The asymptotic expansion is Ei(x)∼ex x/parenleftbigg 1+1! x+2! x2+···/parenrightbigg (6.3.11 ) The lower limit for the use of the asymptotic expansion is approximately |lnEPS|, where EPSis the required relative error. FUNCTION ei(x) INTEGER MAXITREAL ei,x,EPS,EULER,FPMINPARAMETER (EPS=6.e-8,EULER=.57721566,MAXIT=100,FPMIN=1.e-30) Computes the exponential integral Ei(x)forx> 0. Parameters: EPS is the relative error, or absolute error near the zero of Eiatx=0.3725 ; EULER is Euler’s constant γ;MAXIT is the maximum number of iterations allowed; FPMIN is a number near the smallest representable floating-point number. INTEGER kREAL fact,prev,sum,termif(x.le.0.) pause ’bad argument in ei’ if(x.lt.FPMIN)then Special case: avoid failure of convergence test be- cause of underflow. ei=log(x)+EULER else if(x.le.-log(EPS))then Use power series. sum=0. fact=1.do 11k=1,MAXIT fact=fact*x/k term=fact/k sum=sum+termif(term.lt.EPS*sum)goto 1 enddo 11 pause ’series failed in ei’ 1 ei=sum+log(x)+EULER else Use asymptotic series. sum=0. Start with second term. term=1.do 12k=1,MAXIT prev=term term=term*k/x if(term.lt.EPS)goto 2 Since final sum is greater than one, term itself ap- proximates the relative error. if(term.lt.prev)then sum=sum+term Still converging: add new term. else sum=sum-prev Diverging: subtract previous term and exit. goto 2 endif enddo 12 2 ei=exp(x)*(1.+sum)/x endif 6.4IncompleteBetaFunction,Student’sDistribution,F-Distribution,CumulativeBinomialDistribution 219Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).return END CITED REFERENCES AND FURTHER READING: Stegun, I.A., and Zucker, R. 1974, Journal of Research of the National Bureau of Standards , vol. 78B, pp. 199–216; 1976, op. cit., vol. 80B, pp. 291–311. Amos D.E. 1980, ACM Transactions on Mathematical Software , vol. 6, pp. 365–377 [1]; also vol. 6, pp. 420–428. Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), Chapter 5. Wrench J.W. 1952, Mathematical Tablesand Other Aids to Computation , vol. 6, p. 255. [2] 6.4 Incomplete Beta Function, Student’s Distribution, F-Distribution, CumulativeBinomial Distribution The incomplete beta function is defined by Ix(a, b)≡Bx(a, b) B(a, b)≡1 B(a, b)/integraldisplayx 0ta−1(1−t)b−1dt (a, b > 0) (6.4.1 ) It has the limiting values I0(a, b)=0 I1(a, b)=1 ( 6.4.2 ) and the symmetry relation Ix(a, b)=1 −I1−x(b, a)( 6.4.3 ) Ifaandbare both rather greater than one, then Ix(a, b)rises from “near-zero” to “near-unity” quite sharply at about x=a/(a+b). Figure 6.4.1 plots the function for several pairs (a, b). The incomplete beta function has a series expansion Ix(a, b)=xa(1−x)b aB(a, b)/bracketleftBigg 1+∞/summationdisplay n=0B(a+1,n+1 ) B(a+b, n+1 )xn+1/bracketrightBigg , (6.4.4 ) butthisdoesnotprovetobeveryusefulinitsnumericalevaluation. (Note,however, that the beta functions in the coefficients can be evaluated for each value of nwith just the previousvalue and a few multiplies, using equations6.1.9 and 6.1.3.) The continued fraction representationproves to be much more useful, Ix(a, b)=xa(1−x)b aB(a, b)/bracketleftbigg1 1+d1 1+d2 1+···/bracketrightbigg (6.4.5 )