Toward Symbolic Integration of Elliptic Integrals

B. C. CARLSON1 2
Ames Laboratory and Department of Mathematics,
Iowa State University, Ames, Iowa 50011-3020

Abstract

A method is proposed by which elliptic integrals can be integrated symbolically without the kind of information about limits of integration and branch points of the integrand that is required in integral tables using Legendre's integrals. However, it is assumed that when all polynomials in the integrand have been factored symbolically into linear factors, the exponents of all distinct linear factors are known. The recurrence relations are one-parameter relations, all formulas are given explicitly, and the integral is eventually expressed in terms of canonical R-functions, with no increase in their number if neither limit of integration is a branch point of the integrand. It is the use of R-functions rather than Legendre's integrals that makes it possible to carry out the whole process symbolically. If (possibly complex) numerical values of the symbols are known, there are published algorithms for numerical computation of the R-functions.

1  Introduction

Elliptic integrals are usually defined to be integrals of the form
ó
õ 		x

y 
 r(t,s) dt ,
(1)
where s2 is a polynomial of degree three or four in t with simple zeros and where r is a rational function of t and s containing at least one odd power of s. In 1954 P. F. Byrd and M. D. Friedman published the first edition of an extensive table of elliptic integrals [ Byrd and Friedman1971] that became an important source. An attempt to generate in a computer-algebra system a substantial part of this table, as well as extend it further, faces several complications, including the following.
(1) The integrals were first converted, by various substitutions that differ from one integrand to another, to integrals of Jacobian elliptic functions, which were then expressed separately in terms of Legendre's canonical forms. Two Russian tables [ Gradshteyn and Ryzhik1994, Prudnikov et al. 1986] later omitted the intermediate integrals of Jacobian elliptic functions and did not specify methods of integration.
(2) Each interval of integration began or ended at a branch point of the integrand. An integral with neither limit of integration at a branch point would have to be split into two parts, doubling the number of canonical forms, and even this remedy was not always available because of divergence at both neighboring branch points.
(3) The ordering of the free limit of integration and the branch points on the real line had to be specified. In conjunction with (2) this often required listing 8, 18, 36, or even 72 cases of what was essentially the same integral. The formulas for these numerous cases often had little resemblance to one another, differing by an integral of the first kind or an algebraic function or both, and gave no hint of how they might be unified. The problem is most easily seen in the arrangement chosen by .
investigated four possible approaches to symbolic integration of elliptic integrals and found serious difficulties in each case. Besides the method of Byrd and Friedman they considered direct conversion to the canonical forms of Legendre and of Weierstrass. In the fourth method, elaborated earlier by , an elliptic integral was represented as a multivariate hypergeometric R-function before being reduced to canonical forms by recurrence relations. Although the R-functions offered certain advantages, a major difficulty lay in multiparameter recurrence relations. settled on a first stage of reduction by a classical method using one-parameter recurrence relations, to be followed by a second stage of expression in terms of R-functions or alternatively Legendre's integrals. Implementation was not free of troubles.
In a paper not yet published, choose Legendre's canonical forms but make improvements on the classical method. They use Hermite reduction [ Moses1971][ Geddes et al. 1992] to arrive at integrands without multiple poles, and the terms with simple poles are decomposed by an implicit full partial-fraction expansion due to . To reduce to Legendre's canonical forms they finally choose one of 21 transformations according to the zeros of s2 and the limits of integration.
In the present paper we separate two stages of reduction, as did . In the first stage we start from an integrand symbolically factored into possibly complex linear factors, permitting explicit formulas for reduction by one-parameter recurrence relations to a set of basic integrals somewhat different from the classical ones. All coefficients can remain symbolic throughout this reduction, which is relatively simple for integrands that occur frequently in practice. For integrands containing polynomials of high degree, like some of the examples chosen by Labahn and Mutrie, the iteration of recurrence relations could become onerous, and it might be better to begin, as they do, with a preliminary Hermite reduction before applying the present method to the resulting individual terms.
In the second stage the basic integrals are expressed in terms of canonical R-functions, using reduction theorems unknown at the time of Ng and Polajnar's work. These make it possible to avoid the complications in (2) above and unify the numerous cases mentioned in (3) (e.g. see [ Carlson1988]) while still retaining symbolic coefficients, some of which may be complex. Although the 21 transformations used by are avoided, and although there are efficient algorithms [ Carlson1995] for computing numerically the canonical R-functions of complex variables, the possible presence of complex numbers will seem a disadvantage to those who use a computer language without complex arithmetic and to those who demand that a real integrand have an integral that contains no complex numbers. The complex numbers can be eliminated by a Landen transformation (e.g. see [ Carlson1991]), but only at the cost of more complicated formulas. On the other hand an integrand containing only linear factors with symbolic coefficients can be integrated in terms of canonical R-functions with no assumptions about the numerical values of the symbols; for this purpose Legendre's canonical forms are unsuitable.
A simple rationalization procedure, found in every introduction to the subject, allows us to separate from (1.1) the integral of a rational function and assume that r(s,t) = r(t)/s, where r is a rational function of t. (Moreover, if s2 is an even function of t, it is advantageous but not essential to separate the odd part of r, which leads to an elementary integral, and replace t by Öt in the remaining elliptic integral.) Instead of (1.1) we henceforth consider
ó
õ 		x

y 
 P(t) dt
Q(t)s(t)
 = ó
õ 		x

y 
 h
Õ
i=1 
 (ai +bi t)-1/2 n
Õ
j=1 
 (aj +bj t)mj dt ,
(2)
where P and Q are polynomials, h=3 or 4, and the m's are integers. We assume that the integral is well defined (in particular that the open interval of integration contains no finite branch point of the integrand) and that no two of the n linear factors on the right side are proportional. Although the a's and b's, some of which may be complex, can be left as symbols throughout the integration procedure described in this paper, the integer values of the m's are needed at the beginning. In many cases, including most of those in present integral tables, the integrand contains contains only linear and quadratic factors, making the m's obvious. If unfactored polynomials of higher degree are present, we must not only ensure that P and Q have no common divisor (except a constant) but also determine whether they have divisors in common with s2. If so, and if all coefficients are given numerically, the m's can be determined by making square-free factorizations (, § 8.2) of P and Q and finding the greatest common divisor of each square-free factor and s2, which is square-free by assumption.

2  Partial-fraction decomposition

We begin by decomposing into partial fractions the rational function P/Q in the integrand of (1.2), performing this task analytically once and for all rather than requiring the computer to do it in each individual case.
Lemma 1 If mj is an integer and |rzj| < 1 for j = 1,¼,n , then
n
Õ
j=1 
 (1-rzj)mj = ¥
å
s=0 
 Ts rs ,
(3)

T0 = 1 ,       Ts = å
 H1a1¼Hsas
a1 ! ¼as !
  ,       Hs = -1
s
   n
å
j=1 
  mj zjs ,        s ³ 1,
(4)
where the sum defining Ts extends over all nonnegative integers a1 ,¼,as such that a1 + 2a2 + ¼+ sas = s. If mj Î {0,1} for j=1,¼,n then (-1)s Ts is by definition the elementary symmetric function Es of degree s in m1 z1 ,¼,mn zn .
Using the power-series expansion of a logarithm, we see that
ln n
Õ
j=1 
 (1-rzj)mj = n
å
j=1 
 mj ln(1-rzj) = - n
å
j=1 
 mj  ¥
å
s=1 
 rs zjs
s
 = ¥
å
s=1 
 Hs rs .
(5)
To obtain the coefficient Ts of rs in
n
Õ
j=1 
 (1-rzj)mj = exp ¥
å
s=1 
 Hs rs = ¥
Õ
k=1 
 exp(Hk rk) ,
(6)
we omit factors in the infinite product with k > s, multiply the exponential series for exp(H1 r),¼,exp(Hs rs), and pick out the coefficient of rs.
Definition 1 Define
M= n
å
j=1 
 mj ,       B = n
Õ
j=1 
 bjmj,       dij = ai bj -aj bi ,       Di = n
Õ
j ¹ ij=1 
 dji mj ,
(7)

m±s(i) = -1
s
 n
å
j ¹ ij=1 
 mj æ
è 		dij
bj
 ö
ø 		±s

 
  ,       s = 1,2,3,¼,
(8)

C0(i)=1,       C±s(i) = å
 m±1a1(i)¼ m±sas(i)
a1 ! ¼as !
  ,       s = 1,2,3,¼,
(9)
where upper (lower) signs go together and the last sum extends over all nonnegative integers a1,¼,as such that a1 +2a2 +¼+ sas = s. In particular we have
C±1(i)
=
m±1(i) ,       C±2(i) = 1
2
   m2±1(i)+m±2(i) ,      
C±3(i)
=
1
6
  m3±1(i) + m±1(i) m±2(i) + m±3(i) .      
(10)
C±4(i)
=
1
24
  m4±1(i) + 1
2
  m2±1(i) m±2(i) + m±1(i)m±3(i) + 1
2
  m2±2(i) +m±4(i) .
 We note that dij ¹ 0 if i ¹ j because of the assumption that no two linear factors on the right side of (1.2) are proportional.
Theorem 1 The decomposition of P(t)/Q(t) into partial fractions is
P(t)
Q(t)
 = n
Õ
j=1 
 (aj + bj t)mj
=
B bi-M  M
å
q=0 
 CM-q(i) (ai +bi t)q
+
n
å
i=1 
 Di bimi -M -mi
å
q=1 
  Cmi +q(i) (ai+bi t)-q .   
(11)
The polynomial part is independent of the choice of i, and each sum over q is empty if the upper limit is less than the lower limit.
The definition (2.5) of dij implies
aj +bj t = bj
bi
  (ai + bi t) æ
è 		1 - dij
bj(ai +bi t)
 ö
ø 		 ,
(12)
whence
n
Õ
j=1 
 (aj + bj t)mj = B bi-M (ai + bi t)M  n
Õ
j=1 
   æ
è 		1 - dij
bj(ai +bi t)
 ö
ø 		mj

 
 .
(13)
This equation is independent of the choice of i. If |t| is sufficiently large, we can expand the last product by Lemma 2.1 with r = 1/(ai +bi t) and zj = dij/bj . It follows that Hs is m+s(i) as defined by (2.6), the term in Hs with j=i being irrelevant because dii=0. Thus we see that Ts is C+s(i) as defined by (2.7), and (2.11) becomes
P(t)
Q(t)
 = B bi-M ¥
å
s=0 
 Cs(i)(ai +bit)M-s.
(14)
Since this is valid for sufficiently large t, the polynomial part of P/Q consists of the terms with 0 £ s £ M, which are absent unless M ³ 0. By putting s=M-q these become the first term on the right side of (2.9).
To find the rest of (2.9), which contains poles arising from negative m's, we replace (2.10) by
aj + bj t = dji
bi
 æ
è 		1- bj
dij
   (ai +bi t) ö
ø 		 ,    j ¹ i ,
(15)
whence
n
Õ
j=1 
 (aj +bj t)mj = Di bimi -M (ai +bi t)mi  n
Õ
j ¹ ij=1 
   æ
è 		1- bj
dij
  (ai +bi t) ö
ø 		mj

 
 .
(16)
If |ai +bi t| is sufficiently small, we can expand the last product by Lemma 2.1 with r = ai +bi t and zj = bj/dij. It follows that Hs is m-s(i) as defined by (2.6) and that Ts is C-s(i) as defined by (2.7). Hence (2.11) becomes
P(t)
Q(t)
 = Di bimi -M  ¥
å
s=0 
 C-s(i) (ai +bi t)mi + s .
(17)
Since this is valid if |ai +bi t| is sufficiently small, the part of P/Q representing a pole at -ai/bi consists of the terms with 0 £ s £ -mi -1, which are absent unless -mi ³ 1. These can be rewritten by putting s = -mi -q as
Di bimi -M  -mi
å
q=1 
 Cmi +q(i) (ai +bi t)-q.
(18)
Summation over i takes account of all the poles of P/Q and yields the second term on the right side of (2.9).
We shall now apply the partial-fraction decomposition (2.9) to the integral (1.2), designated henceforth by
I(m) = ó
õ 		x

y 
 h
Õ
i=1 
 (ai +bi t)-1/2 n
Õ
j=1 
 (aj +bj t)mj dt ,
(19)
where m=(m1,¼,mn) is an n-tuple of integers. We define n-tuples
e0 = (0,0,¼,0),
e1 = (1,0,¼,0),
(20)
:
en = (0,¼,0,1).
Substitution of (2.9) in (2.17) reduces I(m) to a set of simpler integrals in which at most one component of m is nonzero:
I(m) = B bi-M  M
å
q=0 
 CM-q(i) I(qei) + n
å
i=1 
  Di  bimi -M  -mi
å
q=1 
  Cmi +q(i) I(-qei) .
(21)
The first term on the right side is independent of the choice of i, and each sum over q is empty if the upper limit is less than the lower limit.

3  Recurrence relations and basic integrals

The integrals I(±qei) in (2.19) can be expressed in terms of I(e0) and I(±ei) by using recurrence relations. The main relation involves also an algebraic term of the form A(m)=vm(x)-vm(y), where vm(t) has the form of an integrand in (2.17).
Theorem 2 Define
A(m)=vm(x)-vm(y),
(22)

vm(t)= 1
s(t)
   n
Õ
j=1 
 (aj +bj t)mj,       s(t)= h
Õ
i=1 
   _______
Öai +bi t
 
  .
(23)
For 1 £ j £ n and h=3 or 4 let Ep(j) be the elementary symmetric function of degree p in
d1j
b1
  ,¼, dhj
bh
  ,
(24)
and define
s0(j) = s0 = h
Õ
i=1 
 bi ,       sp(j)=s0 Ep(j) ,
(25)
whence sh(j)=0 if 1 £ j £ h. Then, for any integer q and for 1 £ j £ n,
h
å
r=0 
 æ
è 		q+ r
2
 +1 ö
ø 		 sh-r(j) I((q+r)ej) = bjh-1 A æ
è 		(q+1)ej + h
å
i=1 
 ei ö
ø 		.
(26)

Choose m=(q+1)ej +åhi=1ei , whence
vm(t) = s(t)(aj +bj t)q+1,
(27)

v¢m(t) = s¢(t)(aj +bj t)q+1 +(q+1) bj s(t)(aj +bj t)q .
(28)
From
ai +bi t = bi
bj
   æ
è 		aj +bj t + dij
bi
 ö
ø 		,
it follows that
s2(t)
=
s0
bjh
   h
Õ
i=1 
 æ
è 		aj +bj t + dij
bi
 ö
ø
=
s0
bjh
   h
å
r=0 
 Eh-r(j)(aj +bj t)r
=
1
bjh
   h
å
r=0 
 sh-r(j)(aj +bj t)r
(29)
and therefore
s¢(t)= 1
2bjh-1s(t)
   h
å
r=0 
 rsh-r(j)(aj +bj t)r-1.
(30)
Replacing s by s2 /s in the last term of (3.7) and substituting (3.8) and (3.9), we find
v¢m(t)= 1
bjh-1s(t)
   h
å
r=0 
 æ
è 		r
2
   +q+1 ö
ø 		sh-r(j) (aj +bj t)q+r.
(31)
Integration from y to x proves (3.5).
If 1 £ j £ h one of the quantities listed in (3.3) is 0. Then Eh(j)=sh(j)=0, and the sum over r in (3.5) effectively starts with r=1. 3 Thus the recurrence relation contains h integrals if 1 £ j £ h or h+1 integrals if h+1 £ j £ n. In the first case I(±qej) can be reduced to the h-1 integrals I(-ej), I(e0),¼,I((h-3)ej), but in the second case I((h-2)ej) is added to the list. Further reduction uses the following theorem in the cases q=1 and 2, cases that could alternatively be obtained less directly from (2.19).
Theorem 3 If q is a positive integer, then
biq I(qej) = q
å
r=0 
 (
q
r
 )  bjr djiq-r I(rei).
(32)
 Raising both sides of
bi (aj +bj t) = bj (ai +bi t) + dji
(33)
to the power q and making a binomial expansion of the right side, we find
biq (aj +bj t)q = q
å
r=0 
 (
q
r
 )  bjr (ai +bi t)r djiq-r .
(34)
Multiplying both sides by Õhk=1 (ak +bk t)-1/2 and integrating from y to x, we obtain (3.11).
If we choose i £ h in the first sum in (2.19), the integrals in that sum and in the first h terms of the second sum can be reduced by (3.5) to integrals of the form I(qei),  i £ h, with q=-1,0 if h=3 and q=-1,0,1 if h=4. The remaining terms of the second sum can be reduced by (3.5) to integrals of the form I(qej),  j > h, with q=-1,0,1 if h=3 and q=-1,0,1,2 if h=4. Now I(qej),  q=1,2, can be expressed by (3.11) in terms of I(rei) with i £ h and 0 £ r £ q, and these integrals can be treated in the same way as the integrals in the first sum in (2.19). In this way all the integrals on the right side of (2.19) can be reduced to basic integrals of the form I(-ej),  0 £ j £ n, in the cubic case (h=3) plus I(ei),  1 £ i £ 4 in the quartic case (h=4).
As the integrand in (1.2) becomes more complicated, so of course does this reduction to basic integrals, but for integrands that occur frequently in practice it is quick and simple. To illustrate this, we show how 136 quartic integrals from [ Gradshteyn and Ryzhik1994] are represented in terms of basic integrals by five formulas. Each entry in Table 1 starts with a list of those integrals that have the form I(m) = I(åj=15 mj ej) displayed in the entry and defined in (1.2). Also displayed is the expression for I(m) in terms of the basic integrals I(e0),  I(-e5), and I(±ei),  i=1,¼,4. Equation (2.19) suffices for I(e5) and I(ek -ei) ; recurrence relations are not needed.
Table 1: Some integrals from Gradshteyn and Rhyzhik in terms of basic integrals
§ 3.147, 1-8: I(e0)
§ 3.148, 1-8: I(e5)=[1/(bi)] [b5 I(ei)-di5 I(e0)]
§ 3.149, 1-8; § 3.151, 1-8]: I(-e5)
§ 3.167, 1-32: I(ei)
§ 3.168, 1-72: I(ek -ei) = [1/(bi)] [dki I(-ei) + bk I(e0)]
If the basic integrals are to be expressed in terms of Legendre's canonical forms, each of 136 formulas has to be accompanied by inequalities relating the branch points of the integrand and the interval of integration. In the next Section we shall use R-functions to avoid most of this complexity.

4  Reduction of basic integrals to R-functions

Let z1 , z2 , z3 lie in the complex plane cut along the negative real axis, at most one of them being 0. We define two elliptic integrals that are symmetric in z1 , z2 , z3 :
RF(z1 , z2 , z3) = 1
2
   ó
õ 		¥

0 
 dt
3
Õ
i=1 
   ____
Öt+zi
 
  ,
(35)

RJ(z1 , z2 , z3 , zn) = 3
2
   ó
õ 		¥

0 
 dt
3
Õ
i=1 
   ____
Öt+zi
 
   (t+zn)
  ,    zn ¹ 0 .
(36)
The functions RF and RJ are respectively homogeneous of degree -1/2 and -3/2 in their variables. If zn is real and negative, RJ takes its Cauchy principal value. Useful degenerate cases are
RC(z1 , z2) = RF(z1 , z2 , z2) = 1
2
   ó
õ 		¥

0 
 dt
  ____
Öt+z1
 
  (t+z2)
  ,    z2 ¹ 0 ,
(37)

RD(z1 , z2 , z3) = RJ(z1 , z2 , z3 , z3)= 3
2
   ó
õ 		¥

0 
 dt
2
Õ
i=1 
   ____
Öt+zi
 
   (t+z3)3/2
  ,    z3 ¹ 0 .
(38)
If z2 is real and negative, RC takes its Cauchy principal value,
RC(z1 , -p) =   æ
Ö
z1
z1 +p
 
  RC(z1 +p, p) ,       p > 0 .
(39)
If the z's are distinct, then RF , RD , and RJ are elliptic integrals of the first, second, and third kinds, respectively, and RC is an inverse circular or inverse hyperbolic function.
The following reduction theorem, which is not part of the classical theory of elliptic integrals, replaces a quartic polynomial under the square root by a cubic polynomial in the integrands of RF and RJ. Permuting the zeros of the quartic polynomial simply permutes the zeros of the cubic.
Theorem 4 Let z1 , z2 , z3 , z4 lie in the plane cut along the negative real axis, at most one of them being 0, and take their square roots in the right half-plane. Let {i,j,k,l}={1,2,3,4} , and define
zij=zi -zj ,       Vij=   __
Özi
 
   __
Özj
 
 +   __
Özk
 
   __
Özl
 
  ,
(40)
whence
V2ij - V2ik = zilzjk .
(41)
Assume that Vij has positive real part or is 0. (At most one of the V's is 0 because the sum of any two is not 0.) Then
ó
õ 		¥

0 
 dt
4
Õ
j=1 
   ____
Öt+zj
 
 = 2 RF(V212 , V213 , V223) .
(42)
Let zn be complex and nonzero, and if zn ¹ zi define
zin=zi -zn ,       V2in=V2ij- zikzilzjn
zin
  ,       ain = zn+
  __
Özj
 
   __
Özk
 
   __
Özl
 
  __
Özi
  ,
(43)
where ain is assumed to have positive real part or be 0. Then
ó
õ 		¥

0 
 t+zi
t+zn
   dt
4
Õ
j=1 
   ____
Öt+zj
 
=
2zijzikzil
3zin
  RJ(V212 ,  V213 , V223 , V2in)
      + 2  RC(a2in , b2in) ,       zn ¹ zi ,
(44)

b2in= zn
zi
  V2in ,              a2in = b2in + zjnzknzln
zin
  .
(45)
We omit the proof of this theorem because a very simple proof of (4.8) is given in [ Carlson1998] and the proof of (4.10) is somewhat cumbersome. It splits the integral into two parts and recombines them by the addition theorem for RJ, which contains a term in RC. One can show that (4.10) reduces to (4.8) in the limit as zn® zi .
Because Vij=Vji=Vkl=Vlk (in particular, V23=V14), there are only three distinct V's with subscripts 1,2,3,4.
Corollary 1 Let x > y and {i,j,k,l}={1,2,3,4} . Assume the line segment with endpoints ai +bi x = Xi2 and ai +bi y = Yi2 lies in the complex plane cut along the negative real axis, and take all square roots in the right half-plane. Define
Uij = Xi Xj Yk Yl + Yi Yj Xk Xl
x-y
  ,
(46)
whence
U2ij -U2ik = dil djk ,       dij = ai bj -aj bi .
(47)
Assume that dij ¹ 0 if i ¹ j and that Uij has positive real part or is 0. (At most one of the U's is 0 because of (4.13).) Then
I(e0) = 2 RF(U212 , U213 , U223) ,
(48)
where I(e0) is defined by (2.17) and (2.18).
Let 5 £ n £ n and assume the line segment with endpoints an+bnx = X2n and an+bny = Y2n lies in the complex plane punctured at 0. Define
U2in=U2ij - dik dil djn
din
  ,       din=ai bn-anbi ¹ 0 ,
(49)

Sin= 1
x-y
 æ
è 		Xj Xk Xl
Xi
  Y2n + Yj Yk Yl
Yi
  X2n ö
ø 		,
(50)
and assume Sin has positive real part or is 0. Then
I(ei -en)= 2dij dik dil
3din
  RJ(U212 ,  U213 , U223 , U2in) +2 RC(S2in , Q2in) ,       din ¹ 0 ,
(51)

Q2in= X2nY2n
Xi2 Yi2
  U2in ,          S2in-Q2in= djn dkn dln
din
  .
(52)
 Map the interval of integration in
I(e0)= ó
õ 		x

y 
 dt
4
Õ
j=1 
   _______
Öaj +bj t
 
(53)
onto the positive real line by substituting
t= xs+y
s+1
  ,          s= t-y
x-t
  ,          dt
ds
 = x-y
(s+1)2
  ,
(54)

aj +bj t=Xj2  s+zj
s+1
  ,       zj = Yj2
Xj2
   ,       zij = zi -zj = (x-y)dij
Xi2 Xj2
  .
(55)
The result is
I(e0)
=
N ó
õ 		¥

0 
 ds
4
Õ
j=1 
   ____
Ös+zj
 
  ,       N= x-y
4
Õ
j=1 
 Xj
  ,
=
2N RF(V212 , V213 , V223) ,
=
2 RF(U212 , U213 , U223) ,       Uij = Vij
N
  ,
(56)
where we have used (4.8) and the homogeneity of RF. The last line proves (4.14), while (4.6) and (4.7) become (4.12) and (4.13).
By the same change of integration variable we find
I(ei -en)
=
ó
õ 		x

y 
 ai +bi t
an+bnt
   dt
4
Õ
j=1 
   _______
Öaj +bj t
 
=
L ó
õ 		¥

0 
 s+zi
s+zn
   dt
4
Õ
j=1 
   ____
Ös+zj
 
  ,       L= X2i
X2n
  N ,
=
L é
ë 		2zijzikzil
3zin
  RJ(V212 , V213 ,  V223 , V2in)+ 2 RC(a2in , b2in) ù
û
=
L é
ë 		2X2nN2 dijdikdil
3X2i din
   RJ(V212 , V213 , V223 , V2in) + 2 RC(a2in , b2in) ù
û 		 
=
2dijdikdil
3din
   RJ(U212 , U213 , U223 , U2in) + 2 RC(S2in , Q2in) ,
(57)

Sin = ain
L
  ,             Qin = bin
L
  ,
(58)
where we have used (4.10) and the homogeneity of RJ and RC. The last line of (4.23) proves (4.17), while (4.9) and (4.11) become (4.15), (4.16), and (4.18). To obtain I(-en) from (4.17), we shall want the relation
dkl I(m) = bl I(m+ek) - bk I(m+el) .
(59)
This is proved by multiplying both sides of
dkl = bl (ak +bk t) - bk (al +bl t)
by the integrand of (2.17) and integrating from y to x.
Corollary 2 Under the same conditions as in Corollary 4.1 we have
I(e0)
=
2 RF(U212 , U213 , U223) ,
(60)
din I(-en)
=
2bn
3din
  dij dik dil  RJ(U212 , U213 , U223 , U2in) + 2bn RC(S2in , Q2in)
                                                   -bi I(e0) ,
(61)
dji I(-ei)
=
2
3
  bi djk djl  RD(U2ik , U2jk , U2ij) + 2bi Xj Yj
Xi Yi Uij
 - bj I(e0) ,
(62)
bi I(ei)
=
2
3
  dij dik dli  RJ(U212 , U213 , U223 , U2i0) + 2bi  RC(S2i0 , Q2i0) ,
(63)
bj I(ei)
=
bi I(ej) + dij I(e0) .
(64)
In (4.29) the quantities Ui0 , Si0 , Qi0 are obtained from (4.15), (4.16), and (4.18) by putting an=1 and bn=0, whence Xn=Yn=1 and din=-bi .
For convenient reference, (4.14) is copied here as (4.26).
To prove (4.27) we substitute (4.17) in a special case of (4.25) :
din I(-en) = bn I(ei -en) - bi I(e0) .
(65)
To prove (4.28) we put n =j in (4.27) and note that S2ij-Q2ij=0 implies
RC(S2ij , Q2ij) = RC(Q2ij , Q2ij) = 1
Qij
 = Xi Yi
Xj Yj Uij
  .
Then interchange i and j .
To prove (4.29) from (4.17) we replace n by 0, where a0 = 1 and b0 = 0, whence en=e0 and din=-bi .
Equation (4.30) is a special case of (4.25). The free index i in (4.27) allows four choices for U2in , at most two of which can be real and negative. Thus i can be chosen to avoid a Cauchy principal value of RJ . The free index j in (4.30) achieves the same end in conjunction with (4.29). The free index j in (4.28) can be chosen so that Uij ¹ 0, for at most one of the three U's can be 0 because of (4.13).
In the cubic case (h=3) we put m4 = 0, a4 = 1, and b4 = 0, whence X4 = Y4 = 1 and di4=-bi .

5  Summary of procedure

The parts of this paper that are needed for a symbolic integration program are (1) the reduction formula (2.19), in which various quantities are defined by (2.17), (2.18), and Definition 2.1; (2) the recurrence formulas in Theorems 3.1 and 3.2, along with the instructions in the first paragraph after the proof of Theorem 3.2; (3) the reductions to R-functions in Corollary 4.2, in which various quantities are defined by (4.12), (4.13), (4.15), (4.16), and (4.18), along with the remarks in the last two paragraphs of Section 4; if x or y is infinite, appropriate limits should be taken in (4.12), (4.16), and (4.18). The procedure presented here has not yet been implemented in software. It differs from the method used to construct integral tables like [ Carlson1988] in several ways: it does not depend on human judgment in using recurrence relations, some of the basic integrals are chosen differently, and provisions are made for avoiding Cauchy principal values of RJ and infinite values of RD.

6  An example

The examples in Table 1 are expressed at once in terms of R-functions by (4.26)-(4.30). A more difficult example, not found in [ Gradshteyn and Ryzhik1994], which illustrates the various steps in the first paragraph after the proof of Theorem 3.2, is a cubic case,
I(m)= ó
õ 		x

y 
   æ
Ö
(a1 +b1 t)(a2 +b2 t)
a3 +b3 t
 
   dt
(a5 +b5 t)2
  ,    m=(1,1,0,0,-2).
(66)
We use the subscript 5 and take m4=0 so that we can put a4 = 1 and b4 = 0 when the basic cubic integrals I(-ej), 0 £ j £ 5, are expressed in terms of R-functions.
Applying (2.19) and Definition 2.1, we find
b52 I(m)=b1 b2 I(e0)+(b1 d25 +b2 d15) I(-e5) +d15  d25 I(-2e5) ,
(67)
where the last term needs to be reduced to basic integrals. Putting h=3 and q=-2 in (3.5), we have
2s3(5) I(-2e5) = -s2(5) I(-e5) +s0 I(e5) -2b52 A(-e5 + 3
å
i=1 
 ei  ).
(68)
To deal with I(e5) we use (3.11) to get
bi I(e5) = b5 I(ei) -di5 I(e0) ,       1 £ i £ 3,
(69)
and I(ei) is reduced to basic integrals by putting h=3 and q=-2 in (3.5) . Because s3(i)=0 we find
s0 I(ei)=s2(i) I(-ei)+2bi2 A(-ei + 3
å
i=1 
 ei),       1 £ i £ 3.
(70)
Substituting (6.5) in (6.4), then (6.4) in (6.3), then (6.3) in (6.2), and collecting terms, we arrive at
2b52 I(m) = g0 I(e0)+gi I(-ei)+g5 I(-e5)+ga ,
(71)
where
g0
=
2b1 b2 - s0 d15 d25 di5
bis3(5)
 =2b1 b2 - b1 b2 b3 di5
bi d35
  ,
(72)
gi
=
b5 d15 d25 s2(i)
bi s3(5)
 = b5 dji  dki
d35
  ,    {i,j,k}={1,2,3},
(73)
g5
=
2b1 d25+2b2 d15- d15 d25 s2(5)
s3(5)
 =b1 d25+b2 d15- b3 d15 d25
d35
  ,
(74)
ga
=
2b5 d15 d25
s3(5)
  [bi A(-ei + 3
å
i=1 
 ei)-b5 A(-e5 + 3
å
i=3 
 ei)]
=
2b5
d35
  [bi A(-ei + 3
å
i=1 
 ei) -b5 A(-e5 + 3
å
i=3 
 ei)],
(75)
where we have substituted the values of the s's from (3.4). The next step, simplification of the g's, is perhaps the least straightforward to automate. A glance at (6.7) suggests choosing i=3. The last two terms in (6.9) can be combined by the identity
bi djk + bj dki + bk dij = 0 .
(76)
The last two terms in square brackets in (6.10) can be combined by a companion of (4.25),
dkl A(m) = bl A(m+ek) - bk A(m+el) ,
(77)
whence
ga= 2b5 d5i
d35
  A(-ei -e5 + 3
å
i=1 
 ei) = -2b5 A(e1 +e2 -e5) .
(78)
Thus (6.6) becomes
2b52 I(m)
=
b1 b2 I(e0)+ b5 d13 d23
d35
  I(-e3)
+
æ
è 		b1 d25- b5 d15 d23
d35
 ö
ø 		I(-e5) -2b5 A(e1 +e2 -e5) .
(79)
The basic integrals in this equation are expressed in terms of R-functions by (4.26), (4.28) with i=3 and j chosen as 1, 2, or 4 so that Uij ¹ 0, and (4.27) with n =5 and i chosen as 1, 2, 3, or 4 so that U2i5 is not real and negative. Because we are putting a4 = 1 and b4 = 0, the choice to try first in each case is 4 in order to give the extra term in I(e0) a coefficient 0. In the other coefficients and in the variables of the R-functions, defined by (4.12), (4.15), and (4.18), we note that X4 = Y4 = 1, dj4=-bj, and d45=b5.
It is a matter of taste whether or not to substitute the expressions in terms of R-functions for the basic integrals in (6.14). If one wants to exhibit a result with a single equality sign, the right side becomes unwieldy. For an integral table one would list the basic integrals separately, as we have done here.
As a check we have computed
I(m)= ó
õ 		2

0.5 
   æ
Ö
(0.3 +0.3 t)(0.5 +0.1 t)
0.7 -0.1 t
 
   dt
(0.9 -0.3 t)2
 = 6.2430  9544 ,
(80)
obtained from (6.6) with
                     g0 = 0.03
          I(e0)
=
3.0973  71530
                     g3 = 0.072
          I(-e3)
=
5.3785  84194
                     g5 = -0.18
          I(-e5)
=
6.1542  46948

ga = -2b5 A(e1 +e2 -e5) = -2b5 æ
è 		X1 X2
X3 X52
 - Y1 Y2
Y3 Y52
 ö
ø 		= 1.7513  42421 .
A numerical integration gave I(m)=6.2430  97 with a stated limit of absolute error of 0.0000 4, much larger than the difference from (6.15).
Acknowledgment: I wish to thank George Labahn for sending me a preprint of [ Labahn and Mutrie1997].

References

[ Bronstein and Salvy1993]
Bronstein, M., Salvy, B. (1993). Full partial fraction decomposition of rational functions. In Proc. of ISSAC 93 (ISSAC'93) , pages 157-160, New York. The Assoc. for Computing Machinery.
[ Byrd and Friedman1971]
Byrd, P. F., Friedman, M. D. (1971). Handbook of Elliptic Integrals for Engineers and Physicists . Springer-Verlag, New York, 2nd edition.
[ Carlson1988]
Carlson, B. C. (1988). A table of elliptic integrals of the third kind. Math. Comp. , 51 :267-280, S1-S5.
[ Carlson1991]
Carlson, B. C. (1991). A table of elliptic integrals: one quadratic factor. Math. Comp. , 56 :267-280.
[ Carlson1995]
Carlson, B. C. (1995). Numerical computation of real or complex elliptic integrals. Numer. Algorithms , 10 :13-26.
[ Carlson1998]
Carlson, B. C. (1998). Elliptic integrals: symmetry and symbolic integration. Atti dei Convegni Lincei, Accad. Naz. Lincei . Preprint of a paper presented Dec. 1, 1997, in Turin at the conference on Tricomi's Ideas and Contemporary Applied Mathematics.
[ Geddes et al. 1992]
Geddes, K. O., Czapor, S. R., Labahn, G. (1992). Algorithms for Computer Algebra . Kluwer, Boston.
[ Gradshteyn and Ryzhik1994]
Gradshteyn, I. S., Ryzhik, I. M. (1994). Table of Integrals, Series, and Products . Academic Press, New York, 5th edition.
[ Labahn and Mutrie1997]
Labahn, G., Mutrie, M. (1997). Reduction of elliptic integrals to Legendre normal form. Preprint.
[ Moses1971]
Moses, J. (1971). Symbolic integration: the stormy decade. Comm. ACM , 14 :548-560.
[ Ng1974]
Ng, E. W. (1974). Symbolic integration of a class of algebraic functions. Technical Memorandum 33-713, Jet Propulsion Laboratory.
[ Ng and Polajnar1976]
Ng, E. W., Polajnar, D. (1976). A study of alternative methods for the symbolic calculation of elliptic integrals. In Jenks, R. D., editor, Proc. 1976 ACM Symp. on Symbolic and Algebraic Computation (SYMSAC'76) , pages 372-376, New York. The Assoc. for Computing Machinery.
[ Prudnikov et al. 1986]
Prudnikov, A. P., Brychkov, Y. A., Marichev, O. I. (1986). Integrals and Series , volume 1. Gordon and Breach, New York.

Footnotes:

1This work was supported by the Director of Energy Research, Office of Basic Energy Sciences. The Ames Laboratory is operated for the U. S. Department of Energy by Iowa State University under Contract W-7405-ENG-82.
2E-mail: carlson@decst2.ams.ameslab.gov
3Replacement of r by r+1 and q by q-1 yields Equation (5.12) in [ Carlson1998], an equation that is now seen to be merely a special case of Equation (5.13) and so can be omitted.

File translated from TEX by TTH, version 3.85.
On 26 May 2009, 10:19.