Toward Symbolic Integration of Elliptic Integrals

B. C. CARLSON $^{1}$ $^{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

\int_{y}^{x} r (t, s) dt,

(1)

where

s^{2}

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

s^{2}

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) = ρ (t) / s

, where

ρ

is a rational function of t. (Moreover, if

s^{2}

is an even function of t, it is advantageous but not essential to separate the odd part of

ρ

, which leads to an elementary integral, and replace

t

\sqrt{t}

in the remaining elliptic integral.) Instead of (1.1) we henceforth consider

\int_{y}^{x} \frac{P (t) dt}{Q (t) s (t)} = \int_{y}^{x} Π_{i = 1}^{h} (a_{i} + b_{i} t)^{- 1 / 2} Π_{j = 1}^{n} (a_{j} + b_{j} t)^{m_{j}} 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

s^{2}

. 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

s^{2}

, 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 $m_{j}$ is an integer and $‖ {rz}_{j} ‖ < 1$ for $j = 1, \dots, n$ , then

Π_{j = 1}^{n} (1 - {rz}_{j})^{m_{j}} = \sum_{s = 0}^{\infty} T_{s} r^{s},

(3)

T_{0} = 1, T_{s} = \sum \frac{H_{1}^{α_{1}} \dots H_{s}^{α_{s}}}{α_{1}! \dots α_{s}!}, H_{s} = \frac{- 1}{s} \sum_{j = 1}^{n} m_{j} z_{j}^{s}, s \geq 1,

(4)

where the sum defining

T_{s}

extends over all nonnegative integers

α_{1}, \dots, α_{s}

such that

α_{1} + 2 α_{2} + \dots + s α_{s} = s

. If

m_{j} \in {0, 1}

for

j = 1, \dots, n

then

(- 1)^{s} T_{s}

is by definition the elementary symmetric function

E_{s}

of degree

s

m_{1} z_{1}, \dots, m_{n} z_{n}

Using the power-series expansion of a logarithm, we see that

\ln Π_{j = 1}^{n} (1 - {rz}_{j})^{m_{j}} = \sum_{j = 1}^{n} m_{j} \ln (1 - {rz}_{j}) = - \sum_{j = 1}^{n} m_{j} \sum_{s = 1}^{\infty} \frac{r^{s} z_{j}^{s}}{s} = \sum_{s = 1}^{\infty} H_{s} r^{s} .

(5)

To obtain the coefficient

T_{s}

r^{s}

Π_{j = 1}^{n} (1 - {rz}_{j})^{m_{j}} = \exp \sum_{s = 1}^{\infty} H_{s} r^{s} = Π_{k = 1}^{\infty} \exp (H_{k} r^{k}),

(6)

we omit factors in the infinite product with

k > s

, multiply the exponential series for

\exp (H_{1} r), \dots, \exp (H_{s} r^{s})

, and pick out the coefficient of

r^{s}

Definition 1 Define

M = \sum_{j = 1}^{n} m_{j}, B = Π_{j = 1}^{n} b_{j}^{m_{j}}, d_{ij} = a_{i} b_{j} - a_{j} b_{i}, D_{i} = Π_{\overset{j \neq i}{Π}}^{n}

(7)

μ_{\pm s} (i) = \frac{- 1}{s} \sum_{\sum^{j \neq i}}^{n}

(8)

C_{0} (i) = 1, C_{\pm s} (i) = \sum \frac{μ_{\pm 1}^{α_{1}} (i) \dots μ_{\pm s}^{α_{s}} (i)}{α_{1}! \dots α_{s}!}, s = 1, 2, 3, \dots,

(9)

where upper (lower) signs go together and the last sum extends over all nonnegative integers

α_{1}, \dots, α_{s}

such that

α_{1} + 2 α_{2} + \dots + s α_{s} = s

. In particular we have

\begin{matrix} C_{\pm 1} (i) & = & μ_{\pm 1} (i), C_{\pm 2} (i) = \frac{1}{2} μ_{\pm 1}^{2} (i) + μ_{\pm 2} (i), \\ C_{\pm 3} (i) & = & \frac{1}{6} μ_{\pm 1}^{3} (i) + μ_{\pm 1} (i) μ_{\pm 2} (i) + μ_{\pm 3} (i) . & (10) \\ C_{\pm 4} (i) & = & \frac{1}{24} μ_{\pm 1}^{4} (i) + \frac{1}{2} μ_{\pm 1}^{2} (i) μ_{\pm 2} (i) + μ_{\pm 1} (i) μ_{\pm 3} (i) + \frac{1}{2} μ_{\pm 2}^{2} (i) + μ_{\pm 4} (i) . \end{matrix}

We note that

d_{ij} \neq 0

i \neq 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

\begin{matrix} \frac{P (t)}{Q (t)} = Π_{j = 1}^{n} (a_{j} + b_{j} t)^{m_{j}} & = & B b_{i}^{- M} \sum_{q = 0}^{M} C_{M - q} (i) (a_{i} + b_{i} t)^{q} \\ + & \sum_{i = 1}^{n} D_{i} b_{i}^{m_{i} - M} \sum_{q = 1}^{- m_{i}} C_{m_{i} + q} (i) (a_{i} + b_{i} t)^{- q} . & (11) \end{matrix}

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

d_{ij}

implies

a_{j} + b_{j} t = \frac{b_{j}}{b_{i}} (a_{i} + b_{i} t) (1 - \frac{d_{ij}}{b_{j} (a_{i} + b_{i} t)}),

(12)

whence

Π_{j = 1}^{n} (a_{j} + b_{j} t)^{m_{j}} = B b_{i}^{- M} (a_{i} + b_{i} t)^{M} Π_{j = 1}^{n} {(1 - \frac{d_{ij}}{b_{j} (a_{i} + b_{i} t)})}^{m_{j}} .

(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 / (a_{i} + b_{i} t)

and

z_{j} = d_{ij} / b_{j}

. It follows that

H_{s}

μ_{+ s} (i)

as defined by (2.6), the term in

H_{s}

with

j = i

being irrelevant because

d_{ii} = 0

. Thus we see that

T_{s}

C_{+ s} (i)

as defined by (2.7), and (2.11) becomes

\frac{P (t)}{Q (t)} = B b_{i}^{- M} \sum_{s = 0}^{\infty} C_{s} (i) (a_{i} + b_{i} t)^{M - s} .

(14)

Since this is valid for sufficiently large

t

, the polynomial part of

P / Q

consists of the terms with

0 \leq s \leq M

, which are absent unless

M \geq 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

a_{j} + b_{j} t = \frac{d_{ji}}{b_{i}} (1 - \frac{b_{j}}{d_{ij}} (a_{i} + b_{i} t)), j \neq i,

(15)

whence

Π_{j = 1}^{n} (a_{j} + b_{j} t)^{m_{j}} = D_{i} b_{i}^{m_{i} - M} (a_{i} + b_{i} t)^{m_{i}} Π_{\overset{j \neq i}{Π}}^{n}

(16)

‖ a_{i} + b_{i} t ‖

is sufficiently small, we can expand the last product by Lemma 2.1 with

r = a_{i} + b_{i} t

and

z_{j} = b_{j} / d_{ij}

. It follows that

H_{s}

μ_{- s} (i)

as defined by (2.6) and that

T_{s}

C_{- s} (i)

as defined by (2.7). Hence (2.11) becomes

\frac{P (t)}{Q (t)} = D_{i} b_{i}^{m_{i} - M} \sum_{s = 0}^{\infty} C_{- s} (i) (a_{i} + b_{i} t)^{m_{i} + s} .

(17)

Since this is valid if

‖ a_{i} + b_{i} t ‖

is sufficiently small, the part of

P / Q

representing a pole at

- a_{i} / b_{i}

consists of the terms with

0 \leq s \leq - m_{i} - 1

, which are absent unless

- m_{i} \geq 1

. These can be rewritten by putting

s = - m_{i} - q

D_{i} b_{i}^{m_{i} - M} \sum_{q = 1}^{- m_{i}} C_{m_{i} + q} (i) (a_{i} + b_{i} 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) = \int_{y}^{x} Π_{i = 1}^{h} (a_{i} + b_{i} t)^{- 1 / 2} Π_{j = 1}^{n} (a_{j} + b_{j} t)^{m_{j}} dt,

(19)

where

m = (m_{1}, \dots, m_{n})

is an

n

-tuple of integers. We define

n

-tuples

\begin{matrix} e_{0} = (0, 0, \dots, 0), \\ e_{1} = (1, 0, \dots, 0), & (20) \\ : \\ e_{n} = (0, \dots, 0, 1) . \end{matrix}

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 b_{i}^{- M} \sum_{q = 0}^{M} C_{M - q} (i) I ({qe}_{i}) + \sum_{i = 1}^{n} D_{i} b_{i}^{m_{i} - M} \sum_{q = 1}^{- m_{i}} C_{m_{i} + q} (i) I (- {qe}_{i}) .

(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 (\pm {qe}_{i})

in (2.19) can be expressed in terms of

I (e_{0})

and

I (\pm e_{i})

by using recurrence relations. The main relation involves also an algebraic term of the form

A (m) = v_{m} (x) - v_{m} (y)

, where

v_{m} (t)

has the form of an integrand in (2.17).

Theorem 2 Define

A (m) = v_{m} (x) - v_{m} (y),

(22)

v_{m} (t) = \frac{1}{s (t)} Π_{j = 1}^{n} (a_{j} + b_{j} t)^{m_{j}}, s (t) = Π_{i = 1}^{h} \sqrt{a_{i} + b_{i} t} .

(23)

For

1 \leq j \leq n

and

h = 3

4

let

E_{p} (j)

be the elementary symmetric function of degree

p

\frac{d_{1 j}}{b_{1}}, \dots, \frac{d_{hj}}{b_{h}},

(24)

and define

σ_{0} (j) = σ_{0} = Π_{i = 1}^{h} b_{i}, σ_{p} (j) = σ_{0} E_{p} (j),

(25)

whence

σ_{h} (j) = 0

1 \leq j \leq h

. Then, for any integer

q

and for

1 \leq j \leq n

\sum_{r = 0}^{h} (q + \frac{r}{2} + 1) σ_{h - r} (j) I ((q + r) e_{j}) = b_{j}^{h - 1} A ((q + 1) e_{j} + \sum_{i = 1}^{h} e_{i}) .

(26)

Choose

m = (q + 1) e_{j} + \sum_{i = 1}^{h} e_{i}

, whence

v_{m} (t) = s (t) (a_{j} + b_{j} t)^{q + 1},

(27)

v'_{m} (t) = s' (t) (a_{j} + b_{j} t)^{q + 1} + (q + 1) b_{j} s (t) (a_{j} + b_{j} t)^{q} .

(28)

From

a_{i} + b_{i} t = \frac{b_{i}}{b_{j}} (a_{j} + b_{j} t + \frac{d_{ij}}{b_{i}}),

it follows that

\begin{matrix} s^{2} (t) & = & \frac{σ_{0}}{b_{j}^{h}} Π_{i = 1}^{h} (a_{j} + b_{j} t + \frac{d_{ij}}{b_{i}}) \\ = & \frac{σ_{0}}{b_{j}^{h}} \sum_{r = 0}^{h} E_{h - r} (j) (a_{j} + b_{j} t)^{r} \\ = & \frac{1}{b_{j}^{h}} \sum_{r = 0}^{h} σ_{h - r} (j) (a_{j} + b_{j} t)^{r} & (29) \end{matrix}

and therefore

s' (t) = \frac{1}{2 b_{j}^{h - 1} s (t)} \sum_{r = 0}^{h} r σ_{h - r} (j) (a_{j} + b_{j} t)^{r - 1} .

(30)

Replacing

s

s^{2} / s

in the last term of (3.7) and substituting (3.8) and (3.9), we find

v'_{m} (t) = \frac{1}{b_{j}^{h - 1} s (t)} \sum_{r = 0}^{h} (\frac{r}{2} + q + 1) σ_{h - r} (j) (a_{j} + b_{j} t)^{q + r} .

(31)

Integration from

y

x

proves (3.5).

1 \leq j \leq h

one of the quantities listed in (3.3) is 0. Then

E_{h} (j) = σ_{h} (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 \leq j \leq h

h + 1

integrals if

h + 1 \leq j \leq n

. In the first case

I (\pm {qe}_{j})

can be reduced to the

h - 1

integrals

I (- e_{j}), I (e_{0}), \dots, I ((h - 3) e_{j})

, but in the second case

I ((h - 2) e_{j})

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

b_{i}^{q} I ({qe}_{j}) = \sum_{r = 0}^{q} (\begin{matrix} q \\ r \end{matrix}) b_{j}^{r} d_{ji}^{q - r} I ({re}_{i}) .

(32)

Raising both sides of

b_{i} (a_{j} + b_{j} t) = b_{j} (a_{i} + b_{i} t) + d_{ji}

(33)

to the power

q

and making a binomial expansion of the right side, we find

b_{i}^{q} (a_{j} + b_{j} t)^{q} = \sum_{r = 0}^{q} (\begin{matrix} q \\ r \end{matrix}) b_{j}^{r} (a_{i} + b_{i} t)^{r} d_{ji}^{q - r} .

(34)

Multiplying both sides by

Π_{k = 1}^{h} (a_{k} + b_{k} t)^{- 1 / 2}

and integrating from

y

x

, we obtain (3.11).

If we choose

i \leq 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 ({qe}_{i}), i \leq h

, with

q = - 1, 0

h = 3

and

q = - 1, 0, 1

h = 4

. The remaining terms of the second sum can be reduced by (3.5) to integrals of the form

I ({qe}_{j}), j > h

, with

q = - 1, 0, 1

h = 3

and

q = - 1, 0, 1, 2

h = 4

. Now

I ({qe}_{j}), q = 1, 2,

can be expressed by (3.11) in terms of

I ({re}_{i})

with

i \leq h

and

0 \leq r \leq 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 (- e_{j}), 0 \leq j \leq n,

in the cubic case (

h = 3

) plus

I (e_{i}), 1 \leq i \leq 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 (\sum_{j = 1}^{5} m_{j} e_{j})

displayed in the entry and defined in (1.2). Also displayed is the expression for

I (m)

in terms of the basic integrals

I (e_{0}), I (- e_{5})

, and

I (\pm e_{i}), i = 1, \dots, 4

. Equation (2.19) suffices for

I (e_{5})

and

I (e_{k} - e_{i})

; recurrence relations are not needed.

Table 1: Some integrals from Gradshteyn and Rhyzhik in terms of basic integrals

§ 3.147, 1-8:	$I (e_{0})$
§ 3.148, 1-8:	$I (e_{5}) = \frac{1}{b_{i}} [b_{5} I (e_{i}) - d_{i 5} I (e_{0})]$
§ 3.149, 1-8; § 3.151, 1-8]:	$I (- e_{5})$
§ 3.167, 1-32:	$I (e_{i})$
§ 3.168, 1-72:	$I (e_{k} - e_{i}) = \frac{1}{b_{i}} [d_{ki} I (- e_{i}) + b_{k} I (e_{0})]$

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

z_{1}, z_{2}, z_{3}

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

z_{1}, z_{2}, z_{3}

R_{F} (z_{1}, z_{2}, z_{3}) = \frac{1}{2} \int_{0}^{\infty} \frac{dt}{Π_{i = 1}^{3} \sqrt{t + z_{i}}},

(35)

R_{J} (z_{1}, z_{2}, z_{3}, z_{ν}) = \frac{3}{2} \int_{0}^{\infty} \frac{dt}{Π_{i = 1}^{3} \sqrt{t + z_{i}} (t + z_{ν})}, z_{ν} \neq 0 .

(36)

The functions

R_{F}

and

R_{J}

are respectively homogeneous of degree -1/2 and -3/2 in their variables. If

z_{ν}

is real and negative,

R_{J}

takes its Cauchy principal value. Useful degenerate cases are

R_{C} (z_{1}, z_{2}) = R_{F} (z_{1}, z_{2}, z_{2}) = \frac{1}{2} \int_{0}^{\infty} \frac{dt}{\sqrt{t + z_{1}} (t + z_{2})}, z_{2} \neq 0,

(37)

R_{D} (z_{1}, z_{2}, z_{3}) = R_{J} (z_{1}, z_{2}, z_{3}, z_{3}) = \frac{3}{2} \int_{0}^{\infty} \frac{dt}{Π_{i = 1}^{2} \sqrt{t + z_{i}} (t + z_{3})^{3 / 2}}, z_{3} \neq 0 .

(38)

z_{2}

is real and negative,

R_{C}

takes its Cauchy principal value,

R_{C} (z_{1}, - p) = \sqrt{\frac{z_{1}}{z_{1} + p}} R_{C} (z_{1} + p, p), p > 0 .

(39)

If the

z

's are distinct, then

R_{F}, R_{D}

, and

R_{J}

are elliptic integrals of the first, second, and third kinds, respectively, and

R_{C}

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

R_{F}

and

R_{J}

. Permuting the zeros of the quartic polynomial simply permutes the zeros of the cubic.

Theorem 4 Let $z_{1}, z_{2}, z_{3}, z_{4}$ 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

z_{ij} = z_{i} - z_{j}, V_{ij} = \sqrt{z_{i}} \sqrt{z_{j}} + \sqrt{z_{k}} \sqrt{z_{l}},

(40)

whence

V_{ij}^{2} - V_{ik}^{2} = z_{il} z_{jk} .

(41)

Assume that

V_{ij}

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

\int_{0}^{\infty} \frac{dt}{Π_{j = 1}^{4} \sqrt{t + z_{j}}} = 2 R_{F} (V_{12}^{2}, V_{13}^{2}, V_{23}^{2}) .

(42)

Let

z_{ν}

be complex and nonzero, and if

z_{ν} \neq z_{i}

define

z_{i ν} = z_{i} - z_{ν}, V_{i ν}^{2} = V_{ij}^{2} - \frac{z_{ik} z_{il} z_{j ν}}{z_{i ν}}, α_{i ν} = z_{ν} + \frac{\sqrt{z_{j}} \sqrt{z_{k}} \sqrt{z_{l}}}{\sqrt{z_{i}}},

(43)

where

α_{i ν}

is assumed to have positive real part or be 0. Then

\begin{matrix} \int_{0}^{\infty} \frac{t + z_{i}}{t + z_{ν}} \frac{dt}{Π_{j = 1}^{4} \sqrt{t + z_{j}}} & = & \frac{2 z_{ij} z_{ik} z_{il}}{3 z_{i ν}} R_{J} (V_{12}^{2}, V_{13}^{2}, V_{23}^{2}, V_{i ν}^{2}) \\ + 2 R_{C} (α_{i ν}^{2}, β_{i ν}^{2}), z_{ν} \neq z_{i}, & (44) \end{matrix}

β_{i ν}^{2} = \frac{z_{ν}}{z_{i}} V_{i ν}^{2}, α_{i ν}^{2} = β_{i ν}^{2} + \frac{z_{j ν} z_{k ν} z_{l ν}}{z_{i ν}} .

(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

R_{J}

, which contains a term in

R_{C}

. One can show that (4.10) reduces to (4.8) in the limit as

z_{ν} \to z_{i}

Because

V_{ij} = V_{ji} = V_{kl} = V_{lk}

(in particular,

V_{23} = V_{14}

), 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 $a_{i} + b_{i} x = X_{i}^{2}$ and $a_{i} + b_{i} y = Y_{i}^{2}$ lies in the complex plane cut along the negative real axis, and take all square roots in the right half-plane. Define

U_{ij} = \frac{X_{i} X_{j} Y_{k} Y_{l} + Y_{i} Y_{j} X_{k} X_{l}}{x - y},

(46)

whence

U_{ij}^{2} - U_{ik}^{2} = d_{il} d_{jk}, d_{ij} = a_{i} b_{j} - a_{j} b_{i} .

(47)

Assume that

d_{ij} \neq 0

i \neq j

and that

U_{ij}

has positive real part or is 0. (At most one of the

U

's is 0 because of (4.13).) Then

I (e_{0}) = 2 R_{F} (U_{12}^{2}, U_{13}^{2}, U_{23}^{2}),

(48)

where

I (e_{0})

is defined by (2.17) and (2.18).

Let

5 \leq ν \leq n

and assume the line segment with endpoints

a_{ν} + b_{ν} x = X_{ν}^{2}

and

a_{ν} + b_{ν} y = Y_{ν}^{2}

lies in the complex plane punctured at 0. Define

U_{i ν}^{2} = U_{ij}^{2} - \frac{d_{ik} d_{il} d_{j ν}}{d_{i ν}}, d_{i ν} = a_{i} b_{ν} - a_{ν} b_{i} \neq 0,

(49)

S_{i ν} = \frac{1}{x - y} (\frac{X_{j} X_{k} X_{l}}{X_{i}} Y_{ν}^{2} + \frac{Y_{j} Y_{k} Y_{l}}{Y_{i}} X_{ν}^{2}),

(50)

and assume

S_{i ν}

has positive real part or is 0. Then

I (e_{i} - e_{ν}) = \frac{2 d_{ij} d_{ik} d_{il}}{3 d_{i ν}} R_{J} (U_{12}^{2}, U_{13}^{2}, U_{23}^{2}, U_{i ν}^{2}) + 2 R_{C} (S_{i ν}^{2}, Q_{i ν}^{2}), d_{i ν} \neq 0,

(51)

Q_{i ν}^{2} = \frac{X_{ν}^{2} Y_{ν}^{2}}{X_{i}^{2} Y_{i}^{2}} U_{i ν}^{2}, S_{i ν}^{2} - Q_{i ν}^{2} = \frac{d_{j ν} d_{k ν} d_{l ν}}{d_{i ν}} .

(52)

Map the interval of integration in

I (e_{0}) = \int_{y}^{x} \frac{dt}{Π_{j = 1}^{4} \sqrt{a_{j} + b_{j} t}}

(53)

onto the positive real line by substituting

t = \frac{xs + y}{s + 1}, s = \frac{t - y}{x - t}, \frac{dt}{ds} = \frac{x - y}{(s + 1)^{2}},

(54)

a_{j} + b_{j} t = X_{j}^{2} \frac{s + z_{j}}{s + 1}, z_{j} = \frac{Y_{j}^{2}}{X_{j}^{2}}, z_{ij} = z_{i} - z_{j} = \frac{(x - y) d_{ij}}{X_{i}^{2} X_{j}^{2}} .

(55)

The result is

\begin{matrix} I (e_{0}) & = & N \int_{0}^{\infty} \frac{ds}{Π_{j = 1}^{4} \sqrt{s + z_{j}}}, N = \frac{x - y}{Π_{j = 1}^{4} X_{j}}, \\ = & 2 N R_{F} (V_{12}^{2}, V_{13}^{2}, V_{23}^{2}), \\ = & 2 R_{F} (U_{12}^{2}, U_{13}^{2}, U_{23}^{2}), U_{ij} = \frac{V_{ij}}{N}, & (56) \end{matrix}

where we have used (4.8) and the homogeneity of

R_{F}

. 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

\begin{matrix} I (e_{i} - e_{ν}) & = & \int_{y}^{x} \frac{a_{i} + b_{i} t}{a_{ν} + b_{ν} t} \frac{dt}{Π_{j = 1}^{4} \sqrt{a_{j} + b_{j} t}} \\ = & L \int_{0}^{\infty} \frac{s + z_{i}}{s + z_{ν}} \frac{dt}{Π_{j = 1}^{4} \sqrt{s + z_{j}}}, L = \frac{X_{i}^{2}}{X_{ν}^{2}} N, \\ = & L [\frac{2 z_{ij} z_{ik} z_{il}}{3 z_{i ν}} R_{J} (V_{12}^{2}, V_{13}^{2}, V_{23}^{2}, V_{i ν}^{2}) + 2 R_{C} (α_{i ν}^{2}, β_{i ν}^{2})] \\ = & L [\frac{2 X_{ν}^{2} N^{2} d_{ij} d_{ik} d_{il}}{3 X_{i}^{2} d_{i ν}} R_{J} (V_{12}^{2}, V_{13}^{2}, V_{23}^{2}, V_{i ν}^{2}) + 2 R_{C} (α_{i ν}^{2}, β_{i ν}^{2})] \\ = & \frac{2 d_{ij} d_{ik} d_{il}}{3 d_{i ν}} R_{J} (U_{12}^{2}, U_{13}^{2}, U_{23}^{2}, U_{i ν}^{2}) + 2 R_{C} (S_{i ν}^{2}, Q_{i ν}^{2}), & (57) \end{matrix}

S_{i ν} = \frac{α_{i ν}}{L}, Q_{i ν} = \frac{β_{i ν}}{L},

(58)

where we have used (4.10) and the homogeneity of

R_{J}

and

R_{C}

. 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 (- e_{ν})

from (4.17), we shall want the relation

d_{kl} I (m) = b_{l} I (m + e_{k}) - b_{k} I (m + e_{l}) .

(59)

This is proved by multiplying both sides of

d_{kl} = b_{l} (a_{k} + b_{k} t) - b_{k} (a_{l} + b_{l} t)

by the integrand of (2.17) and integrating from

y

x

Corollary 2 Under the same conditions as in Corollary 4.1 we have

\begin{matrix} I (e_{0}) & = & 2 R_{F} (U_{12}^{2}, U_{13}^{2}, U_{23}^{2}), & (60) \\ d_{i ν} I (- e_{ν}) & = & \frac{2 b_{ν}}{3 d_{i ν}} d_{ij} d_{ik} d_{i ℓ} R_{J} (U_{12}^{2}, U_{13}^{2}, U_{23}^{2}, U_{i ν}^{2}) + 2 b_{ν} R_{C} (S_{i ν}^{2}, Q_{i ν}^{2}) \\ - b_{i} I (e_{0}), & (61) \\ d_{ji} I (- e_{i}) & = & \frac{2}{3} b_{i} d_{jk} d_{j ℓ} R_{D} (U_{ik}^{2}, U_{jk}^{2}, U_{ij}^{2}) + \frac{2 b_{i} X_{j} Y_{j}}{X_{i} Y_{i} U_{ij}} - b_{j} I (e_{0}), & (62) \\ b_{i} I (e_{i}) & = & \frac{2}{3} d_{ij} d_{ik} d_{ℓ i} R_{J} (U_{12}^{2}, U_{13}^{2}, U_{23}^{2}, U_{i 0}^{2}) + 2 b_{i} R_{C} (S_{i 0}^{2}, Q_{i 0}^{2}), & (63) \\ b_{j} I (e_{i}) & = & b_{i} I (e_{j}) + d_{ij} I (e_{0}) . & (64) \end{matrix}

In (4.29) the quantities

U_{i 0}, S_{i 0}, Q_{i 0}

are obtained from (4.15), (4.16), and (4.18) by putting

a_{ν} = 1

and

b_{ν} = 0

, whence

X_{ν} = Y_{ν} = 1

and

d_{i ν} = - b_{i}

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) :

d_{i ν} I (- e_{ν}) = b_{ν} I (e_{i} - e_{ν}) - b_{i} I (e_{0}) .

(65)

To prove (4.28) we put

ν = j

in (4.27) and note that

S_{ij}^{2} - Q_{ij}^{2} = 0

implies

R_{C} (S_{ij}^{2}, Q_{ij}^{2}) = R_{C} (Q_{ij}^{2}, Q_{ij}^{2}) = \frac{1}{Q_{ij}} = \frac{X_{i} Y_{i}}{X_{j} Y_{j} U_{ij}} .

Then interchange

i

and

j

.
To prove (4.29) from (4.17) we replace

ν

by 0, where

a_{0} = 1

and

b_{0} = 0

, whence

e_{ν} = e_{0}

and

d_{i ν} = - b_{i}

.
Equation (4.30) is a special case of (4.25).

The free index

i

in (4.27) allows four choices for

U_{i ν}^{2}

, at most two of which can be real and negative. Thus

i

can be chosen to avoid a Cauchy principal value of

R_{J}

. 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

U_{ij} \neq 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

m_{4} = 0

a_{4} = 1

, and

b_{4} = 0

, whence

X_{4} = Y_{4} = 1

and

d_{i 4} = - b_{i}

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

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

R_{J}

and infinite values of

R_{D}

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) = \int_{y}^{x} \sqrt{\frac{(a_{1} + b_{1} t) (a_{2} + b_{2} t)}{a_{3} + b_{3} t}} \frac{dt}{(a_{5} + b_{5} t)^{2}}, m = (1, 1, 0, 0, - 2) .

(66)

We use the subscript 5 and take

m_{4} = 0

so that we can put

a_{4} = 1

and

b_{4} = 0

when the basic cubic integrals

I (- e_{j}), 0 \leq j \leq 5

, are expressed in terms of

R

-functions.

Applying (2.19) and Definition 2.1, we find

b_{5}^{2} I (m) = b_{1} b_{2} I (e_{0}) + (b_{1} d_{25} + b_{2} d_{15}) I (- e_{5}) + d_{15} d_{25} I (- 2 e_{5}),

(67)

where the last term needs to be reduced to basic integrals. Putting

h = 3

and

q = - 2

in (3.5), we have

2 σ_{3} (5) I (- 2 e_{5}) = - σ_{2} (5) I (- e_{5}) + σ_{0} I (e_{5}) - 2 b_{5}^{2} A (- e_{5} + \sum_{i = 1}^{3} e_{i}) .

(68)

To deal with

I (e_{5})

we use (3.11) to get

b_{i} I (e_{5}) = b_{5} I (e_{i}) - d_{i 5} I (e_{0}), 1 \leq i \leq 3,

(69)

and

I (e_{i})

is reduced to basic integrals by putting

h = 3

and

q = - 2

in (3.5) . Because

σ_{3} (i) = 0

we find

σ_{0} I (e_{i}) = σ_{2} (i) I (- e_{i}) + 2 b_{i}^{2} A (- e_{i} + \sum_{i = 1}^{3} e_{i}), 1 \leq i \leq 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

2 b_{5}^{2} I (m) = g_{0} I (e_{0}) + g_{i} I (- e_{i}) + g_{5} I (- e_{5}) + g_{a},

(71)

where

\begin{matrix} g_{0} & = & 2 b_{1} b_{2} - \frac{σ_{0} d_{15} d_{25} d_{i 5}}{b_{i} σ_{3} (5)} = 2 b_{1} b_{2} - \frac{b_{1} b_{2} b_{3} d_{i 5}}{b_{i} d_{35}}, & (72) \\ g_{i} & = & \frac{b_{5} d_{15} d_{25} σ_{2} (i)}{b_{i} σ_{3} (5)} = \frac{b_{5} d_{ji} d_{ki}}{d_{35}}, {i, j, k} = {1, 2, 3}, & (73) \\ g_{5} & = & 2 b_{1} d_{25} + 2 b_{2} d_{15} - \frac{d_{15} d_{25} σ_{2} (5)}{σ_{3} (5)} = b_{1} d_{25} + b_{2} d_{15} - \frac{b_{3} d_{15} d_{25}}{d_{35}}, & (74) \\ g_{a} & = & \frac{2 b_{5} d_{15} d_{25}}{σ_{3} (5)} [b_{i} A (- e_{i} + \sum_{i = 1}^{3} e_{i}) - b_{5} A (- e_{5} + \sum_{i = 3}^{3} e_{i})] \\ = & \frac{2 b_{5}}{d_{35}} [b_{i} A (- e_{i} + \sum_{i = 1}^{3} e_{i}) - b_{5} A (- e_{5} + \sum_{i = 3}^{3} e_{i})], & (75) \end{matrix}

where we have substituted the values of the

σ

'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

b_{i} d_{jk} + b_{j} d_{ki} + b_{k} d_{ij} = 0 .

(76)

The last two terms in square brackets in (6.10) can be combined by a companion of (4.25),

d_{kl} A (m) = b_{l} A (m + e_{k}) - b_{k} A (m + e_{l}),

(77)

whence

g_{a} = \frac{2 b_{5} d_{5 i}}{d_{35}} A (- e_{i} - e_{5} + \sum_{i = 1}^{3} e_{i}) = - 2 b_{5} A (e_{1} + e_{2} - e_{5}) .

(78)

Thus (6.6) becomes

\begin{matrix} 2 b_{5}^{2} I (m) & = & b_{1} b_{2} I (e_{0}) + \frac{b_{5} d_{13} d_{23}}{d_{35}} I (- e_{3}) \\ + & (b_{1} d_{25} - \frac{b_{5} d_{15} d_{23}}{d_{35}}) I (- e_{5}) - 2 b_{5} A (e_{1} + e_{2} - e_{5}) . & (79) \end{matrix}

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

U_{ij} \neq 0

, and (4.27) with

ν = 5

and

i

chosen as 1, 2, 3, or 4 so that

U_{i 5}^{2}

is not real and negative. Because we are putting

a_{4} = 1

and

b_{4} = 0

, the choice to try first in each case is 4 in order to give the extra term in

I (e_{0})

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

X_{4} = Y_{4} = 1

d_{j 4} = - b_{j}

, and

d_{45} = b_{5}

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) = \int_{0.5}^{2} \sqrt{\frac{(0.3 + 0.3 t) (0.5 + 0.1 t)}{0.7 - 0.1 t}} \frac{dt}{(0.9 - 0.3 t)^{2}} = 6.2430 9544,

(80)

obtained from (6.6) with

\begin{matrix} g_{0} = 0.03 & I (e_{0}) & = & 3.0973 71530 \\ g_{3} = 0.072 & I (- e_{3}) & = & 5.3785 84194 \\ g_{5} = - 0.18 & I (- e_{5}) & = & 6.1542 46948 \end{matrix}

g_{a} = - 2 b_{5} A (e_{1} + e_{2} - e_{5}) = - 2 b_{5} (\frac{X_{1} X_{2}}{X_{3} X_{5}^{2}} - \frac{Y_{1} Y_{2}}{Y_{3} Y_{5}^{2}}) = 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].

Footnotes:

^{1}

This 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.

^{2}

E-mail: carlson@decst2.ams.ameslab.gov

^{3}

Replacement of

r

r + 1

and

q

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 T_EX by T_TM, version 3.01.
On 5 Jul 2001, 00:33.

Toward Symbolic Integration of Elliptic Integrals

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

Abstract

1 Introduction

2 Partial-fraction decomposition

3 Recurrence relations and basic integrals

4 Reduction of basic integrals to R-functions

5 Summary of procedure

6 An example

Footnotes:

B. C. CARLSON $^{1}$ $^{2}$
Ames Laboratory and Department of Mathematics,
Iowa State University, Ames, Iowa 50011-3020

4 Reduction of basic integrals to $R$ -functions