elliptical integrals

2.9 Incomplete Elliptic Integrals A. Purpose An integral of the form Z R t;P(t)1=2 dt (1) in which P( t) is a polynomial of the third or fourth de- gree that has no multiple roots, and R is a rational func- tion oftand P(t)1=2, is either elementary, or is an el- liptic integral . It is always possible to express integrals of the form of Eq. (1) linearly in terms of elementary functions and three elliptic integrals of canonical form . These functions are described more completely in [1] and [2]. Several canonical forms have been proposed, but the most widely used are due to Jacobi, Legendre and Carl- son. In each of Eqs. (2){(4) we present rst Jacobi's and then Legendre's form of the canonical elliptic integrals: F(';k) =Zy 0 1t2
1=2 1k2t2
1=2dt =Z' 0 1k2sin2
1=2d(2) E(';k) =Zy 0 1t2
1=2 1k2t2
1=2dt =Z' 0 1k2sin2
1=2d(3) (';;k ) =Zy 0 12t2
1 1t2
1=2 1k2t2
1=2dt =Z' 0 12sin2
1 1k2sin2
1=2d(4) in whichy= sin'. If'is equal to =2, the integrals are said to be complete , otherwise they are incomplete . Carlson's forms of the canonical elliptic integrals are RD(a;b;c ) =3 2Z1 0(t+a)1=2(t+b)1=2(t+c)3=2dt (5) in whichaandbare nonnegative such that a+b >0 andcis positive; if either aorbis zero, the integral is complete , otherwise it is incomplete , RF(a;b;c ) =1 2Z1 0(t+a)1=2(t+b)1=2(t+c)1=2dt (6)in whicha,bandcare nonnegative and at most one of them is zero; if one of a,borcis zero, the integral is complete , otherwise it is incomplete , and RJ(a;b;c;r ) =3 2Z1 0(t+r)1(t+a)1=2(t+b)1=2(t+c)1=2dt (7) in whicha,bandcare nonnegative, and at most one of them is zero, and ris nonzero; if one of a,borcis zero, the integral is complete , otherwise it is incomplete . Notice that RD(a,b,c) =RJ(a,b,c,c). But the neces- sity to compute RD(a,b,c) arises frequently in practice, and a procedure especially tailored to compute RD(a,b, c) is more ecient than computing RJ(a,b,c,c). The functionRC(a,b) =RF(a,b,b) is elementary, but also appears frequently. A procedure is provided to compute RC(a,b): Identify a, bandcsuch thatabc, and assume a<c . Then c3=2RD(a;b;c ) =3 k2sin3'[F(';k)E(';k)] (8) c1=2RF(a;b;c ) =F(';k) sin'(9) c3=2RJ(a;b;c;r ) =3 2sin3'[(';;k )F(';k)] (10) where cos2'=a=c,k2= (cb)=(ca) and2= (cr)=(ca): The subprograms described in this chapter evaluate the canonical forms of incomplete elliptic integrals, using ei- ther the Legendre or the Carlson parameterization. B. Usage B.1 Program Prototype, Single Precision, Leg- endre's Form, E and F REAL PHI, K, F, E INTEGER IERR Assign values to PHI and K. CALL SELEFI (PHI, K, F, E, IERR) B.1.a Argument Denitions PHI [in] Argument, ', of the elliptic integral. Require jPHIj=2: K[in] Modulus, k. RequirejKj1:0: c 1997 Calif. Inst. of Technology, 2010 Math  a la Carte, Inc. June 17, 2010 Incomplete Elliptic Integrals 2.9{1 F[out] F(',k) with'given by PHI and kgiven by K. E[out] E(',k) with'given by PHI and kgiven by K. IERR [out] status indicator: 0 = no errors 1 = Magnitude of argument too large, jPHIj>= 2: 2 = Magnitude of Modulus too large, jKj>1:0: 3 =jPHIj==2 andjKj= 1, F is innite. B.2 Program Prototype, Single Precision, Leg- endre's Form,  REAL PHI, K2, ALPHA2, PI INTEGER IERR Assign values to PHI, K2 and ALPHA2. CALL SELPII (PHI, K2, ALPHA2, PI, IERR) B.2.a Argument Denitions PHI [in] Argument, ', of the elliptic integral. Require jPHIj=2: K2 [in] Square of the modulus, k2. Requirek2sin2' 1:0. See Section E. ALPHA2 [in] Characteristic, 2. Require2sin2' 1:0. See Section E. PI [out] (';;k ), with'given by PHI, 2given by ALPHA2 and k2given by K2. IERR [out] status indicator. If IERR = 0, there were no errors. Other values are produced by procedures SRFVAL and SRJVAL (see Sections B.5 and B.6) which are used in computing ( ;;k ): B.3 Program Prototype, Single Precision, Carl- son's Form, RC REAL X, Y, RC INTEGER IERR Assign values to X and Y. CALL SRCVAL (X, Y, RC, IERR) B.3.a Argument Denitions X, Y [in] Arguments of the elliptic integral. Require X 0, Y6= 0. See Section E. RC [out] The computed value of RC(X, Y). IERR [out] Status indicator: 0 = no errors 1 = X<0:0 or Y = 0.0. 2 = X +jYjtoo small (See Section E). 3 = X orjYjor X +jYjtoo large (See Section E). 4 = Y<0 andjYjtoo large and X too small (See Section E).B.4 Program Prototype, Single Precision, Carl- son's Form, RD REAL X, Y, Z, RD INTEGER IERR Assign values to X, Y and Z. CALL SRDVAL (X, Y, Z, RD, IERR) B.4.a Argument Denitions X, Y, Z [in] Arguments of the elliptic integral. Re- quire X0, Y0, X + Y>0, Z>0. See Section E. RD [out] The computed value of RD(X, Y, Z). IERR [out] Status indicator: 0 = no errors 1 = X<0:0 or Y<0:0 or Z<0:0: 2 = X + Y too small or Z too small (See Section E). 3 = X or Y or Z too large (See Section E). B.5 Program Prototype, Single Precision, Carl- son's Form, RF REAL X, Y, Z, RF INTEGER IERR Assign values to X, Y and Z. CALL SRFVAL (X, Y, Z, RF, IERR) B.5.a Argument Denitions X, Y, Z [in] Arguments of the elliptic integral. Re- quire X0, Y0, Z0, at most one of X, Y or Z equal zero. See Section E. RF [out] The computed value of RF(X, Y, Z). IERR [out] Status indicator: 0 = no errors 1 = X<0:0 or Y<0:0 or Z<0:0: 2 = X + Y or X + Z or Y + Z too small (See Section E). 3 = X or Y or Z too large (See Section E). B.6 Program Prototype, Single Precision, Carl- son's Form, RJ REAL X, Y, Z, R, RJ INTEGER IERR Assign values to X, Y, Z and R. CALL SRJVAL (X, Y, Z, R, RJ, IERR) 2.9{2 Incomplete Elliptic Integrals June 17, 2010 B.6.a Argument Denitions X, Y, Z, R [in] Arguments of the elliptic integral. Re- quire X0, Y0, Z0, at most one of X, Y or Z equal zero, R6= 0. See Section E. RJ [out] The computed value of RJ(X, Y, Z, R). IERR [out] Status indicator: 0 = no errors 1 = X<0:0 or Y<0:0 or Z<0:0 or R = 0:0: 2 = X + Y or X + Z or Y + Z or jRjtoo small (See Section E). 3 = X or Y or Z or jRjtoo large (See Section E). B.7 Modications for Double Precision For double precision usage, change the REAL type state- ments to DOUBLE PRECISION and change the sub- program names SELEFI, SELPII, SRCVAL, SRDVAL, SRFVAL and SRJVAL to DELEFI, DELPII, DRCVAL, DRDVAL, DRFVAL and DRJVAL, respectively. C. Examples and Remarks C.1 Related Functions Logarithms, inverse circular functions and inverse hy- perbolic functions can be expressed in terms of RC, see [9, pp. 163, 186]: (lnx)=(x1) =RC((1 2+1 2x)2;x); x> 0; (sin1x)=x=RC(1x2;1);1x1; (sinh1x)=x=RC(1 +x2;1);1<x<1; (cos1x)=(1x2)1 2=RC(x2;1);0x1; (cosh1x)=(x21)1 2=RC(x2;1); x1; (tan1x)=x=RC(1;1 +x2);1<x<1; (tanh1x)=x=RC(1;1x2);1<x< 1; cot1x=RC(x2;x2+ 1); 0x<1; coth1x=RC(x2;x21); x> 1: The rst seven of these allow computing nearly indeter- minate forms with more accuracy than would be possible using the na ve formulation. Heuman's lambda function [3] is a variant of Legendre's third integral: 1cos2sin2
1=2 cos2sincos(;;' ) = sin'RF(cos2';1sin2sin2';1) +sin2sin3' 3 1cos2sin2
RJ cos2'; 1sin2sin2';1;1sin2sin2' 1cos2sin2 (11) 20(;) = (;;= 2) = sin RF(0;cos2;1)1 3sin2RD(0;cos2;1) RF cos2;1cos2sin2;1
1 3cos2sin3 RF(0;cos2;1)RD(cos2;1cos2sin2;1) (12) The variants of Legendre's integrals used by Bulirsch in [4] and [5] are el1(x;kc) =xRF 1;1 +k2 cx2;1 +x2
; (13) el2(x;kc;a;b) =axR F 1;1 +k2 cx2;1 +x2
+1 3(ba)x3RD 1;1 +k2 cx2;1 +x2
(14) ele3(x;kc;p) =xRF 1;1 +k2 cx2;1 +x2
+1 3(1p)x3RJ 1;1 +k2 cx2;1 +x2;1 +px2
(15) cel(kc;p;a;b ) =aRF(0;k2 c;1) =1 3(bpa)RJ 0;k2 c;1;p
(16) C.2 Which Procedure Should Be Used? Several factors in uence the choice of procedure. If one needs to write a simple program and use it once, one should probably choose the procedure that evaluates the functions in the form most similar to the way the prob- lem is posed. If one needs to write a program that will have substantial use, one should usually prefer SE- LEFI to SRDVAL and SRFVAL, as the former is up to 30 times faster than the latter two. An exception to this rule occurs if one needs to compute RD(a;b;c ) with c<max(a;b), in which case the parameters for SELEFI will be out of range. If accuracy is an issue but speed is not, one may prefer SRDVAL and SRFVAL to SELEFI, at least for computing F( ';k):(See testing in Section D below). SRCVAL is somewhat slower, on an IBM PC/AT with a numeric data processor, than using the equivalent For- tran intrinsic functions. This is no surprise, as most of the intrinsic functions are implemented by hardware. But the inverse hyperbolic functions are not. SRCVAL is roughly the same speed as the procedures in Chap- ter 2.1. As mentioned above, it may be advantageous to use SRCVAL to compute nearly indeterminate forms. SELPII is implemented by using SRJVAL and SRFVAL (see Eqs. (9) and (10) above). Thus, there is no special June 17, 2010 Incomplete Elliptic Integrals 2.9{3 advantage in speed or accuracy to one or the other. The sole criterion is how closely the forms of the functions evaluated directly by the procedures match the forms of the functions the user needs to evaluate. D. Functional Description D.1 Properties of the Functions The rst form given in Eqs. (2){(4) is the Jacobi or alge- braic form. When expressed in this form Eq. (2) is nite for all real and complex y, including1, has a simple pole of order 1 for y=1, and is logarithmically innite fory= 1=2: D.2 Method of Computation The procedure SELEFI is based upon a procedure ELLPI developed by Allan V. Hershey and modied by Alfred H. Morris, described in [6]. The procedure uses se- ries expansions due to DiDonato and Hershey, described in [7]. The procedure SELPII is based upon a proce- dure EPI developed by Alfred H. Morris, described in [5]. It computes ( ';k2;2) using Eqs. (9) and (10), as computed by SRFVAL and SRJVAL. The procedures SRCVAL, SRDVAL, SRFVAL and SRJVAL are based on procedures developed by B. C. Carlson and Elaine M. Notis, described in [8] and [9]. All of the referenced procedures were revised to be consistent with low level modules and naming conventions of MATH77. D.3 Testing The single precision programs for E( ';k), F(';k), RD(a;b;c ) andRF(a;b;c ) were tested on an IBM PC/AT (using IEEE arithmetic) by comparison to dou- ble precision results, as described below. The relative precision of IEEE single precision arithmetic is = 2230:119106: The accuracy of procedure SELEFI was assessed by comparing its results to double precision results ob- tained by applying Eqs. (8) and (9), with RD(a;b;c ) andRF(a;b;c ) evaluated by DRDVAL and DRFVAL, respectively. The accuracy of procedures SRDVAL and SRFVAL was assessed by comparing their results to dou- ble precision results obtained by applying Eqs. (8) and (9), with E( ';k) and F(';k) evaluated by DELEFI. To test SELEFI, the rectangular region 0 '=2 0k1 of the'kplane was divided into 2000 re- gions, and a point was randomly selected in each region. To test SRDVAL and SRFVAL, the argument cwas set to 1.0, the rectangular region 0 a < 10b <1 of theabplane was divided into 2000 regions, and a point was randomly selected in each region. The maxi- mum relative and absolute errors are summarized in the following table.Max. Rel. Max. Abs. Function Error Error E(';k) :82 : 98 F(';k) 5 :24 15:91 RD(a;b;1) 1:20 3:01 RF(a;b;1) 1:35 2:55 Errors in F( ';k) increase as the arguments approach the innite singularity at '==2 andk= 1: References 1. Milton Abramowitz and Irene A. Stegun, Handbook of Mathematical Functions ,Applied Mathematics Series 55 , National Bureau of Standards (1966) Chap- ter 17, 587{626. 2. Paul F. Byrd and Morris D. Friedman, Handbook of Elliptic Integrals for Engineers and Scientists , Springer Verlag, Berlin (1971). 3. H. Kuki, Tables of complete elliptic integrals ,J. Math. and Physics 20 (1941) 127{206. 4. Roland Bulirsch, Numerical calculation of elliptic in- tegrals and elliptic functions ,Numerische Mathe- matik 7 (1965) 78{90. 5. Roland Bulirsch, Numerical calculation of elliptic in- tegrals and elliptic functions ,Numerische Mathe- matik 13 (1969) 305{315. 6. Alfred. H. Morris, Jr., NSWC Library of Mathematics Subroutines . Technical Report NSWCDD/TR-92/425, Naval Surface Warfare Center, Dahlgren, VA 22448-5000 USA (Jan. 1993) 107{110. 7. Armido R. DiDonato and Allan V. Hershey, New for- mulas for computing incomplete elliptic integrals of the rst and second kind ,J. ACM 6 (1959) 512{526. 8. B. C. Carlson, Computing elliptic integrals by dupli- cation ,Numerische Mathematik 33 (1979) 1{16. 9. B. C. Carlson and Elaine M. Notis, Algorithm 577: Algorithms for incomplete elliptic integrals [S21] , ACM Trans. on Math. Software 7 , 3 (Sept. 1981) 398{403. 10. B. C. Carlson, Special Functions of Applied Mathematics , Academic Press, New York (1977). E. Error Procedures and Restrictions The procedure SELEFI requires j'j=2, andjkj1. Procedure SELPII computes ( ';k2;2) from RJ(a;b;c;r ) andRF(a;b;c ) using Eqs. (9) and (10). The initial values for the arguments are a= cos2', b= 1k2sin2',r= 12sin2', andc= max(a;b;r ). Then a,bandrare replaced by ac,bcandrc, respectively. SELPII requires j'j=2. Restrictions onk2and2are enforced indirectly by restrictions on 2.9{4 Incomplete Elliptic Integrals June 17, 2010 a,bandcimposed by SRFVAL and SRJVAL, described below. The ranges for 'andk2can be extended using formu- lae 113.01, 113.02, 114.01, 115.01, 115.02, 160.02, 161.02 and 162.02 from [2], or formulae 17.4.1 through 17.4.18 from [1]. Denote the largest representable magnitude by , and the smallest nonzero representable magnitude by !. General restrictions on the arguments to procedures SR- CVAL, SRDVAL, SRFVAL and SRJVAL were described above in Section B. SRCVAL requires X + jYj5!;X =5,jYj =5, and, if Y<2:236=p!it requires X(! )2=25. Denote the machine round-o level by , that is,is the smallest positive number such that the representa- tion of 1 + is dierent from 1. Let "be the solution of the equation = 3"6(1")3=2, D= 2 2=3and !D="!2=3=10. SRDVAL requires X + Y !Dand Z!D;X D;Y Dand Z D. SRFVAL requires X + Y 5!;X + Z5!, Y + Z 5!;X =5;Y =5 and Z =5. Let J= ( =5)1=3=5 and!J= (5!)1=3. SRJVAL requires X + Y !J;Y + Z!J;X + Z!J, jRj!J;X J;Y J;Z JandjRj J. The accessible ranges of the arguments may be extended beyond the ranges admissible in the procedures by using the homogeneity of the functions: RF(ka;kb;kc ) =k1=2RF(a;b;c );and RJ(ka;kb;kc;kr ) =k3=2RJ(a; b; c; r ): If any of the restrictions above is violated, all proce- dures return an error indicator in the argument named IERR, and invoke the error message processor (see Chap- ter 19.2) with LEVEL = 0. The procedure ERMSET (see Chapter 19.2) may be used to aect the default er- ror processing action.F. Supporting Information The source language for these subroutines is ANSI For- tran 77. The procedures SELEFI and SELPII were written by W. V. Snyder in December 1990, based on earlier pro- cedures described by Alfred H. Morris, Naval Surface Warfare Center, Dahlgren, VA in [5]. The procedures SRCVAL, SRDVAL, SRFVAL and SRJVAL were writ- ten by W. V. Snyder in December 1990, based on earlier procedures described by Carlson and Notis in [9]. Entry Required Files DELEFI AMACH, DELEFI, DERM1, DERV1, DLNREL, ERFIN, ERMSG DELPII AMACH, DELPII, DERM1, DERV1, DRCVAL, DRFVAL, DRJVAL, ERFIN, ERMSG DRCVAL AMACH, DERM1, DERV1, DRCVAL, ERFIN, ERMSG DRDVAL AMACH, DERM1, DERV1, DRDVAL, ERFIN, ERMSG DRFVAL AMACH, DERM1, DERV1, DRFVAL, ERFIN, ERMSG DRJVAL AMACH, DERM1, DERV1, DRCVAL, DRFVAL, DRJVAL, ERFIN, ERMSG SELEFI AMACH, ERFIN, ERMSG, SELEFI, SERM1, SERV1, SLNREL SELPII AMACH, ERFIN, ERMSG, SELPII, SERM1, SERV1, SRCVAL, SRFVAL, SRJVAL SRCVAL AMACH, ERFIN, ERMSG, SERM1, SERV1, SRCVAL SRDVAL AMACH, ERFIN, ERMSG, SERM1, SERV1, SRDVAL SRFVAL AMACH, ERFIN, ERMSG, SERM1, SERV1, SRFVAL SRJVAL AMACH, ERFIN, ERMSG, SERM1, SERV1, SRCVAL, SRFVAL, SRJVAL June 17, 2010 Incomplete Elliptic Integrals 2.9{5 DRSELI program DRSELI c>>19941019 DRSELI Krogh Changes to use M77CON c>>19920309 DRSELI W V Snyder Create separate s i n g l e and double demos . c>>19911004 DRSELI W V Snyder JPL Original code . cS r e p l a c e s "?": DR?ELI ,?RCVAL,?ELEFI,? ELPII ,?RDVAL,?RFVAL,?RJVAL c c Demonstration d r i v e r f o r incomplete e l l i p t i c i n t e g r a l procedures . c real ALPHA2, E, F, K, K2, PHI , PI , R, RC, RD, RF, RJ real SINPHI , T, U, X, Y, Z integer IERR c c Compute arc sine x using ASIN and RC, f o r x = 0.5 c print , ' I d e n t i t i e s from write up : ' x = 0.5 e0 c a l l s r c v a l ( 1 . 0 e0 xx , 1 . 0 e0 , rc , i e r r ) i f( i e r r . eq . 0 ) then t = asin ( x ) xrc print ' ( ' ' ASIN ( 0 . 5 ) 0.5RC(1 0.52 ,1) =' ' , g15 . 8 ) ' , t else print ' ( ' ' SRCVAL returns e r r o r s i g n a l ' ' , i 1 ) ' , i e r r end i f c c Evaluate i d e n t i t i e s given by equations (8 10) in the write up c with k 2 = 1/2 , sin ( phi ) 2 = 1/4 , alpha 2 = 1/2 , c = 1. c From t h i s , we have a = 3/4 , b = r = 7/8. c alpha2 = 0.5 e0 k = s q r t ( 0 . 5 e0 ) k2 = 0.5 e0 s i n p h i = 0.5 e0 phi = asin ( s i n p h i ) r = 0.875 e0 x = 0.75 e0 y = 0.875 e0 z = 1.0 e0 c a l l s e l e f i ( phi , k , f , e , i e r r ) i f( i e r r . ne . 0 ) then print ' ( ' ' SELEFI returns e r r o r s i g n a l ' ' , i 1 ) ' , i e r r go to 99 end i f c a l l s e l p i i ( phi , k2 , alpha2 , pi , i e r r ) i f( i e r r . ne . 0 ) then print ' ( ' ' SELPII returns e r r o r s i g n a l ' ' , i 1 ) ' , i e r r go to 99 end i f c a l l s r d v a l (x , y , z , rd , i e r r ) i f( i e r r . ne . 0 ) then print ' ( ' ' SRDVAL returns e r r o r s i g n a l ' ' , i 1 ) ' , i e r r go to 99 end i f c a l l s r f v a l (x , y , z , rf , i e r r ) i f( i e r r . ne . 0 ) then print ' ( ' ' SRFVAL returns e r r o r s i g n a l ' ' , i 1 ) ' , i e r r go to 99 2.9{6 Incomplete Elliptic Integrals June 17, 2010 end i f c a l l s r j v a l (x , y , z , r , rj , i e r r ) i f( i e r r . ne . 0 ) then print ' ( ' ' SRJVAL returns e r r o r s i g n a l ' ' , i 1 ) ' , i e r r go to 99 end i f u = s q r t ( z 3)rd t = 3.0 e0 / ( k2 s i n p h i 3)( fe ) r = (u t ) / u print ' ( ' ' Equation ( 8 ) , (LHS RHS)/LHS =' ' , g15 . 8 ) ' , r u = s q r t ( z ) r f t = f / s i n p h i r = (u t ) / u print ' ( ' ' Equation ( 9 ) , (LHS RHS)/LHS =' ' , g15 . 8 ) ' , r u = s q r t ( z 3)r j t = 3 / ( alpha2 s i n p h i 3)( pif ) r = (u t ) / u print ' ( ' ' Equation (10) , (LHS RHS)/LHS =' ' , g15 . 8 ) ' , r c 99 stop end ODSELI I d e n t i t i e s from write up : ASIN ( 0 . 5 ) 0.5RC(1 0.52 ,1) = 0.59604645E 07 Equation ( 8 ) , (LHS RHS)/LHS = 0.18963517E 05 Equation ( 9 ) , (LHS RHS)/LHS = 0.0000000 Equation (10) , (LHS RHS)/LHS = 0.23314153E 05 June 17, 2010 Incomplete Elliptic Integrals 2.9{7