Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 63
Текст из файла (страница 63)
Here we simply implement a routinebased on the formula.It is first convenient to shift the summation index to center it approximately onthe maximum of the exponential term. Define n0 to be the even integer nearest tox/h, and x0 ≡ n0 h, x ≡ x − x0 , and n ≡ n − n0 , so that2N1 e−(x −n h),F (x) ≈ √π n =−N n + n0n odd(6.10.4)260Chapter 6.Special Functionswhere the approximate equality is accurate when h is sufficiently small and N issufficiently large. The computation of this formula can be greatly speeded up ifwe note that22e−(x −n h) = e−x e−(n h)2#e2x h$n.(6.10.5)The first factor is computed once, the second is an array of constants to be stored,and the third can be computed recursively, so that only two exponentials need be2evaluated. Advantage is also taken of the symmetry of the coefficients e−(n h) bybreaking the summation up into positive and negative values of n separately.In the following routine, the choices h = 0.4 and N = 11 are made.
Becauseof the symmetry of the summations and the restriction to odd values of n, the limitson the do loops are 1 to 6. The accuracy of the result in this float version is about2 × 10−7 . In order to maintain relative accuracy near x = 0, where F (x) vanishes,the program branches to the evaluation of the power series [2] for F (x), for |x| < 0.2.#include <math.h>#include "nrutil.h"#define NMAX 6#define H 0.4#define A1 (2.0/3.0)#define A2 0.4#define A3 (2.0/7.0)float dawson(float x)'Returns Dawson’s integral F (x) = exp(−x2 ) 0x exp(t2 )dt for any real x.{int i,n0;float d1,d2,e1,e2,sum,x2,xp,xx,ans;static float c[NMAX+1];static int init = 0;Flag is 0 if we need to initialize, else 1.if (init == 0) {init=1;for (i=1;i<=NMAX;i++) c[i]=exp(-SQR((2.0*i-1.0)*H));}if (fabs(x) < 0.2) {Use series expansion.x2=x*x;ans=x*(1.0-A1*x2*(1.0-A2*x2*(1.0-A3*x2)));} else {Use sampling theorem representation.xx=fabs(x);n0=2*(int)(0.5*xx/H+0.5);xp=xx-n0*H;e1=exp(2.0*xp*H);e2=e1*e1;d1=n0+1;d2=d1-2.0;sum=0.0;for (i=1;i<=NMAX;i++,d1+=2.0,d2-=2.0,e1*=e2)sum += c[i]*(e1/d1+1.0/(d2*e1));√ans=0.5641895835*SIGN(exp(-xp*xp),x)*sum;Constant is 1/ π.}return ans;}6.11 Elliptic Integrals and Jacobian Elliptic Functions261Other methods for computing Dawson’s integral are also known [2,3] .CITED REFERENCES AND FURTHER READING:Rybicki, G.B.
1989, Computers in Physics, vol. 3, no. 2, pp. 85–87. [1]Cody, W.J., Pociorek, K.A., and Thatcher, H.C. 1970, Mathematics of Computation, vol. 24,pp. 171–178. [2]McCabe, J.H. 1974, Mathematics of Computation, vol. 28, pp. 811–816. [3]6.11 Elliptic Integrals and Jacobian EllipticFunctionsElliptic integrals occur in many applications, because any integral of the form&R(t, s) dt(6.11.1)where R is a rational function of t and s, and s is the square root of a cubic orquartic polynomial in t, can be evaluated in terms of elliptic integrals.
Standardreferences [1] describe how to carry out the reduction, which was originally doneby Legendre. Legendre showed that only three basic elliptic integrals are required.The simplest of these is& xdt((6.11.2)I1 =(a1 + b1 t)(a2 + b2 t)(a3 + b3 t)(a4 + b4 t)ywhere we have written the quartic s2 in factored form. In standard integral tables [2],one of the limits of integration is always a zero of the quartic, while the other limitlies closer than the next zero, so that there is no singularity within the interval.
Toevaluate I1 , we simply break the interval [y, x] into subintervals, each of whicheither begins or ends on a singularity. The tables, therefore, need only distinguishthe eight cases in which each of the four zeros (ordered according to size) appears asthe upper or lower limit of integration. In addition, when one of the b’s in (6.11.2)tends to zero, the quartic reduces to a cubic, with the largest or smallest singularitymoving to ±∞; this leads to eight more cases (actually just special cases of the firsteight). The sixteen cases in total are then usually tabulated in terms of Legendre’sstandard elliptic integral of the 1st kind, which we will define below.
By a change ofthe variable of integration t, the zeros of the quartic are mapped to standard locationson the real axis. Then only two dimensionless parameters are needed to tabulateLegendre’s integral. However, the symmetry of the original integral (6.11.2) underpermutation of the roots is concealed in Legendre’s notation. We will get back toLegendre’s notation below. But first, here is a better way:Carlson [3] has given a new definition of a standard elliptic integral of the first kind,&dt1 ∞(RF (x, y, z) =(6.11.3)2 0(t + x)(t + y)(t + z)262Chapter 6.Special Functionswhere x, y, and z are nonnegative and at most one is zero. By standardizing the rangeof integration, he retains permutation symmetry for the zeros.
(Weierstrass’ canonical formalso has this property.) Carlson first shows that when x or y is a zero of the quartic in(6.11.2), the integral I1 can be written in terms of RF in a form that is symmetric underpermutation of the remaining three zeros. In the general case when neither x nor y is azero, two such RF functions can be combined into a single one by an addition theorem,leading to the fundamental formula222, U13, U14)I1 = 2RF (U12(6.11.4)whereUij = (Xi Xj Yk Ym + Yi Yj Xk Xm )/(x − y)(6.11.5)Xi = (ai + bi x)1/2 ,(6.11.6)Yi = (ai + bi y)1/2and i, j, k, m is any permutation of 1, 2, 3, 4. A short-cut in evaluating these expressions is22U13= U12− (a1 b4 − a4 b1 )(a2 b3 − a3 b2 )22U14= U12− (a1 b3 − a3 b1 )(a2 b4 − a4 b2 )(6.11.7)The U ’s correspond to the three ways of pairing the four zeros, and I1 is thus manifestlysymmetric under permutation of the zeros.
Equation (6.11.4) therefore reproduces all sixteencases when one limit is a zero, and also includes the cases when neither limit is a zero.Thus Carlson’s function allows arbitrary ranges of integration and arbitrary positions ofthe branch points of the integrand relative to the interval of integration. To handle ellipticintegrals of the second and third kind, Carlson defines the standard integral of the third kind as&dt3 ∞(RJ (x, y, z, p) =(6.11.8)2 0 (t + p) (t + x)(t + y)(t + z)which is symmetric in x, y, and z.
The degenerate case when two arguments are equalis denotedRD (x, y, z) = RJ (x, y, z, z)(6.11.9)and is symmetric in x and y. The function RD replaces Legendre’s integral of the secondkind. The degenerate form of RF is denotedRC (x, y) = RF (x, y, y)(6.11.10)It embraces logarithmic, inverse circular, and inverse hyperbolic functions.Carlson [4-7] gives integral tables in terms of the exponents of the linear factors ofthe quartic in (6.11.1). For example, the integral where the exponents are ( 12 , 12 ,− 12 ,− 32 )can be expressed as a single integral in terms of RD ; it accounts for 144 separate cases inGradshteyn and Ryzhik [2]!Refer to Carlson’s papers [3-7] for some of the practical details in reducing ellipticintegrals to his standard forms, such as handling complex conjugate zeros.Turn now to the numerical evaluation of elliptic integrals.
The traditional methods [8]are Gauss or Landen transformations. Descending transformations decrease the modulusk of the Legendre integrals towards zero, increasing transformations increase it towardsunity. In these limits the functions have simple analytic expressions. While these methodsconverge quadratically and are quite satisfactory for integrals of the first and second kinds,they generally lead to loss of significant figures in certain regimes for integrals of the thirdkind. Carlson’s algorithms [9,10] , by contrast, provide a unified method for all three kindswith no significant cancellations.The key ingredient in these algorithms is the duplication theorem:RF (x, y, z) = 2RF (x + λ, y + λ, z + λ)x+λ y+λ z+λ,,= RF444(6.11.11)2636.11 Elliptic Integrals and Jacobian Elliptic Functionswhereλ = (xy)1/2 + (xz)1/2 + (yz)1/2(6.11.12)This theorem can be proved by a simple change of variable of integration [11].
Equation(6.11.11) is iterated until the arguments of RF are nearly equal. For equal arguments we haveRF (x, x, x) = x−1/2(6.11.13)When the arguments are close enough, the function is evaluated from a fixed Taylor expansionabout (6.11.13) through fifth-order terms. While the iterative part of the algorithm is onlylinearly convergent, the error ultimately decreases by a factor of 46 = 4096 for each iteration.Typically only two or three iterations are required, perhaps six or seven if the initial valuesof the arguments have huge ratios.
We list the algorithm for RF here, and refer you toCarlson’s paper [9] for the other cases.Stage 1: For n = 0, 1, 2, . . . computeµn = (xn + yn + zn )/3Xn = 1 − (xn /µn ),Yn = 1 − (yn /µn ),Zn = 1 − (zn /µn )n = max(|Xn |, |Yn |, |Zn |)If n < tol go to Stage 2; else computeλn = (xn yn )1/2 + (xn zn )1/2 + (yn zn )1/2xn+1 = (xn + λn )/4,yn+1 = (yn + λn )/4,zn+1 = (zn + λn )/4and repeat this stage.Stage 2: ComputeE2 = Xn Yn − Zn2 ,RF = (1 −110 E2E3 = Xn Yn Zn+114 E3+2124 E2−1/2344 E2 E3 )/(µn )In some applications the argument p in RJ or the argument y in RC is negative, and theCauchy principal value of the integral is required. This is easily handled by using the formulasRJ (x, y,z, p) =[(γ − y)RJ (x, y, z, γ) − 3RF (x, y, z) + 3RC (xz/y, pγ/y)] /(y − p)(6.11.14)whereγ≡y+(z − y)(y − x)y−p(6.11.15)is positive if p is negative, andRC (x, y) =xx−y1/2RC (x − y, −y)(6.11.16)The Cauchy principal value of RJ has a zero at some value of p < 0, so (6.11.14) will givesome loss of significant figures near the zero.264Chapter 6.Special Functions#include <math.h>#include "nrutil.h"#define ERRTOL 0.08#define TINY 1.5e-38#define BIG 3.0e37#define THIRD (1.0/3.0)#define C1 (1.0/24.0)#define C2 0.1#define C3 (3.0/44.0)#define C4 (1.0/14.0)float rf(float x, float y, float z)Computes Carlson’s elliptic integral of the first kind, RF (x, y, z).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
