SURFACE AREA AND CAPACITY OF ELLIPSOIDS IN n DIMENSIONS

GARRY J. TEE

Department of Mathematics, University of Auckland Auckland, New Zealand

2004–4–5

Abstract

The surface area of a general n-dimensional ellipsoid is represented as an Abelian integral, which can readily be evaluated numerically. If there are only 2 values for the semi-axes then the area is expressed as an elliptic integral, which reduces in most cases to elementary functions.

The capacity of a general n-dimensional ellipsoid is represented as a hyperelliptic integral, which can readily be evaluated numerically. If no more than 2 lengths of semi-axes occur with odd multiplicity, then the capacity is expressed in terms of elementary functions. If only 3 or 4 lengths of semi-axes occur with odd multiplicity, then the capacity is expressed as an elliptic integral.

Keywords: ellipsoid, n dimensions, surface area, capacity, Legendre, elliptic integral, hyperelliptic integral, Abelian integral.

AMS subject classification: primary 41A63 41A55, secondary 41-03.

1 Spheroids and ellipses

Consider an ellipsoid centred at the coordinate origin, with rectangular Cartesian coordinate axes along the semi-axes a, b, c,

\frac{x^{2}}{a^{2}} + \frac{y^{2}}{b^{2}} + \frac{z^{2}}{c^{2}} = 1.  (1)

1.1 Surface area of spheroid

In 1714, Roger Cotes found the surface area for ellipsoids of revolution [Cotes], called spheroids.

For the case in which two axes are equal b = c, the surface is generated by rotation around the x–axis of the half–ellipse \frac{x^{2}}{a^{2}} + \frac{y^{2}}{b^{2}} = 1 with y ≥ 0. On that half-ellipse, dy/dx = −b^{2}x/(a^{2}y), and hence the surface area of the spheroid is

A = 2 ∫_{0}^{a} 2πy\sqrt{1 + \frac{b^{4}x^{2}}{a^{4}y^{2}}} dx = 4π ∫_{0}^{a} \sqrt{y^{2} + \frac{b^{4}}{a^{4}} x^{2}} dx

= 4πb ∫_{0}^{a} \sqrt{1 − \frac{x^{2}}{a^{2}} + \frac{b^{2}}{a^{2}} \frac{x^{2}}{a^{2}}} dx

= 4πab ∫_{0}^{1} \sqrt{1 − (1 − \frac{b^{2}}{a^{2}} ) u^{2}} du = 4πab ∫_{0}^{1} \sqrt{1 − δu^{2}} du, (2)

where u = x/a and δ = 1 − b^{2}/a^{2}. Therefore, the surface areas for prolate spheroids (a > b), spheres (a = b) and oblate spheroids (a < b) are:

A = \{ \begin{matrix}2πb (a × \frac{arcsin √δ}{√δ} + b) & (prolate), \\ 2πb(a + b) = 4πa^{2} & (sphere), \\ 2πb (a × \begin{matrix}arcsinh & \sqrt{−δ} \\ √ & \overline{−δ} \end{matrix} + b) & (oblate) . \end{matrix}  (3)

Neither the hyperbolic functions nor their inverses had then been invented, and Cotes gave a logarithmic formula for the oblate spheroid [Cotes, pp. 169-171]. In modern notation [Cotes, p.50],

A = π [2a^{2} + b^{2} \frac{1}{\sqrt{−δ}} log (\frac{1 + \sqrt{−δ}}{1 − \sqrt{−δ}} )] .  (4)

For |δ| ≪ 1, either use the power series for (arcsin x)/x to get

A = 2πb (a [1 + \frac{1}{6}δ + \frac{3}{40}δ^{2} + \frac{5}{112}δ^{3} + ⋯] + b) ,  (5)

or else expand the integrand in (2) as a power series in u^{2} and integrate that term by term:

A = 4πab ∫_{0}^{1} (1 − δu^{2})^{1/2}du

= 4πab ∫_{0}^{1} (1 − \frac{1}{2}δu^{2} + \frac{\frac{1}{2} \frac{−1}{2}}{2!} δ^{2}u^{4} − \frac{\frac{1}{2} \frac{−1}{2} \frac{−3}{2}}{3!} δ^{3}u^{6}

+ \frac{\frac{1}{2} \frac{−1}{2} \frac{−3}{2} \frac{−5}{2}}{4!} δ^{4}u^{8} − \frac{\frac{1}{2} \frac{−1}{2} \frac{−3}{2} \frac{−5}{2} \frac{−7}{2}}{5!} δ^{5}u^{10} + ⋯) du

= 4πab (1 − \frac{1}{2} \frac{δ}{3} − \frac{1}{8} \frac{δ^{2}}{5} − \frac{1}{16} \frac{δ^{3}}{7} − \frac{5}{128} \frac{δ^{4}}{9} − \frac{7}{256} \frac{δ^{5}}{11} − ⋯) .  (6)

1.2 Circumference of ellipse

In 1742, Colin MacLaurin constructed a definite integral for the circumference of an ellipse [MacLaurin]. Consider an ellipse with semi-axes a and b, with Cartesian coordinates along the axes:

\frac{x^{2}}{a^{2}} + \frac{y^{2}}{b^{2}} = 1 ,  (7)

On that ellipse, 2x dx/a^{2} + 2y d y/b^{2} = 0, and hence dy/dx = −b^{2}x/(a^{2}y), and the circumference is 4 times the ellipse quadrant with x ≥ 0 and y ≥ 0. That quadrant has arclength

I = ∫_{0}^{a} \sqrt{1 + \frac{b^{4}x^{2}}{a^{4}y^{2}}} dx = ∫_{0}^{a} \sqrt{1 + \frac{b^{2}(x/a)^{2}}{a^{2}(y/b)^{2}}} dx

= ∫_{0}^{a} \sqrt{1 + \frac{(b/a)^{2}(x/a)^{2}}{1 − (x/a)^{2}}} dx .  (8)

Substitute z = x/a, and the circumference becomes

4I = 4a ∫_{0}^{1} \sqrt{1 + \frac{(b/a)^{2}z^{2}}{1 − z^{2}}} dz = 4a ∫_{0}^{1} \sqrt{\frac{1 − mz^{2}}{1 − z^{2}}} dz ,  (9)

where

m = 1 − \frac{b^{2}}{a^{2}} .  (10)

With a ≥ b this gives 0 ≤ m < 1.

That integral could not be expressed finitely in terms of standard functions. Many approximations for the circumference L(a, b) of an ellipse have been

published, and some of those give very close upper or lower bounds for L(a, b) [Barnard, Pearce & Schovanec]. A close approximation was given by Thomas Muir in 1883:

L(a, b) ≈ M(a, b) =^{def} 2π (\frac{a^{3/2} + b^{3/2}}{2} )^{2/3} .  (11)

That is a very close lower bound for all values of m ∈ (0, 1). Indeed, [Barnard, Pearce & Schovanec, (2)]:

0.00006m^{4} < \frac{L(a, b) − M(a, b)}{a} < 0.00666m^{4} .  (12)

1.3 Legendre on elliptic integrals

Adrien-Marie Legendre (1742-1833) worked on elliptic integrals for over 40 years, and summarized his work in [Legendre 1825]. He investigated systematically the integrals of the form ∫ R(t, y) dt, where R is a general rational function and y^{2} = P (t), where P is a general polynomial of degree 3 or 4. Legendre called them “fonctions élliptique”, because the formula (9) is of that form — now they are called elliptic integrals. He shewed how to express any such integral in terms of elementary functions, supplemented by 3 standard types of elliptic integral.

Each of Legendre’s standard integrals has 2 (or 3) parameters, including x = sin ϕ. Notation for those integrals varies considerably between various authors. Milne-Thomson’s notation for Legendre’s elliptic integrals [Milne-Thomson, §17.2] uses the parameter m, where Legendre (and many other authors) had used k^{2}.

Each of the three kinds is given as two integrals. In each case, the second form is obtained from the first by the substitutions t = sin θ and x = sin ϕ.

The Incomplete Elliptic Integral of the First Kind is:

F (ϕ|m) =^{def} ∫_{0}^{ϕ} \frac{dθ}{\sqrt{1 − m sin^{2} θ}} = ∫_{0}^{x} \frac{dt}{\sqrt{(1 − t^{2})(1 − mt^{2})}} . (13)

The Incomplete Elliptic Integral of the Second Kind is:

E(ϕ|m) =^{def} ∫_{0}^{ϕ} \sqrt{1 − m sin^{2} θ} dθ = ∫_{0}^{x} \sqrt{\frac{1 − mt^{2}}{1 − t^{2}}} dt .  (14)

That can be rewritten as

∫_{0}^{x} \frac{1 − mt^{2}}{\sqrt{(1 − t^{2})(1 − mt^{2})}} dt ,  (15)

which is of the form ∫ R(t, y) dt, where y^{2} = (1 − t^{2})(1 − mt^{2}).

The Incomplete Elliptic Integral of the Third Kind is:

Π(n; ϕ|m) =^{def} ∫_{0}^{ϕ} \frac{dθ}{(1 − n sin^{2} θ)\sqrt{1 − m sin^{2} θ}}

= ∫_{0}^{x} \frac{dt}{(1 − nt^{2})\sqrt{(1 − t^{2})(1 − mt^{2})}} .  (16)

The special cases for which ϕ = \frac{1}{2}π (and x = 1) are found to be particularly important, and they are called the Complete Elliptic Integrals [Milne-Thomson, §17.3].

The Complete Elliptic Integral of the First Kind is:

K(m) =^{def} F (\frac{1}{2}π|m) =^{def} ∫_{0}^{π/2} \frac{dθ}{\sqrt{1 − m sin^{2} θ}}

= ∫_{0}^{1} \frac{dt}{\sqrt{(1 − t^{2})(1 − mt^{2})}} .  (17)

The Complete Elliptic Integral of the Second Kind is:

E(m) =^{def} E (\frac{1}{2}π|m) = ∫_{0}^{π/2} \sqrt{1 − m sin^{2} θ} dθ = ∫_{0}^{1} \sqrt{\frac{1 − mt^{2}}{1 − t^{2}}} dt . (18)

The complete elliptic integrals K(m) and E(m) can efficiently be computed to high precision, by constructing arithmetic-geometric means [Milne-Thomson. §17.6.3 & 17.6.4].

2 Surface Area of 3-dimensional ellipsoid

For a surface defined by z = z(x, y) in rectangular Cartesian coordinates xyz, the standard formula for surface area is:

Area = ∫ ∫ \sqrt{1 + (\frac{∂z}{∂x})^{2} + (\frac{∂z}{∂y} )^{2}} dx dy.  (19)

On the ellipsoid (1),

\frac{∂z}{∂x} = \frac{−c^{2}x}{a^{2}z} ,  \frac{∂z}{∂y} = \frac{−c^{2}y}{b^{2}z} .  (20)

Consider the octant for which x, y, z are all non–negative. Then the surface area for that octant is

S = ∫_{0}^{a} ∫_{0}^{b\sqrt{1−x^{2}/a^{2}}}\sqrt{1 + \frac{c^{4}x^{2}}{a^{4}z^{2}} + \frac{c^{4}y^{2}}{b^{4}z^{2}}} dy dx

= ∫_{0}^{a} ∫_{0}^{b\sqrt{1−x^{2}/a^{2}}}\sqrt{\frac{\frac{z^{2}}{c^{2}} + \frac{c^{2}}{a^{2}} \frac{x^{2}}{a^{2}} + \frac{c^{2}}{b^{2}} \frac{y^{2}}{b^{2}}}{z^{2}/c^{2}}} dy dx

= ∫_{0}^{a} ∫_{0}^{b\sqrt{1−x^{2}/a^{2}}}\sqrt{\frac{1 − \frac{x^{2}}{a^{2}} − \frac{y^{2}}{b^{2}} + \frac{c^{2}}{a^{2}} \frac{x^{2}}{a^{2}} + \frac{c^{2}}{b^{2}} \frac{y^{2}}{b^{2}}}{1 − \frac{x^{2}}{a^{2}} − \frac{y^{2}}{b^{2}}}} dy dx

= ∫_{0}^{a} ∫_{0}^{b\sqrt{1−x^{2}/a^{2}}}\sqrt{\frac{1 − (1 − \frac{c^{2}}{a^{2}} ) \frac{x^{2}}{a^{2}} − (1 − \frac{c^{2}}{b^{2}} ) \frac{y^{2}}{b^{2}}}{1 − \frac{x^{2}}{a^{2}} − \frac{y^{2}}{b^{2}}}} dy dx . (21)

Hence, if two semi-axes (a and b) are fixed and the other semi-axis c increases, then the surface area increases.

Denote

δ = 1 − \frac{c^{2}}{a^{2}} ,  ϵ = 1 − \frac{c^{2}}{b^{2}} ,  (22)

and then (21) becomes

S = ∫_{0}^{a} ∫_{0}^{b\sqrt{1−x^{2}/a^{2}}}\sqrt{\frac{1 − δ \frac{x^{2}}{a^{2}} − ϵ\frac{y^{2}}{b^{2}}}{1 − \frac{x^{2}}{a^{2}} − \frac{y^{2}}{b^{2}}}} dy dx .  (23)

For a general ellipsoid, the coordinate axes can be named so that a ≥ b ≥ c > 0, and then 1 > δ ≥ ϵ ≥ 0.

2.1 Legendre’s series expansion for ellipsoid area

In 1788, Legendre converted this double integral to a convergent series [Legendre 1788] [Legendre 1825, pp. 350–351].

Replace the variables of integration (x, y) by (ϕ, θ), where cos θ = z/c, so that \frac{x^{2}}{a^{2}} + \frac{y^{2}}{b^{2}} = sin^{2} θ, or

\frac{x^{2}}{(a sin θ)^{2}} + \frac{y^{2}}{(b sin θ)^{2}} = 1 .  (24)

Then let cos ϕ = y/(b sin θ) so that sin ϕ = x/(a sin θ), or

x = a sin θ sin ϕ, y = b sin θ cos ϕ.  (25)

Differentiating x in (25) with respect to ϕ (with constant θ), we get that dx = a cos ϕ sin θ dθ; and differentiating the equation \frac{x^{2}}{a^{2}} + \frac{y^{2}}{b^{2}} = sin^{2} θ with respect to θ (with constant x), we get that 2y dy = 2b^{2} sin θ cos θ dθ. Thus the element of area in (23) becomes

dx dy = ab sin θ cos θ dθ dϕ,  (26)

and the area S of the ellipsoid octant becomes [Legendre 1825, p.350]

∫_{θ=0}^{π/2} ∫_{ϕ=0}^{π/2} \sqrt{\frac{1 − δ sin^{2} θ sin^{2} ϕ − ϵ (sin^{2} θ − \frac{x^{2}}{a^{2}} )}{1 − sin^{2} θ}} ab sin θ cos θ dϕ dθ

= ab ∫_{ϕ=0}^{π/2} ∫_{θ=0}^{π/2} sin θ \sqrt{1 − (δ sin^{2} ϕ + ϵ cos^{2} ϕ) sin^{2} θ} dθ dϕ . (27)

Thus,

S = ab ∫_{ϕ=0}^{π/2} ∫_{θ=0}^{π/2} sin θ \sqrt{1 − p sin^{2} θ} dθ dϕ ,  (28)

where p is a function of ϕ:

p = δ sin^{2} ϕ + ϵ cos^{2} ϕ = ϵ + (δ − ϵ) sin^{2} ϕ .  (29)

Hence, as ϕ increases from 0 to \frac{1}{2}π, p increases from ϵ ≥ 0 to δ < 1. Define

I(m) =^{def} ∫_{0}^{π/2} sin θ \sqrt{1 − m sin^{2} θ} dθ .  (30)

Clearly, I(m) is a decreasing function of m (for m ≤ 1). That integral can be expressed explicitly. For m ∈ (0, 1),

I(m) = \frac{1}{2} + \frac{1 − m}{4\sqrt{m}} log (\frac{1 + \sqrt{m}}{1 − \sqrt{m}}) .  (31)

Expand the integrand in (30) as a power series in p and integrate for θ from 0 to \frac{1}{2}π, to get a series expansion for I(p):

I(p) = ∫_{θ=0}^{π/2} sin θ (1 − p sin^{2} θ)^{1/2} dθ

= ∫_{0}^{π/2} sin θ (1 − \frac{1}{2}p sin^{2} θ + \frac{\frac{1}{2} \frac{−1}{2}}{2!} p^{2} sin^{4} θ − \frac{\frac{1}{2} \frac{−1}{2} \frac{−3}{2}}{3!} p^{3} sin^{6} θ + ⋯) dθ

= ∫_{0}^{π/2} sin θ dθ − \frac{1}{2}p ∫_{0}^{π/2} sin^{3} θ dθ − \frac{1·1}{2·4}p^{2} ∫_{0}^{π/2} sin^{5} θ dθ

− \frac{1·1·3}{2·4·6}p^{3} ∫_{0}^{π/2} sin^{7} θ dθ − \frac{1·1·3·5}{2·4·6·8}p^{4} ∫_{0}^{π/2} sin^{9} θ dθ − · · · . (32)

Define

s_{k} =^{def} ∫_{0}^{π/2} sin^{k} θ dθ = \{ \begin{matrix}\frac{π · 1 · 3 · 5 … (k − 1)}{2 · 2 · 4 · 6 … k} , & (even k ≥ 2) , \\ \frac{2 · 4 · 6 … (k − 1)}{3 · 5 · 7 … k} , & (odd k ≥ 3) , \end{matrix} (33)

with s_{0} = \frac{1}{2}π and s_{1} = 1. For all k > −1 [Dwight, §854.1],

s_{k} = \frac{\sqrt{π} Γ(\frac{1}{2}(1 + k))}{2Γ(1 + \frac{1}{2}k)} .  (34)

In particular,

s_{1} = 1, s_{3} = \frac{2}{3} , s_{5} = \frac{2·4}{3·5} , s_{5} = \frac{2·4·6}{3·5·7} , …,  (35)

and hence

I(p) = 1 − \frac{1}{1·3}p − \frac{1}{3·5}p^{2} − \frac{1}{5·7}p^{3} − \frac{1}{7·9}p^{4} − · · · .  (36)

Therefore, the surface area of the ellipsoid is

A = 8ab ∫_{0}^{π/2} I(p) dϕ
= 8ab ∫_{0}^{π/2} (1 − \frac{1}{1·3}p − \frac{1}{3·5}p^{2} − \frac{1}{5·7}p^{3} − \frac{1}{7·9}p^{4} − ⋯) dϕ
= 4πab(1 − \frac{1}{1·3}P_{1} − \frac{1}{3·5}P_{2} − \frac{1}{5·7}P_{3} − \frac{1}{7·9}P_{4} − ⋯) dϕ ,  (37)

where

P_{1} = \frac{2}{π} ∫_{0}^{π/2}(δ sin^{2} ϕ + ϵ cos^{2} ϕ) dϕ = \frac{1}{2}δ + \frac{1}{2}ϵ,
P_{2} = \frac{2}{π} ∫_{0}^{π/2}(δ sin^{2} ϕ + ϵ cos^{2} ϕ)^{2} dϕ = \frac{1·3}{2·4}δ^{2} + \frac{1·1}{2·2}δϵ + \frac{1·3}{2·4}ϵ^{2},
P_{3} = \frac{2}{π} ∫_{0}^{π/2}(δ sin^{2} ϕ + ϵ cos^{2} ϕ)^{3} dϕ
= \frac{1·3·5}{2·4·6}δ^{3} + \frac{1·3·1}{2·4·2}δ^{2}ϵ + \frac{1·1·3}{2·2·4}δϵ^{2} + \frac{1·3·5}{2·4·6}ϵ^{3}, et cetera. (38)

Legendre gave [Legendre 1825, p.51] a generating function for the P_{k}:

\frac{1}{\sqrt{(1 − δz)(1 − ϵz)}} = (1 − δz)^{−1/2}(1 − ϵz)^{−1/2}

= (1 + \frac{1}{2}δz + \frac{1·3}{2·4}δ^{2}z^{2} + \frac{1·3·5}{2·4·6}δ^{3}z^{3} + ⋯) (1 + \frac{1}{2}ϵz + \frac{1·3}{2·4}ϵ^{2}z^{2} + ⋯)
= 1 + P_{1}z + P_{2}z^{2} + P_{3}z^{3} + P_{4}z^{4} + · · · .  (39)

Infinite series had been used by mathematicians since the 13th century in India and later in Europe, but very little attention had been given to convergence. Consequently much nonsense had been published, resulting from the use of infinite series which did not converge. From 1820 onwards, Cauchy developed the theory of infinite series, and he stressed the importance of convergence [Grabiner, Chapter 4]. In 1825, Legendre carefully explained that his series (37) for the area does converge [Legendre 1825, p.351].

All terms after the first in Legendre’s series (37) are negative, and hence the partial sums of that series decrease monotonically towards the surface area.

I have searched many books on elliptic integrals and elliptic functions, and I have not found any later reference to Legendre’s series (37) for the surface area of a general ellipsoid.

Derrick H. Lehmer stated (in 1950) a different infinite series for the surface area, in terms of the eccentricities

α = \sqrt{1 − \frac{b^{2}}{a^{2}}} . β = \sqrt{1 − \frac{c^{2}}{a^{2}}} .  (40)

The surface area is

S(a, b, c) = 4πab [1 − \frac{1}{6}(α^{2} + β^{2}) − \frac{1}{120}(3α^{4} + 2α^{2}β^{2} + 3β^{4}) − ⋯]

= 4πab ∑_{ν=0}^{∞} \frac{(αβ)^{ν}}{1 − 4ν^{2}} P_{ν} (\frac{α^{2} + β^{2}}{2αβ} ) , (41)

where P_{ν}(x) is the Legendre polynomial of degree ν [Lehmer, (6)].

Philip Kuchel and Brian Bulliman constructed (in 1988) a more complicated series expansion for the area [Kuchel & Bulliman].

2.2 Bounds for ellipsoid area

As ϕ increases from 0 to \frac{1}{2}π, then δ sin^{2} ϕ+ϵ cos^{2} ϕ = (δ −ϵ) sin^{2} ϕ+ϵ increases from ϵ to δ. Hence, for all values of ϕ, the integrand in (27) lies between the upper and lower bounds

sin θ \sqrt{1 − ϵ sin^{2} θ} ≥

sin θ \sqrt{1 − (δ sin^{2} ϕ + ϵ cos^{2} ϕ) sin^{2} θ} ≥ sin θ \sqrt{1 − δ sin^{2} θ} . (42)

Accordingly, for all values of ϕ, the integral over θ in (27) lies between the upper and lower bounds

I(ϵ) ≥ ∫_{θ=0}^{π/2} sin θ \sqrt{1 − (δ sin^{2} ϕ + ϵ cos^{2} ϕ) sin^{2} θ} dθ ≥ I(δ) .  (43)

How close are those bounds? As λ increases from 0 to 1, I(λ) decreases from ∫_{θ=0}^{π/2} sin θ dθ = 1 to ∫_{θ=0}^{π/2} sin θ\sqrt{1 − sin^{2} θ} dθ = ∫_{θ=0}^{π/2} sin θ cos θ dθ = \frac{1}{2} . Therefore, the upper bound in (43) is not more than twice the upper bound.

Integrating these bounds over ϕ = 0 to ϕ = \frac{1}{2}π, we get upper and lower bounds for the double integral in (27):

\frac{π}{2} ≥ ∫_{ϕ=0}^{π/2} ∫_{θ=0}^{π/2} sin θ \sqrt{1 − (δ sin^{2} ϕ + ϵ cos^{2} ϕ) sin^{2} θ} dθ dϕ ≥ \frac{π}{4} . (44)

Thus, we get the following bounds for the surface area A = 8S of a general ellipsoid with semi-axes a ≥ b ≥ c:

4πab ≥ 8ab I(ϵ) ≥ A ≥ 8abI(δ) ≥ 2πab.  (45)

The extreme upper bound 4πab is attained for the sphere with a = b = c, and the extreme lower bound 2πab is attained with c = 0, when the ellipsoid collapses into a double-sided ellipse with semi-axes a and b.

2.3 Legendre’s explicit formula for ellipsoid area

In 1825, Legendre constructed [Legendre 1825, pp.352–359] an explicit expression for the area of a general ellipsoid, in terms of Incomplete Elliptic Integrals. In Milne-Thomson’s notation, with δ and ϵ as in (22),

A = 2πc^{2} + \frac{2πab}{\sqrt{δ}} [(1 − δ)F (\sqrt{δ}|ϵ/δ) + δ E(\sqrt{δ}|ϵ/δ)] .  (46)

Note that Legendre’s formula does not hold for a sphere, and for a near-sphere some rapidly-convergent series should be used for F/\sqrt{δ}, or else Legendre’s series (37) should be used for the area.

Legendre’s proof of his elliptic integral formula for area is long and complicated, and that formula has been very little known.

In 1953, Frank Bowman published [Bowman] an obscure derivation of the formula for the area (without mentioning Legendre). In 1958, Albert Eagle used his interesting version of elliptic functions to derive the formula for the area [Eagle, p.281 (12)]. He commented [Eagle, §10.31\frac{1}{2}]:

The formula (12) giving the solution of this really difficult problem was actually given by Legendre in his Traité des Fonctions elliptiques in 1832,^{1} obtained by means which only involved working with trig elliptic integrals, and not with elliptic functions as we now understand them. I think it is only appropriate that I should add a word here of admiration for the extraordinary mathematical ability of Legendre for solving such a difficult problem by means which it would be hopelessly beyond my abilities to accomplish. And yet I could, I suppose, if I tried, re-write all the steps of my solution in terms of the old trig integrals. But I am not making the attempt to do so!

I should also like to add a word of tribute to Legendre for the immense industry and labour he put into the calculation of his extensive tables of the First and Second Incomplete Elliptic Integrals. He laboured on elliptic integrals and their calculation for something like 40 years. And before his death in 1833, at the age of 81,^{2} he was a competent enough mathematician to realize that the new ideas that had very recently been put forth by Abel and Jacobi were right; and that the integrals he had spent his life considering were only the inverses of the functions that he should have been considering.

But the elliptic integral formula remained almost unknown. For example, in 1979 Stuart P. Keller asserted that “Except for the special cases of the sphere, the prolate spheroid and the oblate spheroid, no closed form expression exists for the surface area of the ellipsoid. This situation arises because of the fact that it is impossible to carry out the integration in the expression for the surface area in closed form for the most general case of three unequal axes”[Keller, p.310].

In 1989 Derek Lawden published a clear proof [Lawden, pp.100–102] of the formula (without mentioning Legendre), and in 1994 Leo Maas derived the elliptic integral formula [Maas], which he credits to Legendre. In 1999 an incorrect version of the formula was published [Wolfram, p.976], without proof or references.

1Rather, in 1825.

2Rather, at the age of 91.

The surface area of an ellipsoid was represented by Bille Carlson in terms of his very complicated function R_{G} [Carlson 1977, p.271]. In 2003, Edward Neumann gave a much simpler representation [Neumann, (1.4)] of Carlson’s function:

R_{G}(x, y, z) = \frac{1}{4} ∫_{0}^{∞} [\frac{\frac{x}{t + x} + \frac{y}{t + y} + \frac{z}{t + z}}{\sqrt{(t + x)(t + y)(t + z)}}] t dt .  (47)

But neither representation provides any clear way of computing the area.

3 Ellipsoid in n dimensions

Consider an n-dimensional ellipsoid (n ≥ 3) centred at the coordinate origin, with rectangular Cartesian coordinate axes along the semi-axes a_{1}, a_{2}, …, a_{n},

\frac{x_{1}^{2}}{a_{1}^{2}} + \frac{x_{2}^{2}}{a_{2}^{2}} + ⋯ + \frac{x_{n}^{2}}{a_{n}^{2}} = 1.  (48)

The area of the (n − 1)-dimensional surface is given by a generalization of (23):

A = 2^{n} ∫_{x_{1}=0}^{a_{1}} ∫_{x_{2}=0}^{a_{2}\sqrt{1− \frac{x^{2}_{1}}{a^{2}_{1}}}}⋯ ∫_{x_{n−1}=0}^{a_{n−1}\sqrt{1− \frac{x^{2}_{1}}{a^{2}_{1}} − \frac{x^{2}_{2}}{a^{2}_{2}} −···− \frac{x{2 \atop n−2}}{a{2 \atop n−2}}}}

\sqrt{\frac{1 − δ_{1} \frac{x_{1}^{2}}{a_{1}^{2}} − δ_{2} \frac{x_{2}^{2}}{a_{2}^{2}} − ⋯ − δ_{n−1} \frac{x_{n−1}^{2}}{a_{n−1}^{2}}}{1 − \frac{x_{1}^{2}}{a_{1}^{2}} − \frac{x_{2}^{2}}{a_{2}^{2}} − ⋯ − \frac{x_{n−1}^{2}}{a_{n−1}^{2}}}} dx_{n−1} ⋯ dx_{2} dx_{1}, (49)

where (cf. (22)):

δ_{i} =^{def} 1 − \frac{a_{n}^{2}}{a_{i}^{2}} ,  (i = 1, 2, …, n − 1).  (50)

As in 3 dimensions, if any semi-axis increases while the others are fixed, then the area increases.

Label the axes so that a_{1} ≥ a_{2} ≥ ⋯ ≥ a_{n}, and then 1 > δ_{1} ≥ δ_{2} ≥ ⋯ ≥ δ_{n−1} ≥ 0.

Scale each variable as x_{i} = a_{i}z_{i}, so that (52) reduces to integration over the unit n-sphere

z_{1}^{2} + z_{2}^{2} + ⋯ + z_{n}^{2} = 1.  (51)

Thus, the area of the ellipsoid is

A = 2^{n}a_{1}a_{2} ⋯ a_{n−1} ∫_{z_{1}=0}^{1} ∫_{z_{2}=0}^{\sqrt{1−z_{1}^{2}}}⋯ ∫_{z_{n−1}=0}^{\sqrt{1−z_{1}^{2}−z_{2}^{2}−···−z_{n−2}^{2}}}

\sqrt{\frac{1 − δ_{1}z_{1}^{2} − δ_{2}z_{2}^{2} − ⋯ − δ_{n−1}z_{n−1}^{2}}{1 − z_{1}^{2} − z_{2}^{2} − ⋯ − z_{n−1}^{2}}} dz_{n−1} ⋯ dz_{2} dz_{1} . (52)

3.1 Bounds for the area

The unit ball in n dimensions has volume [Smith & Vamanamurthy]

Ω_{n} = \frac{π^{n/2}}{Γ (1 + \frac{1}{2}n)} .  (53)

Thus,

Ω_{2} = 2π, Ω_{3} = \frac{4}{3}π, Ω_{4} = \frac{1}{2}π^{2}, Ω_{5} = \frac{8}{15}π^{2}, Ω_{6} = \frac{1}{6}π^{3} et cetera. (54)

For each j = 1, 2, ⋯, n, scaling the j-th semi-axis from 1 to a_{j} multiplies the volume by a_{j}, and hence the ellipsoid (48) has the volume

V = Ω_{n}a_{1}a_{2}a_{3} ⋯ a_{n} .  (55)

3.1.1 Carlson’s bounds

The surface area of an n-dimensional ellipsoid was represented by Bille Carlson in terms of his very complicated function R [Carlson 1966, (4.1)], which is defined by an integral representation [Carlson 1966, (2.1)]. From that, he proved the inequalities

(∏_{i=1}^{n} a_{i})^{−1/n} < \frac{1}{n} ∑_{i=1}^{n} \frac{1}{a_{i}} < \frac{A}{nV} < \sqrt{\frac{1}{n} ∑_{i=1}^{n} \frac{1}{a_{i}^{2}}} ,  (56)

where “a_{1}, ⋯, a_{n} are positive, finite and not all equal” [Carlson 1966, (4.3)]. Therefore,

Ω_{n}a_{1}a_{2} ⋯ a_{n} (\frac{1}{a_{1}} + \frac{1}{a_{2}} + ⋯ + \frac{1}{a_{n}} ) < A

< Ω_{n}a_{1}a_{2} ⋯ a_{n}\sqrt{\frac{n}{a_{1}^{2}} + \frac{n}{a_{2}^{2}} + ⋯ + \frac{n}{a_{n}^{2}}} . (57)

3.1.2 Refined bounds

In terms of b_{i} = 1/a_{i}, (56) becomes

n \sqrt{∏_{i=1}^{n} b_{i}} < \frac{1}{n} ∑_{i=1}^{n} b_{i} < \frac{A}{nV} < \sqrt{\frac{1}{n} ∑_{i=1}^{n} b_{i}^{2}} .  (58)

Thus, A/(nV) is bounded by the root-mean-square r = (b_{1}^{2} + b_{2}^{2} + ⋯ + b_{n}^{2})/n)^{1/2} of the b_{i} and their arithmetic mean m = (b_{1} + b_{2} + ⋯ + b_{n})/n, which is greater than the geometric mean of the b_{i} .

For positive numbers b_{i},

(b_{1} + b_{2} + ⋯ + b_{n})^{2} = b_{1}^{2} + b_{2}^{2} + ⋯ + b_{n}^{2} + ∑_{i≠j} b_{i}b_{j} > b_{1}^{2} + b_{2}^{2} + ⋯ + b_{n}^{2} , (59)

and so

∑_{i=1}^{n} b_{i} > \sqrt{∑_{i=1}^{n} b_{i}^{2}} ,  (60)

and hence

\frac{1}{n} ∑_{i=1}^{n} b_{i} > \frac{1}{\sqrt{n}}\sqrt{\frac{1}{n} ∑_{i=1}^{n} b_{i}^{2}} .  (61)

Therefore,

\sqrt{n} > \frac{\sqrt{\frac{1}{n} ∑_{i=1}^{n} b_{i}^{2}}}{\frac{1}{n} ∑_{i=1}^{n} b_{i}} .  (62)

Thus, Carlson’s inequalities give simple upper and lower bounds for the area of an n-dimensional ellipsoid, whose ratio is less than √n.

The variance of the inverse semi-axes is

v =^{def} \frac{1}{n} [(b_{1} − m)^{2} + (b_{2} − m)^{2} + ⋯ + (b_{n} − m)^{2}]

= \frac{1}{n} [b_{1}^{2} + b_{2}^{2} + ⋯ + b_{n}^{2} − 2m(b_{1} + b_{2} + ⋯ + b_{n}) + nm^{2}]

= \frac{1}{n} [b_{1}^{2} + b_{2}^{2} + ⋯ + b_{n}^{2} − nm^{2}] = r^{2} − m^{2} .  (63)

Therefore, the ratio of Carlson’s upper and lower bounds for A/(nV) is expressed in terms of the mean m of the inverse semi-axes and their standard deviation s = √v:

\frac{\sqrt{\frac{1}{n} ∑_{i=1}^{n} b_{i}^{2}}}{\frac{1}{n} ∑_{i=0}^{n} b_{i}} = \frac{r}{m} = \sqrt{\frac{r^{2}}{m^{2}}} = \sqrt{1 + \frac{v}{m^{2}}} = \sqrt{1 + (\frac{s}{m})^{2}} .  (64)

3.1.3 Examples

For a near-sphere with a_{1}/a_{n} not much larger than 1, (64) shews that the ratio of Carlson’s upper and lower bounds is much closer to 1 than to the √n in (62).

For example, consider ellipsoids in which the semi-axes are in arithmetic progression, from a_{1} = 1 to a_{n}. Table 1 gives, for each n and a_{n}, Ω_{n} and Carlson’s lower and upper bounds for the area and their ratio, followed by √n.

Table 1. Lower and upper bounds for surface area

n a_{n}  Ω_{n}  Low Bound Up Bound  Ratio √n
4  2 4.9348 × 10^{0}  6.251 × 10^{1}  6.465 × 10^{1}  1.0343 2
10  2 2.5502 × 10^{0}  8.124 × 10^{2}  8.327 × 10^{2}  1.0250 3.162
20  13 2.5807 × 10^{−2}  1.961 × 10^{14}  2.708 × 10^{14}  1.3813 4.472
256  3 1.1195 × 10^{−152} 1.547 × 10^{−78} 1.627 × 10^{−78} 1.0515 16
256  36 1.1195 × 10^{−152} 4.996 × 10^{147} 8.277 × 10^{147} 1.6568 16
256 100 1.1195 × 10^{−152} 3.428 × 10^{254} 7.804 × 10^{254} 2.2768 16

In each of these examples, including that with a_{256} = 100, the ratio of the bounds is substantially smaller than √n.

Bounds for area of spheroids If the semi-axes of an ellipsoid are a_{1} = a_{2} = ⋯ = a_{n−1} and a_{n} ≠ a_{1} (with n > 2), then we may call that a spheroid, oblate if a_{1} > a_{n} and prolate if a_{1} < a_{n}. That spheroid has volume (55)

V = Ω_{n}a_{1}^{n−1}a_{n} .  (65)

Carlson’s bounds in (57) become

m = \frac{(n − 1)b_{1} + b_{n}}{n} = b_{1} (1 + \frac{(b_{n}/b_{1}) − 1}{n}) ,

r = \sqrt{\frac{(n − 1)b_{1}^{2} + b_{n}^{2}}{n}} = b_{1}\sqrt{1 + \frac{(b_{n}/b_{1})^{2} − 1}{n}} .  (66)

Hence, for a sequence of spheroids with fixed semi-axes a_{1} and a_{n} as n → ∞,

m = \frac{1}{a_{1}} (1 + \frac{(a_{1}/a_{n}) − 1}{n}) and r = \frac{1}{a_{1}} (1 + O (n^{−1})) . (67)

Therefore, the area is asymptotically

A ∼ nV/a_{1} = nΩ_{n}a_{1}^{n−2}a_{n} .  (68)

However, Carlson’s representation of the area does not provide any clear way of computing that area.

4 Reduction to Abelian Integral

The infinite series for surface area of a 3-dimensional ellipsoid has been generalized to n-dimensional ellipsoids by using n-dimensional spherical coordinates — but that proves to be very complicated [Lehmer, (11)].

For 3-dimensional ellipsoids, Lawden used an ingenious method to reduce the double integral (23) to a single integral [Lawden, pp.100-102]. For n-dimensional ellipsoids, Lawden’s method is here generalized to reduce the (n−1)-dimensional integral (52) to an Abelian integral. That Abelian integral can readily be computed.

The multiple integral in (52) can be rewritten as

A =

2a_{1}a_{2} ⋯ a_{n−1} ∫_{z_{1}=−1}^{1} ∫_{z_{2}=−\sqrt{1−z_{1}^{2}}}^{\sqrt{1−z_{1}^{2}}} ⋯  ∫_{z_{n−1}=−\sqrt{1−z_{1}^{2}−z_{2}^{2}−···−z_{n−2}^{2}}}^{\sqrt{1−z_{1}^{2}−z_{2}^{2}−···−z_{n−2}^{2}}}

h^{−1} dz_{n−1} ⋯ dz_{2} dz_{1} ,  (69)

where

h =^{def} \sqrt{\frac{1 − z_{1}^{2} − z_{2}^{2} − ⋯ − z_{n−1}^{2}}{1 − δ_{1}z_{1}^{2} − δ_{2}z_{2}^{2} − ⋯ − δ_{n−1}z_{n−1}^{2}}} ,  (70)

with h^{−1} ≥ 1. In the (n − 1)-dimensional subspace with z_{n} = 0, the value of h ranges from 0 on that (n − 1)-sphere z_{1}^{2} + z_{2}^{2} + ⋯ + z_{n−1}^{2} = 1 to h = 1 at the centre.

At any point in the unit (n − 1)-ball which is the domain of integration in (52), any (n − 1)-dimensional element

dv = dz_{1} dz_{2} ⋯ dz_{n−1}  (71)

gets enlarged by the factor 1/h, to give an element of surface area of the ellipsoid:

a_{1}a_{2} ⋯ a_{n−1}h^{−1} dz_{1} dz_{2} ⋯ dz_{n−1} = h^{−1} dx_{1} dx_{2} ⋯ dx_{n−1} . (72)

Thus

A = 2a_{1}a_{2} ⋯ a_{n−1} ∫ \frac{dv}{h} ,  (73)

integrated over the unit (n − 1)-ball.

The points on the unit n-sphere with that value of h satisfy (70), and hence

h^{2}(1 − δ_{1}z_{1}^{2} − δ_{2}z_{2}^{2} − ⋯ − δ_{n−1}z_{n−1}^{2}) = 1 − z_{1}^{2} − z_{2}^{2} − ⋯ − z_{n−1}^{2} . (74)

Therefore, those points on the n-sphere project onto points in the (n−1)-sphere (with z_{n} = 0), forming the (n − 1)-dimensional ellipsoid

\frac{z_{1}^{2}}{(\frac{1 − h^{2}}{1 − δ_{1}h^{2}} )} + \frac{z_{2}^{2}}{(\frac{1 − h^{2}}{1 − δ_{2}h^{2}} )} + ⋯ + \frac{z_{n−1}^{2}}{(\frac{1 − h^{2}}{1 − δ_{n−1}h^{2}} )} = 1 .  (75)

(Ellipsoids of this type with different values of h are not confocal.)

This ellipsoid with parameter h has the semi-axes

f_{j} = \sqrt{\frac{1 − h^{2}}{1 − δ_{j}h^{2}}}  (j = 1, 2, ⋯, n − 1),  (76)

and hence its (n − 1)-dimensional measure is

w(h) = Ω_{n−1}g(h),  (77)

where

g(h) =^{def} f_{1}f_{2} ⋯ f_{n−1} = \frac{(1 − h^{2})^{(n−1)/2}}{\sqrt{(1 − δ_{1}h^{2})(1 − δ_{2}h^{2}) ⋯ (1 − δ_{n−1}h^{2})}} . (78)

Consider the (n − 1)-dimensional measure dw of the shell bounded by the (n − 1)-dimensional ellipsoids (75) with parameters h and h + dh. This is the projection of a strip of the n-dimensional sphere (51) (with parameter ranging from h to h + dh), and hence the element of area for that strip is dA = 2a_{1}a_{2} ⋯ a_{n−1}dv/h. It follows that the surface area between the point z_{n} = 1 (with h = 1) and the unit (n − 1)-sphere where z_{n} = 0 (with h = 0) is given by

S = ∫ \frac{dw}{h} = Ω_{n−1} ∫ \frac{dg(h)}{h} = Ω_{n−1} ∫_{1}^{0} \frac{1}{h} \frac{dg(h)}{dh} dh .  (79)

Denote z = −h^{2} and define

r(z) =^{def} (1 + δ_{1}z)(1 + δ_{2}z) ⋯ (1 + δ_{n−1}z),  (80)

and

p(z) =^{def} f_{1}^{2}f_{2}^{2} ⋯ f_{n−1}^{2} = \frac{(1 + z)^{n−1}}{r(z)} ,  (81)

so that

g(h) = \sqrt{p( − h^{2})},  (82)

and

\frac{dg}{dh} = \frac{dg}{dp} · \frac{dp}{dz} · \frac{dz}{dh} = \frac{−h}{\sqrt{p}} \frac{dp}{dz} .  (83)

Hence,

\frac{1}{h} \frac{dg(h)}{dh} = \frac{−1}{\sqrt{p(z)}} \frac{dp(z)}{dz} .  (84)

It follows from (81) and (80) that

\frac{dp}{dz} = \frac{(n − 1)(1 + z)^{n−2}}{r} − \frac{(1 + z)^{n−1}}{r^{2}} \frac{dr}{dz} ,  (85)

and

\frac{dr}{dz} = δ_{1}(1 + δ_{2}z)(1 + δ_{3}z) ⋯ (1 + δ_{n−1}z)

+(1 + δ_{1}z)δ_{2} ⋯ (1 + δ_{n−1}z) + (1 + δ_{1}z)(1 + δ_{2}z) ⋯ (1 + δ_{n−2}z)δ_{n−1}

= r (\frac{δ_{1}}{1 + δ_{1}z} + \frac{δ_{2}}{1 + δ_{2}z} + ⋯ + \frac{δ_{n−1}}{1 + δ_{n−1}z} ) ,  (86)

and hence

\frac{dp}{dz} = \frac{(1 + z)^{n−2}}{r(z)} (n − 1 − (1 + z) ∑_{i=1}^{n−1} \frac{δ_{i}}{1 + δ_{i}z} )

= \frac{(1 + z)^{n−2}}{r(z)} ∑_{i=1}^{n−1} (1 − \frac{δ_{i}(1 + z)}{1 + δ_{i}z}) = \frac{(1 + z)^{n−2}}{r(z)} ∑_{i=1}^{n−1} \frac{1 − δ_{i}}{1 + δ_{i}z} . (87)

It follows from (84) that

\frac{1}{h} \frac{dg(h)}{dh} = \frac{−1}{\sqrt{p(z)}} \frac{dp(z)}{dz}
= \frac{−(1 + z)^{(n−3)/2}}{\sqrt{r(z)}} (\frac{1 − δ_{1}}{1 + δ_{1}z} + ⋯ + \frac{1 − δ_{n−1}}{1 + δ_{n−1}z} )
= \frac{−(1 − h^{2})^{(n−3)/2}}{\sqrt{(1 − δ_{1}h^{2})(1 − δ_{2}h^{2}) ⋯ (1 − δ_{n−1}h^{2})}} ×
(\frac{1 − δ_{1}}{1 − δ_{1}h^{2}} + \frac{1 − δ_{2}}{1 − δ_{2}h^{2}} + ⋯ + \frac{1 − δ_{n−1}}{1 − δ_{n−1}h^{2}} ) .  (88)

This last expression is continuous on [0,1] for n ≥ 3, and hence it follows from (79) that the surface area of the ellipsoid is given explicitly as

A = 2a_{1}a_{2} ⋯ a_{n−1}Ω_{n−1} ∫_{1}^{0} \frac{1}{h} \frac{dg}{dh} dh

= \frac{2a_{1} ⋯ a_{n−1}π^{(n−1)/2}}{Γ (\frac{n+1}{2})} ∫_{0}^{1} \sqrt{\frac{(1 − h^{2})^{n−3}}{(1 − δ_{1}h^{2})(1 − δ_{2}h^{2}) ⋯ (1 − δ_{n−1}h^{2})}} ×

(\frac{1 − δ_{1}}{1 − δ_{1}h^{2}} + \frac{1 − δ_{2}}{1 − δ_{2}h^{2}} + ⋯ + \frac{1 − δ_{n−1}}{1 − δ_{n−1}h^{2}} ) dh ,  (89)

which is an Abelian integral.

For n = 2 this gives the circumference of the ellipse

\frac{x_{1}^{2}}{a_{1}^{2}} + \frac{x_{2}^{2}}{a_{2}^{2}} = 1 ,  (90)

with δ_{1} = 1 − (a_{2}/a_{1})^{2}, as

A = \frac{4a_{2}^{2}}{a_{1}} ∫_{0}^{1} \frac{dh}{(1 − δ_{1}h^{2})\sqrt{(1 − h^{2})(1 − δ_{1}h^{2})}} ,  (91)

which is an elliptic integral. Indeed, it was shewn in (9) that the circumference of that ellipse is 4a_{1}E(δ_{1}).

For n = 3, writing the semi-axes as a, b, c as in (1), and writing δ and ϵ as in (22), this becomes:

A = \frac{2πc^{2}}{ab} ∫_{0}^{1} \frac{a^{2} + b^{2} − (a^{2} + b^{2} − 2c^{2})h^{2}}{((1 − δh^{2})(1 − ϵh^{2}))^{3/2}} dh ,  (92)

and this reduces to the formula (46) in standard elliptic intehrals.

For odd n = 2k + 1 (with k > 0) the integrand in (89) is smooth, with the polynomial factor (1−h^{2})^{(n−3)/2} = (1−h^{2})^{k−1}, and hence this Abelian integral can readily be evaluated by Romberg integration. But for even n = 2k (with k > 1), the integrand has the factor (1 − h^{2})^{(n−3)/2} = [(1 + h)(1 − h)]^{k−3/2}. The factor (1 − h)^{k−3/2} makes the (k − 1)-th derivative of the integrand unbounded as h ↗ 1, so that direct numerical integration would be inefficient.

In order to prevent that near-singularity in the integrand for even n ≥ 4, make the substitution h = 1 − x^{2} so that 1 − h^{2} = x^{2}(2 − x^{2}), and (89) becomes the Abelian integral

A = \frac{4a_{1}a_{2} ⋯ a_{n−1}π^{(n−1)/2}}{Γ (\frac{n+1}{2})} ∫_{0}^{1} x^{n−2} ×

\sqrt{\frac{(2 − x^{2})^{n−3}}{(1 − δ_{1}(1 − x^{2})^{2})(1 − δ_{2}(1 − x^{2})^{2}) ⋯ (1 − δ_{n−1}(1 − x^{2})^{2})}} ×

(\frac{1 − δ_{1}}{1 − δ_{1}(1 − x^{2})^{2}} + ⋯ + \frac{1 − δ_{n−1}}{1 − δ_{n−1}(1 − x^{2})^{2}} ) dx .  (93)

This integrand is smooth for all n ≥ 2, and hence this Abelian integral can readily be computed by Romberg integration.

4.1 Examples

A program has been written (in THINK Pascal) to compute A by this formula. Table 2 gives the actual areas for those ellipsoids for which upper and lower bounds had been given in Table 1.

Table 2. Areas with lower and upper bounds

n a_{n} Lower Bound Computed Area  Upper Bound
4  2 6.251 × 10^{1}  6.39222065129 × 10^{1}  6.465 × 10^{1}
10  2 8.124 × 10^{2}  8.29183467191 × 10^{2}  8.327 × 10^{2}
20 13 1.961 × 10^{14}  2.58913095038 × 10^{14} 2.708 × 10^{14}
256 3 1.547 × 10^{−78} 1.62645053874 × 10^{−78} 1.627 × 10^{−78}
256 36 4.996 × 10^{147}  8.18623897626 × 10^{147} 8.277 × 10^{147}
256 100 3.428 × 10^{254}  7.56178774520 × 10^{254} 7.804 × 10^{254}

5 Equal semi-axes

If any of the semi-axes a_{i} share common values, then the corresponding pairs of factors under the square root sign in (89) (or in (93)) can be taken outside the square root sign. If the value of a_{j} occurs with odd multiplicity 2m+1 (amongst a_{1}, ⋯, a_{n−1}), then a single factor (1 − δ_{j}h^{2}) will remain in the denominator under the square root, with the factor (1 − δ_{j}h^{2})^{m} in the denominator of the integrand.

In particular, for odd n = 2k + 1, if only one value of a_{j} has odd multiplicity then select that as a_{n}. The square root term then reduces to 1, since all factors in the denominator, and also (1−h^{2})^{2k−2} in the numerator, can be taken outside the square root sign. In such a case the integrand in (89) is a rational function of h^{2}, and thatcan be integrated in terms of elementary functions. If only 3 values of a_{j} have odd multiplicities then select the smallest as a_{n} and the others as a_{1} and a_{2}, so that the integrand becomes a rational function of h^{2} divided by \sqrt{(1 − δ_{1}h^{2})(1 − δ_{2}h^{2})}. Thus, in that case the integral is an elliptic integral.

For even n, if all values of a_{j} have even multiplicity then select the smallest value as a_{n}, so that δ_{n−1} = 0. After all pairs of factors have been taken outside the square root the integrand becomes a rational function of h^{2} multiplied by \sqrt{1 − h^{2}}, and that can be integrated in terms of elementary functions by the substitution h = sin θ. If only 2 values of a_{j} have odd multiplicities then select the smaller as a_{n} and the other as a_{1}, so that the integrand becomes a rational function of h^{2} multiplied by \sqrt{(1 − h^{2})/(1 − δ_{1}h^{2})}. Thus, in that case the integral is an elliptic integral.

5.1 Biaxial ellipsoids

Consider biaxial ellipsoids of the type B_{p,q}, with p semi-axes equal to β and q semi-axes equal to γ, so that spheroids are of type B_{n−1,1}. Denote n = p + q, and then

β = a_{1} = a_{2} = ⋯ = a_{p} ≠ a_{p+1} = a_{p+2} = ⋯ = a_{n} = γ , (94)

and

δ_{1} = δ_{2} = ⋯ = δ_{p} = 1 − \frac{γ^{2}}{β^{2}} , δ_{p+1} = δ_{p+2} = ⋯ = δ_{n−1} = 0 . (95)

Hereafter, abbreviate δ_{1} to δ. Thus, 1 − δ = γ^{2}/β^{2}.

If either (or both) of β and γ has even multiplicity p = 2m then choose that as a_{1}, so that the ellipsoid is of type B_{2m,q}. With n = 2m + q. the formula (89) becomes

A = \frac{2β^{2m}γ^{q+1}π^{(n−1)/2}}{Γ (\frac{n+1}{2})} ∫_{0}^{1} \sqrt{\frac{(1 − h^{2})^{2m+q−3}}{(1 − h^{2})^{2m}}} (\frac{2m}{β^{2}(1 − δh^{2})} + \frac{q − 1}{γ^{2}}) dh

= \frac{2β^{2m}γ^{q−1}π^{(n−1)/2}}{Γ (\frac{n+1}{2})} ∫_{0}^{1} \frac{(1 − h^{2})^{(n−3)/2}}{(1 − δh^{2})^{m}} (\frac{2mγ^{2}}{β^{2}(1 − δh^{2})} + q − 1) dh. (96)

For odd q = 2k +1, the factor (1−h^{2})^{(n−3)/2} becomes the polynomial factor (1 − h^{2})^{m+k−1}. Thus, the area for ellipsoids of type B_{2m,2k+1} is

A = \frac{4β^{2m}γ^{2k}π^{m+k}}{(m + k)!} ∫_{0}^{1} \frac{(1 − h^{2})^{m+k−1}}{(1 − δh^{2})^{m}} (\frac{mγ^{2}}{β^{2}(1 − δh^{2})} + k) dh.  (97)

The integrand is a rational function of h^{2}, and hence this integral can be expressed in terms of elementary functions. The detailed working is given in Appendix A.

For even q = 2k,

(1 − h^{2})^{(n−3)/2} = \frac{(1 − h^{2})^{m+k−1}}{\sqrt{1 − h^{2}}} .  (98)

Thus, for biaxial ellipsoids of type B_{2m,2k} , the formula (89) reduces to

A = \frac{2β^{2m}γ^{2k}π^{m+k−1}√π}{Γ (m + k + \frac{1}{2})} ×

∫_{0}^{1} \frac{(1 − h^{2})^{m+k−1}}{(1 − δh^{2})^{m}\sqrt{1 − h^{2}}} (\frac{2mγ^{2}}{β^{2}(1 − δh^{2})} + 2k − 1) dh. (99)

The integrand is a rational function of h^{2} and \sqrt{1 − h^{2}}, which can be integrated in terms of elementary functions. The detailed working is given in Appendix B.

Thus, for biaxial ellipsoids with one (or both) semi-axes having even multiplicity, the surface area can be expressed in terms of elementary functions.

If neither of β and γ have even multiplicity then both have odd multiplicities 2m + 1 and 2k + 1, and n = 2k + 2m + 2 is even. Thus, for ellipsoids of type B_{2m+1,2k+1} the formula (89) becomes

A = \frac{2β^{2m+1}γ^{2k+2}π^{(2m+2k+1)/2}}{Γ (\frac{2m+2k+3}{2})} ×

∫_{0}^{1} \sqrt{\frac{(1 − h^{2})^{2m+2k−1}}{(1 − δh^{2})^{2m+1}}} (\frac{2m + 1}{β^{2}(1 − δh^{2})} + \frac{2k}{γ^{2}} ) dh

= \frac{2β^{2m+1}γ^{2k}π^{m+k}√π}{Γ (m + k + 1 + \frac{1}{2})} ×

∫_{0}^{1} \frac{(1 − h^{2})^{m+k}}{(1 − δh^{2})^{m}\sqrt{(1 − δh^{2})(1 − h^{2})}} (\frac{(2m + 1)γ^{2}}{β^{2}(1 − δh^{2})} + 2k) dh, (100)

which is an elliptic integral.

6 Capacity

James Clerk Maxwell gave (in 1873) the first extensive mathematical treatment of electrostatic capacity [Maxwell], in his epoch-making monograph^{3} A Treatise on Electricity and Magnetism. He defined Capacity of a Conductor in §50 [Maxwell v.1, p.48]:

If one conductor is insulated while all the surrounding conductors are kept at the zero potential by being put in communication with the earth, and if the conductor, when charged with a quantity E of electricity, has a potential V, the ratio of E to V is called the Capacity of the conductor. If the conductor is completely enclosed within a conducting vessel without touching it, then the charge on the inner conductor will be equal and opposite to the charge on the inner surface of the outer conductor, and will be equal to the capacity of the inner conductor multiplied by the difference of the potentials of the two conductors.

Maxwell then considered the General Theory of a System of Conductors: “Let A_{1}, A_{2}, … A_{n} be any number of conductors of any form. Let the charge or total quantity of electricity on each of them be E_{1}, E_{2}, … E_{n}, and let their potentials be V_{1}, V_{2}, … V_{n} respectively.” [Maxwell v.1, p.89]. He explained how

we should obtain n equations^{4} of the form [Maxwell v.1, p.90].

E_{1} = q_{11}V_{1} … + q_{1s}V_{s} … + q_{1n}V_{n},
− − − − − − − −
E_{r} = q_{r1}V_{1} … + q_{rs}V_{s} … + q_{rn}V_{n},
− − − − − − − −
E_{n} = q_{n1}V_{1} … + q_{ns}V_{s} … + q_{nn}V_{n}.  (101)

The coefficients in these equations ⋯ may be called Coefficients of Induction.

Of these q_{11} is numerically equal to the quantity of electricity on A_{1} when A_{1} is at potential unity and all the other bodies are at potential zero. This is called the Capacity of A_{1}. It depends on the form and position of all the conductors in the system.

Maxwell remarked that “The mathematical determination of the coefficients ⋯ of capacity from the known forms and positions of the conductors is in general difficult.” [Maxwell v.1, p.90]. He explained how to compute an upper bound for the capacity [Maxwell v.1, p.117], and he explained the method of J. W. Strutt (later Lord Rayleigh) for calculating upper and lower bounds for capacity “in the very important case in which the electrical action is entirely between two conducting surfaces S_{1} and S_{2}, of which S_{2} completely surrounds S_{1} and is kept at potential zero” [Maxwell v.1, p.118].

3My copy was formerly in the library of the Wellington Philosophical Society, and then in the library of the Royal Society of New Zealand.

4This is a noteworthy early example of a system of linear algebraic equations, written with double and single subscripts. GJT

6.1 Capacity of spheres

Maxwell considered “two concentric spherical surfaces of radii a and b, of which b be the greater, be maintained at potentials A and B respectively,” and he shewed that the potential at radius r between those spheres is of the form

V = C_{1} + C_{2}r^{−1} .  (102)

From this he found the capacity of the enclosed sphere:

C = \frac{ab}{b − a} .  (103)

[Maxwell v.1, p.153].

Fix the radius a of the inner sphere and let the radius b of the outer sphere increase without bound. As r ↗ ∞ the potential is

V = C_{1} + C_{2}r^{−1} → C_{1} .  (104)

and as b ↗ ∞ the capacity converges to the limit

C = \frac{a}{1 − (a/b)} ↘ a .  (105)

[Maxwell v.1, p.153].

Thus, we say that a sphere of radius a, in a 3-dimensional universe otherwise empty, has capacity equal to a.

And similar results hold (110) in space of more than 3 dimensions.

6.2 Capacity of coaxial cylinders

Maxwell considered Two Infinite Coaxal Cylindric Surfaces. “Let the radius of the outer surface of a conducting cylinder be a, and let the radius of the inner surface of a hollow cylinder, having the same axis with the first, be b. Let their potentials be A and B respectively”. In any plane orthogonal to the cylinder axis the potential satisfies the 2-dimensional Laplace equation, and Maxwell shewed that the potential at radius r between those cylinders is of the form

V = C_{1} + C_{2} log r .  (106)

From this he found that the capacity of a length l of the inner cylinder is

C = \frac{l}{2 log(b/a)} .  (107)

[Maxwell v.1, p.155].

Fix the radius a of the inner cylinder and let the radius b of the outer cylinder increase without bound. As the radius r ↗ ∞ the potential (106) increases (or decreases) unboundedly, and as b ↗ ∞ the capacity (107) of a length l of the inner cylinder converges to 0.

Thus, in a 2-dimensional potential field, the concept of capacity does not apply to an isolated finite body in an otherwise empty universe.

Accordingly, for the capacity of isolated figures we shall consider figures with dimension n ≥ 3.

7 Ellipsoid capacity in n dimensions

Consider the n-dimensional ellipsoid (48).

7.1 Bounds for the capacity

The electrostatic capacity C of an n-dimensional ellipsoid was represented by Bille Carlson in terms of his very complicated function R [Carlson 1966, (4.2)], which is defined by an integral representation [Carlson 1966, (2.1)]. From that, he obtained the inequalities

(∏_{i=1}^{n} a_{i})^{1−2/n} < \frac{C}{n − 2} < (\frac{1}{n} ∑_{i=1}^{n} a_{i})^{n−2} ,  (108)

where a_{1}, ⋯, a_{n} are positive, finite and not all equal [Carlson 1966, (4.9)]. Taking the (n − 2)-th roots in (108), this becomes

(∏_{i=1}^{n} a_{i})^{1/n} < (\frac{C}{n − 2})^{1/(n−2)} < \frac{1}{n} ∑_{i=1}^{n} a_{i} .  (109)

Thus, [C/(n − 2)]^{1/(n−2)} lies between the geometric mean of the semi-axes a_{i} and their arithmetic mean.

In particular, for an n-sphere of radius r,

C = (n − 2)r^{n−2} .  (110)

Carlson’s inequalities give simple upper and lower bounds for the capacity of an n-dimensional ellipsoid. But if a_{n} ↘ 0 (with fixed a_{1}, a_{2}, …, a_{n−1}), the geometric mean converges to 0 but the arithmetic mean converges to (a_{1} + ⋯ + a_{n−1})/n. Hence, for each n there is no upper bound to the ratio of Carlson’s upper and lower bounds for capacity.

Carlson’s representation of the capacity does not provide any clear way of computing that capacity.

8 Capacity as Hyperelliptic Integral

For a 3-dimensional ellipsoid, the capacity C can be expressed by a well-known elliptic integral [Pólya & Szegő, (6.4)]:

\frac{1}{C} = \frac{1}{2} ∫_{0}^{∞} \frac{du}{\sqrt{(a_{1}^{2} + u)(a_{2}^{2} + u)(a_{3}^{2} + u)}} .  (111)

For an n-dimensional ellipsoid with n ≥ 3, this generalizes to the integral

\frac{1}{C} = \frac{1}{2} ∫_{0}^{∞} \frac{du}{\sqrt{(a_{1}^{2} + u)(a_{2}^{2} + u) ⋯ (a_{n}^{2} + u)}} ,  (112)

which is an elliptic integral if n = 3 or 4, and is a hyperelliptic integral if n ≥ 5.

As a check on the coefficient of the integral, applying this formula to the n-sphere of radius r, we get

\frac{1}{C} = \frac{1}{2} ∫_{0}^{∞}(r^{2} + u)^{−n/2} du = \frac{1}{2} [\frac{(r^{2} + u)^{1−n/2}}{1 − n/2} ]_{0}^{∞} = \frac{r^{2−n}}{n − 2} , (113)

in agreement with (110).

8.1 Evaluation of integral

The hyperelliptic integrand in (112) is not suitable for numerical integration. Define

μ_{j} =^{def} 1 − \frac{a_{j}^{2}}{a_{1}^{2}} ,  (114)

so that each μ_{j} < 1. If the ellipsoid semi-axes are labelled so that

a_{1} ≥ a_{2} ≥ a_{3} ≥ ⋯ ≥ a_{n} > 0,  (115)

then

0 = μ_{1} ≤ μ_{2} ≤ μ_{3} ≤ ⋯ ≤ μ_{n} < 1 .  (116)

Then, with y = u/a_{1}^{2} we get

\frac{1}{C} = \frac{1}{2a_{1}^{n}} ∫_{0}^{∞} [(1 + \frac{u}{a_{1}^{2}} ) (\frac{a_{2}^{2}}{a_{1}^{2}} + \frac{u}{a_{1}^{2}} ) ⋯ (\frac{a_{n}^{2}}{a_{1}^{2}} + \frac{u}{a_{1}^{2}} )]^{−1/2} du

= \frac{1}{2a_{1}^{n−2}} ∫_{0}^{∞} [(1 + y) (\frac{a_{2}^{2}}{a_{1}^{2}} + y) ⋯ (\frac{a_{n}^{2}}{a_{1}^{2}} + y)]^{−1/2} dy
= \frac{1}{2a_{1}^{n−2}} ∫_{0}^{∞} [(1 + y) (1 + y − μ_{2}) ⋯ (1 + y − μ_{n})]^{−1/2} dy . (117)

Now substitute 1 + y = 1/x^{2}, so that dy = −2dx/x^{3}, and we get the hyperelliptic integral

\frac{1}{C} = \frac{1}{a_{1}^{n−2}} ∫_{0}^{1} \frac{x^{n−3}}{\sqrt{(1 − μ_{2}x^{2})(1 − μ_{3}x^{2}) ⋯ (1 − μ_{n}x^{2})}} dx .  (118)

For even n = 2q ≥ 4 we could substitute 1 + y = 1/z, so that dy = −dz/z^{2}, and we get the simpler hyperelliptic integral

\frac{1}{C} = \frac{1}{2a_{1}^{2q−2}} ∫_{0}^{1} \frac{z^{q−2}}{\sqrt{(1 − μ_{2}z)(1 − μ_{3}z) ⋯ (1 − μ_{2q}z)}} dz .  (119)

8.1.1 Capacity of 3-dimensional ellipsoid

In particular, for a 3-dimensional ellipsoid (not a sphere), substituting z = x\sqrt{μ_{3}},

\frac{1}{C} = \frac{1}{a_{1}} ∫_{0}^{1} \frac{dx}{\sqrt{(1 − μ_{2}x^{2})(1 − μ_{3}x^{2})}}

= \frac{1}{a_{1}} ∫_{0}^{√\overline{μ_{3}}} \frac{dz}{\sqrt{(1 − (\frac{μ_{2}}{μ_{3}} ) z^{2}) (1 − z^{2})}}

= \frac{1}{a_{1}\sqrt{μ_{3}}} ∫_{0}^{√\overline{μ_{3}}} \frac{dz}{\sqrt{(1 − (\frac{μ_{2}}{μ_{3}} ) z^{2}) (1 − z^{2})}}

= \frac{1}{\sqrt{a_{1}^{2} − a_{3}^{2}}} ∫_{0}^{√\overline{μ_{3}}} \frac{dz}{\sqrt{(1 − mz^{2})(1 − z^{2})}} = \frac{K(\sqrt{μ_{3}} | m)}{\sqrt{a_{1}^{2} − a_{3}^{2}}} , (120)

where m = μ_{2}/μ_{3} and K(x|m) is Legendre’s Incomplete Elliptic Integral of the First Kind (13). Therefore, the capacity is

C = \frac{\sqrt{a_{1}^{2} − a_{3}^{2}}}{K (arcsin \sqrt{1 − \frac{a_{3}^{2}}{a_{1}^{2}}} | \frac{a_{1}^{2} − a_{2}^{2}}{a_{1}^{2} − a_{3}^{2}} )} .  (121)

8.2 Examples

The integrand in (118) is smooth on the interval [0,1], and hence the integral can be evaluated readily by Romberg integration, to give the capacity C.

A program has been written (in THINK Pascal) to compute C by this formula. For the same ellipsoids whose surface areas have been calculated in Table 2, Table 3 gives the computed capacity, with Carlson’s lower and upper bounds.

Table 3. Capacities with lower and upper bounds

n a_{n} Lower Bound  Capacity  Upper Bound
4  2 4.216 × 10^{0}  4.40659279 × 10^{0}  4.5
10  2 1.697 × 10^{2}  1.81082801 × 10^{2}  2.050 × 10^{2}
20 13 8.854 × 10^{14}  2.14341967 × 10^{15} 2.931 × 10^{16}
256 3 6.832 × 10^{73} 8.30807716 × 10^{73} 7.352 × 10^{78}
256 36 1.991 × 10^{298}  1.15971902 × 10^{299} 1.847 × 10^{324}
256 100 4.269 × 10^{404}  5.99217774 × 10^{405} 1.099 × 10^{435}

9 Capacity for repeated semi-axes

If any of the semi-axes a_{1}, a_{2}, …, a_{n} share common values, then the corresponding pairs of factors under the square root sign in (112) can be taken outside the square root sign. If the value of a_{j} occurs with odd multiplicity 2m + 1, then the factor a_{j}^{2} + u will remain in the denominator under the square root, with the factor (a_{j}^{2} + u)^{m} in the denominator of the integrand.

Hence, if no value of a_{j} has odd multiplicity then the integrand in (112) is the reciprocal of a monic polynomial in u of degree n/2. If only one value of a_{j} has odd multiplicity then label it as a_{1}, and the integrand in (112) is

\frac{1}{P(u)\sqrt{a_{1}^{2} + u}} ,

where P is a monic polynomial of degree (n − 1)/2. If only 2 values of a_{j} have odd multiplicity then label them as a_{1} and a_{2}, and the integrand in (112) is

\frac{1}{Q(u)\sqrt{(a_{1}^{2} + u)(a_{2}^{2} + u)}} ,

where Q is a monic polynomial of degree (n − 2)/2. Therefore, if no more than 2 values of a_{j} have odd multiplicity then the hyperelliptic integral giving the capacity can be reduced to elementary functions of the semi-axes.

If only 3 values of a_{j} have odd multiplicity then label them as a_{1}, a_{2} and a_{3}, and the integrand in (112) is

\frac{1}{R(u)\sqrt{(a_{1}^{2} + u)(a_{2}^{2} + u)(a_{3}^{2} + u)}} ,

where R is a monic polynomial of degree (n − 3)/2. If only 4 values of a_{j} have odd multiplicity then label them as a_{1}, a_{2}, a_{3} and a_{4}, and the integrand in (112) is

\frac{1}{S(u)\sqrt{(a_{1}^{2} + u)(a_{2}^{2} + u)(a_{3}^{2} + u)(a_{4}^{2} + u)}} ,

where S is a monic polynomial of degree (n − 4)/2. Therefore, if only 3 or 4 values of a_{j} have odd multiplicity then the hyperelliptic integral (112) giving the capacity reduces to an elliptic integral of the semi-axes.

In particular, the capacity of a general ellipsoid in 3 or 4 dimensions is given by an elliptic integral, which reduces to elementary functions if any semi-axis is repeated.

9.1 Capacity of biaxial ellipsoids

In particular, consider biaxial ellipsoids of type B_{p,q}, which are defined as ellipsoids with p semi-axes equal to β and q semi-axes equal to γ, where (unlike (94)) we take β > γ. Denote n = p + q, and then

β = a_{1} = a_{2} = ⋯ = a_{p} > a_{p+1} = a_{p+2} = ⋯ = a_{n} = γ , (122)

and

μ_{1} = ⋯ = μ_{p} = 0, μ_{p+1} = ⋯ = μ_{n} = 1 − \frac{γ^{2}}{β^{2}} .  (123)

Hereafter, abbreviate μ_{n} to μ, where 0 ≤ μ < 1.

Then, (118) becomes

\frac{1}{C} = \frac{1}{β^{p+q−2}} ∫_{0}^{1} \frac{x^{p+q−3}}{\sqrt{(1 − μx^{2})^{q}}} dx .  (124)

Substitute x√μ = z and then z = sin θ, and this becomes

\frac{1}{C} = \frac{1}{(β√μ)^{p+q−2}} ∫_{0}^{√μ} \frac{z^{p+q−3}}{(1 − z^{2})^{q/2}} dz

= \frac{1}{(β^{2} − γ^{2})^{(p+q−3)/2}} ∫_{0}^{arcsin√μ} \frac{sin^{p+q−3} θ}{cos^{q−1} θ} dθ .  (125)

This can be integrated explicitly in terms of elementary functions — the full working is given in Appendix C.

9.2 Evaluation of integrals

A general elliptic integral can be reduced to expressions involving Legendre’s standard elliptic integrals. Legendre’s Incomplete Elliptic Integrals of the First and Second Kinds can be computed efficiently [Reinsch & Raab], and the Complete Integrals K(m) and E(m) can be computed very efficiently [Milne-Thomson §17.6.3 & §17.6.4].

But for a general elliptic integral (e.g. (100)), reduction to the standard elliptic integrals is an extremely complicated operation [Milne-Thomson §17.1]. It is usually simplest to evaluate elliptic integrals directly by Romberg integration, possibly after performing some substitution (e.g. x = sin θ) to make the integrand smooth.

And similarly for complicated integrals which could be expressed in terms of elementary functions, as in the Appendices.

Acknowledgement. I thank Reinhard Klette and Igor Rivin, for bringing to my attention the problem of estimating surface area of ellipsoids, and Matti Vuorinen for providing me with some useful modern references.

References

Roger W. Barnard, Kent Pearce & Lawrence Schovanec, Inequalities for the perimeter of an ellipse, Journal of Mathematical Analysis and Applications 260 (2001), 295-306. Mathematical Reviews 2002c:33012.

Frank Bowman, Introduction to Elliptic Functions, with Applications, Methuen, London, 1953. Mathematical Reviews 0058760 (15420d).

Bille C. Carlson, Some inequalities for hypergeometric functions, Proceedings of the American Mathematical Society 17 (1966), 32–39. Mathematical Reviews 0188497 (32 #5935).

Bille C. Carlson, Special Functions of Applied Mathematics, Academic Press, New York, 1977. Mathematical Reviews 0590943 (58 #28707).

Roger Cotes, Logometria, Philosophical Transactions of the Royal Society of London Vol. 29, No.338, 1714. English translation from Latin by Ronald Gowing in Roger Cotes — Natural Philosopher, Cambridge University Press, Cambridge, 1983, Appendix 1. Mathematical Reviews 87b:01033.

Herbert Bristol Dwight, Tables of Integrals and Other Mathematical Data, 3rd edition, The Macmillan Company, New York, 1957.

Albert Eagle, The Elliptic Functions As They Should Be. An account, with applications, of the functions in a new canonical form, Galloway & Porter, Ltd; Cambridge, 1958. Mathematical Reviews 0093599 (20 #123).

Judith V. Grabiner, The Origins of Cauchy’s Rigorous Calculus, M.I.T; Cambridge Mass; 1981. Mathematical Reviews 82d:01013.

Stuart P. Keller, On the surface area of the ellipsoid, Mathematics of Computation vol.33 (January 1979), 310–314. Mathematical Reviews 80k:65022.

Philip W. Kuchel & Brian T. Bulliman, A series expression for the surface area of an ellipsoid and its application to the computation of the surface area of avian erythrocytes, Journal of Theoretical Biology vol.134 (1988), 113-123.

Derek F. Lawden, Elliptic Functions and Applications, Springer-Verlag, Berlin, 1989. Mathematical Reviews 90h:33001.

Adrien-Marie Legendre, Sur les intégrales doubles, Mémoires de l’Académie des Sciences (1788). (not seen).

Adrien-Marie Legendre, Traite des Fonctions Élliptiques, tome 1, Huzard-Courchier, Paris 1825.

Derrick H. Lehmer, Approximations to the area of an n-dimensional ellipsoid, Canadian Journal of Mathematics 2 (1950), 267-282. Mathematical Reviews 0036013(12,45g).

Leo R. M. Maas, On the surface area of an ellipsoid and related integrals of elliptic integrals, Journal of Computational and Applied Mathematics, vol.51 (1994), 237–249. Mathematical Reviews 95i:33056.

Colin MacLaurin, Treatise of Fluxions, 2 volumes, Edinburgh, 1742. James Clerk Maxwell, A Treatise on Electricity and Magnetism (2

volumes), Clarendon Press, Oxford, 1873.

Louis Melville Milne-Thomson, Jacobian Elliptic Functions and Theta Functions (Chapter 16) & Elliptic Integrals (Chapter 17), in: Handbook of Mathematical Functions, with Formulas, Graphs and Mathematical Tables, ed. by Milton Abramowitz & Irene A. Stegun, National Bureau of Standards Applied Mathematics Series • 55, U.S. Government Printing Office, Washington, 1964 (and many later corrected editions). Mathematical Reviews 0167642 (29 #4914).

Edward Neumann, Bounds for symmetric elliptic integrals, Journal of Approximation Theory vol.122 (2003), 249-259.

George Pólya & Gábor Szegő, Inequalities for the capacity of a condenser, American Journal of Mathematics 67 (1945), 1-33. Mathematical Reviews 001187111 (6,227e).

Klaus-Dieter Reinsch & Winfried Raab, Elliptic integrals of the first and second kind — comparisons of Bulirsch’s and Carlson’s algorithms for numerical calculation, Special Functions, World Science Publishing, River Edge, 2000, 293-308 (not seen). Mathematical Reviews 2001k:33044.

David J. Smith & Mavina K. Vamanamurthy, How small is a unit ball?, Mathematics Magazine Vol.62 No.2 (April 1989), 103-107. Mathematical Reviews 90e:51029.

Stephen Wolfram, The Mathematica Book (4th edition), Cambridge University Press, Cambridge, 1999.

Appendices

A Surface area of biaxial ellipsoids B_{2m,2k+1}

Here, m ≥ 1 and k ≥ 0.

Denote

J_{m,k} = ∫_{0}^{1} \frac{(1 − h^{2})^{m+k−1}}{(1 − δh^{2})^{m+1}} dh ,  (126)

where δ = 1 − (γ/β)^{2}, and then the integral in (97) can be represented as

m(γ/β)^{2}J_{m,k} + kJ_{m−1,k+1} .  (127)

The integrand in (126) is

L = (\frac{1 − h^{2}}{1 − δh^{2}} )^{m+1} (1 − h^{2})^{k−2} .  (128)

This will be converted to partial fractions.

Denote x = 1/(1 − δh^{2}), so that

\frac{1 − h^{2}}{1 − δh^{2}} = \frac{1 + (δ − 1)x}{δ} = \frac{1 + u}{δ} ,  (129)

where u = (δ − 1)x, and

1 − h^{2} = (1 − 1/δ)(1 + 1/u) = \frac{(δ − 1)(1 + u)}{δu} .  (130)

Thus, the integrand is

L = (\frac{1 + u}{δ})^{m+1} (\frac{δ − 1}{δ})^{k−2} (\frac{1 + u}{u})^{k−2}
= \frac{(δ − 1)^{k−2}(1 + u)^{m+k−1}}{δ^{m+k−1}u^{k−2}}
= \frac{(δ − 1)^{k−2}}{δ^{m+k−1}} u^{2−k} ∑_{q=0}^{m+k−1} \binom{m + k − 1}{q}u^{q}

= \frac{(δ − 1)^{k−2}}{δ^{m+k−1}} ∑_{q=0}^{m+k−1} \binom{m + k − 1}{q}u^{q−k+2}

= \frac{1}{δ^{m+k−1}} ∑_{q=0}^{m+k−1} \binom{m + k − 1}{q}\frac{(δ − 1)^{q}}{(1 − δh^{2})^{q−k+2}}

= \frac{1}{δ^{m+k−1}} ∑_{q=0}^{m+k−1} \binom{m + k − 1}{q}(δ − 1)^{q}(1 − δh^{2})^{k−q−2} . (131)

This partial fraction expansion of L can then be integrated, to give

J_{m,k} = \frac{1}{δ^{m+k−1}} ∑_{q=0}^{m+k−1} \binom{m + k − 1}{q}(δ − 1)^{q} ∫_{0}^{1}(1 − δh^{2})^{k−q−2} dh . (132)

If k > 1 and q < k−1 then k−q−2 ≥ 0, and so the integrand is a polynomial in h^{2} of degree k − q − 2; but for q ≥ k − 1 the integrand is the reciprocal of a polynomial in h^{2} of degree q + k − 2. Thus, we need to consider

P_{g} =^{def} ∫_{0}^{1}(1 − δh^{2})^{g} dh ,  (133)

for integer g.

Denote X = 1 − δh^{2}, so that dX = −2δh dh.

Integrating by parts,

P_{g} = ∫_{0}^{1} X^{g} dh = hX^{g}|_{0}^{1} − ∫ h d(X^{g})
= (1 − δ)^{g} − g ∫ hX^{g−1}dX = (1 − δ)^{g} + 2g ∫_{0}^{1} δh^{2}X^{g−1} dh
= (1 − δ)^{g} + 2g ∫_{0}^{1} [X^{g−1} − (1 − δh^{2})X^{g−1}] dh
= (1 − δ)^{g} + 2g ∫_{0}^{1} [X^{g−1} − X^{g}] dh
= (1 − δ)^{g} + 2gP_{g−1} − 2gP_{g} ,  (134)

and hence

(2g + 1)P_{g} = 2gP_{g−1} + (1 − δ)^{g} .  (135)

Thus we get the reduction formulæ

P_{g} = \frac{2gP_{g−1} + (1 − δ)^{g}}{2g + 1}  (136)

for all g, and

P_{g−1} = \frac{(2g + 1)P_{g} − (1 − δ)^{g}}{2g}  (137)

for all nonzero g.

For integrating the polynomials with g ≥ 0, start with P_{0} = 1, and then apply the reduction formula (136) successively for g = 1, 2, 3, ⋯, to generate the integrals P_{g} for g ≥ 0.

If β > γ > 0 then 0 < δ < 1, but if γ > β > 0 then δ < 0. For integrating the rational functions with g < 0, start with

P_{−1} = ∫_{0}^{1} \frac{dh}{1 − δh^{2}} = \{ \begin{matrix}\frac{1}{2√δ} log (\frac{1 + √δ}{1 − √δ} ) & (β > γ, δ > 0), \\ \frac{arctan(\sqrt{−δ})}{\sqrt{−δ}} & (β < γ, δ < 0) . \end{matrix}  (138)

Then apply the reduction formula (137) successively for g = −1, −2, −3, ⋯, to generate the integrals P_{−2}, P_{−3}, P_{−4}, · · · .

Thus, the integral J_{m,k} has been constructed in terms of elementary functions.

Then the integral in (97) can be computed (127) from J_{m,k} and J_{m−1,k+1}, to give explicitly the surface area for a biaxial ellipsoid of type B_{2m,2k+1}, in terms of elementary functions of its semi-axes β and γ.

B Surface area of biaxial ellipsoids B_{2m,2k}

Here, m ≥ 1 and k ≥ 1.

Denote

H_{m,k} = ∫_{0}^{1} \frac{(1 − h^{2})^{m+k−1}}{(1 − δh^{2})^{m+1}} \frac{dh}{\sqrt{1 − h^{2}}} ,  (139)

where δ = 1 − (γ/β)^{2}, and then the integral in (99) can be represented as

2m(γ/β)^{2}H_{m,k} + (2k − 1)H_{m−1,k+1} .  (140)

This integrand for H_{m,k} is L/\sqrt{1 − h^{2}}, where L is the integrand (128) for J_{m,k}. Corresponding to the sum (131) of integrals for J_{m,k}, we get the sum of integrals

H_{m,k} =

\frac{1}{δ^{m+k−1}} ∑_{q=0}^{m+k−1} \binom{m + k − 1}{q}(δ − 1)^{q} ∫_{0}^{1} \frac{dh}{(1 − δh^{2})^{q−k+2}\sqrt{1 − h^{2}}} . (141)

Corresponding to P_{j} in (133), we consider

U_{j} =^{def} ∫_{0}^{1} \frac{dh}{(1 − δh^{2})^{j}\sqrt{1 − h^{2}}} ,  (142)

and the substitution h = sin ψ converts this to

U_{j} = ∫_{0}^{π/2} \frac{dψ}{(1 − δ sin^{2} ψ)^{j}} .  (143)

In (141), if k > 1 and q < k − 1 then q − k + 2 ≤ 0, and so the integrand is a polynomial in sin^{2} ψ. For q = 0 to k − 2, the integral in (141) is

U_{q−k+2} = ∫_{0}^{π/2}(1 − δ sin^{2} ψ)^{k−q−2} dψ

= ∫_{0}^{π/2} ∑_{r=0}^{k−q−2} \binom{k − q − 2}{r}(−δ)^{r} sin^{2r} ψ dψ

= ∑_{r=0}^{k−q−2} \binom{k − q − 2}{r}(−δ)^{r} ∫_{0}^{π/2} sin^{2r} ψ dψ

= ∑_{r=0}^{k−q−2} \binom{k − q − 2}{r}(−δ)^{r}s_{2r} ,  (144)

(cf. (33)). In particular, U_{0} = \frac{1}{2}π.

But for q ≥ k − 1 the integrand is the reciprocal of a polynomial in sin^{2} ψ, and so we need to consider U_{j} in (142) for positive integer j.

For j ≥ 1, from (143) we get

U_{j} = ∫_{0}^{π/2} \frac{dψ}{(1 − δ + δ cos^{2} ψ)^{j}} = ∫_{0}^{π/2} \frac{sec^{2j} ψ}{((1 − δ) sec^{2} ψ + δ)^{j}} dψ . (145)

Substitute f = tan ψ so that df = sec^{2} ψ dψ = (1 + f^{2}) dψ, and we get

U_{j} = ∫_{0}^{∞} \frac{(1 + f^{2})^{j}}{(1 + f^{2})[(1 − δ)(1 + f^{2}) + δ]^{j}} df

= ∫_{0}^{∞} \frac{(1 + f^{2})^{j−1}}{(1 + (1 − δ)f^{2})^{j}} df = ∫_{0}^{∞} \frac{(1 + f^{2})^{j−1}}{(1 + ηf^{2})^{j}} df .  (146)

where

η =^{def} 1 − δ = (\frac{γ}{β} )^{2} > 0 .  (147)

For j = 1, substitute x = f√η, and then

U_{1} = ∫_{0}^{∞} \frac{df}{1 + ηf^{2}} = \frac{1}{√η} ∫_{0}^{∞} \frac{dx}{1 + x^{2}} = \frac{πβ}{2γ} .  (148)

The integrand in (146) corresponds to L in (128) for m = j − 1 and k = 1, with −h^{2} replaced by f^{2} and δ replaced by η. From (131), the partial fraction expansion of the integrand in (144) is

\frac{(1 + f^{2})^{j−1}}{(1 + ηf^{2})^{j}} = \frac{1}{η^{j−1}} ∑_{q=0}^{j−1} \binom{j − 1}{q}\frac{(η − 1)^{q}}{(1 + ηf^{2})^{q+1}} .  (149)

Therefore, for j > 1,

U_{j} = ∫_{0}^{∞} \frac{1}{η^{j−1}} ∑_{q=0}^{j−1} \binom{j − 1}{q}\frac{(η − 1)^{q}}{(1 + ηf^{2})^{q+1}} df

= η^{1−j} ∑_{q=0}^{j−1} \binom{j − 1}{q}(η − 1)^{q} ∫_{0}^{∞} \frac{df}{(1 + ηf^{2})^{q+1}} .  (150)

Substitute x = f√η and then x = tan θ, and we get

∫_{0}^{∞} \frac{df}{(1 + ηf^{2})^{q+1}} = \frac{1}{√η} ∫_{0}^{∞} \frac{dx}{(1 + x^{2})^{q+1}} = \frac{β}{γ} ∫_{0}^{∞} \frac{dx}{(1 + x^{2})^{q+1}}
= \frac{β}{γ} ∫_{0}^{π/2} cos^{2q} θ dθ = \frac{β}{γ} ∫_{0}^{π/2} sin^{2q} θ dθ
= \frac{βs_{2q}}{γ} = \frac{β}{γ} · \frac{π}{2} · \frac{1}{2} · \frac{3}{4} · \frac{5}{6} ⋯ \frac{(2q − 1)}{(2q)} ,  (151)

in view of (33).

Thus, the integral H_{m,k} has been constructed in terms of elementary functions.

Then the integral in (99) can be computed (140) from H_{m,k} and H_{m−1,k+1}, to give explicitly the surface area for a biaxial ellipsoid of type B_{2m,2k}, in terms of elementary functions of its semi-axes β and γ.

C Capacity for biaxial ellipsoids

The indefinite integral

G_{m,k} =^{def} ∫ \frac{sin^{m} x}{cos^{k} x} dx .  (152)

has the reduction formula [Dwight, §452.9] for k ≠ 1:

G_{m,k} = \frac{sin^{m−1} x}{(k − 1) cos^{k−1} x} − (\frac{m − 1}{k − 1}) G_{m−2,k−2} .  (153)

Applied to the definite integral

I_{m,k} =^{def} ∫_{0}^{arcsin√μ} \frac{sin^{m} x}{cos^{k} x} dx .  (154)

this becomes

I_{m,k} = \frac{μ^{(m−1)/2}}{(k − 1)(1 − μ)^{(k−1)/2}} − (\frac{m − 1}{k − 1}) I_{m−2,k−2} .  (155)

C.1 Odd q = 2j + 1

Denote m = p + q − 3, and let q = 2j + 1 so that m = p + 2j − 2. The formula (125) requires the integral I_{m,q−1} = I_{p+2j−2,2j}. We shall compute I_{p,2}, and then apply the reduction formula (155) j − 1 times, to compute successively I_{p+2,4}, I_{p+4,6}, I_{p+6,8}, …, I_{p+2j−2,2j}.

As our starting value, we require

I_{p,2} = ∫_{0}^{w} \frac{sin^{p} x}{cos^{2} x} dx .  (156)

For a prolate spheroid with the larger semi-axis β having multiplicity p = 1, this becomes

I_{p,2} = ∫_{0}^{w} \frac{sin x dx}{cos^{2} x} = \frac{1}{cos x}|_{0}^{w} = \frac{1}{cos w} − 1 = \frac{1}{\sqrt{1 − μ}} − 1 . (157)

But for p > 1 we need to compute

I_{p−2,0} = ∫_{0}^{w} sin^{p−2} x dx ,  (158)

and then apply the reduction formula once to get I_{p,2}, and then proceed as above to compute I_{p+2j−2,2j}.

C.1.1 Integrating powers of sines

For non-negative integer h, define

ς_{h} = ∫_{0}^{w} sin^{h} x dx ,  (159)

where w = arcsin √μ. In particular,

ς_{0} = w = arcsin √μ,
ς_{1} = 1 − cos w = 1 − \sqrt{1 − μ} = \frac{μ}{1 + \sqrt{1 − μ}} .  (160)

For h > 1, integrating (159) by parts, we get the reduction formula

ς_{h} = (\frac{h − 1}{h}) ς_{h−2} − \frac{sin^{h−1} w cos w}{h}
= (\frac{h − 1}{h}) ς_{h−2} − \frac{μ^{(h−1)/2}\sqrt{1 − μ}}{h} .  (161)

Therefore,

ς_{2} = \frac{1}{2} (ς_{0} − √μ\sqrt{1 − μ}) = \frac{1}{2} (arcsin √μ − \sqrt{μ(1 − μ)}) ,
ς_{3} = \frac{1}{3} (2ς_{1} − μ\sqrt{1 − μ}) = \frac{1}{1×3} (2 − (2 + μ)\sqrt{1 − μ}) ,
ς_{4} = \frac{1}{4} (3ς_{2} − μ√μ\sqrt{1 − μ}) = \frac{1}{2×4} (3 arcsin √μ − (3 + 2μ)\sqrt{μ(1 − μ)}) ,
ς_{5} = \frac{1}{5} (4ς_{3} − μ^{2}\sqrt{1 − μ}) = \frac{1}{3×5} (8 − (8 + 4μ + 3μ^{2})\sqrt{1 − μ})) ,
ς_{6} = \frac{1}{6} (5ς_{4} − μ^{2}√μ\sqrt{1 − μ})
= \frac{1}{4×6} (15 arcsin √μ − (15 + 10μ + 4μ^{2})\sqrt{μ(1 − μ)}) , et cetera. (162)

Thus, the definite integral ς_{h} in (158) can be expressed in term of elementary functions, for all integers h ≥ 0.

C.2 Even q = 2j

Denote m = p + q − 3, and let q = 2j so that m = p + 2j − 3. The formula (125) requires the integral I_{m,q−1} = I_{p+2j−3,2j−1}. We shall compute I_{p−1,1}, and then apply the reduction formula (155) j − 1 times, to compute successively I_{p+1,3}, I_{p+3,5}, I_{p+5,7}, …, I_{p+2j−3,2j−1}.

As our starting value, we require

I_{p−1,1} = ∫_{0}^{w} \frac{sin^{p−1} x}{cos x} dx .  (163)

C.2.1 Odd p = 2g + 1

For odd p = 2g + 1, this becomes

I_{2g,1} = ∫_{0}^{w} sin^{2g} x \frac{dx}{cos x} = ∫_{0}^{w}(1 − cos^{2} x)^{g} \frac{dx}{cos x}

= ∫_{0}^{w} ∑_{r=0}^{g} \binom{g}{r}(− cos^{2} x)^{r} \frac{dx}{cos x} = ∑_{r=0}^{g} \binom{g}{r}(−1)^{r} ∫_{0}^{w} cos^{2r} x \frac{dx}{cos x}

= ∫_{0}^{w} \frac{dx}{cos x} + ∑_{r=1}^{g} \binom{g}{r}(−1)^{r} ∫_{0}^{w} cos^{2r−1} x dx .  (164)

Now,

∫_{0}^{w} \frac{dx}{cos x} = ∫_{0}^{w} sec x dx = log (\frac{1 + tan(w/2)}{1 − tan(w/2)}) ,  (165)

and

(\frac{1 + tan(w/2)}{1 − tan(w/2)})^{2} = (\frac{cos(w/2) + sin(w/2)}{cos(w/2) − sin(w/2)})^{2}

= \frac{1 + sin w}{1 − sin w} = \frac{1 + √μ}{1 − √μ} .  (166)

Therefore,

∫_{0}^{w} \frac{dx}{cos x} = \frac{1}{2} log (\frac{1 + √μ}{1 − √μ}) .  (167)

Next, substitute x = \frac{1}{2}π − y, and we get

∫_{0}^{w} cos^{k} x dx = ∫_{(π/2)−w}^{π/2} sin^{k} y dy

= ∫_{0}^{π/2} sin^{k} x dx − ∫_{0}^{(π/2)−w} sin^{k} x dx = s_{k} − λ_{k} . (168)

Here, s_{k} is the important special case (34) of ς_{k} with μ = 1.

For non-negative integer h, define

λ_{h} = ∫_{0}^{v} sin^{h} x dx ,  (169)

where

v = \frac{1}{2}π − w = arccos √μ = arcsin \sqrt{1 − μ} ,  (170)

so that

λ_{0} = v = arcsin \sqrt{1 − μ},
λ_{1} = 1 − cos v = 1 − √μ .  (171)

The reduction formula (161) for ς_{h} is to be converted by interchanging sin w (= √μ) and cos w (= \sqrt{1 − μ}), giving the reduction formula for λ_{h} with h > 1

λ_{h} = (\frac{h − 1}{h}) λ_{h−2} − \frac{cos^{h−1} w sin w}{h}
= (\frac{h − 1}{h}) λ_{h−2} − \frac{(1 − μ)^{(h−1)/2}√μ}{h} .  (172)

Staring from λ_{0} and λ_{1} as given in (170), this reduction formula can generate successively λ_{2}, λ_{3}, λ_{4}, · · · .

Thus, we compute the required integral

I_{2g,1} = \frac{1}{2} log (\frac{1 + √μ}{1 − √μ}) + ∑_{r=1}^{g} \binom{g}{r}(−1)^{r}(s_{2r−1} − λ_{2r−1}) .  (173)

And then we apply the reduction formula (155) j − 1 times, to compute successively I_{p+1,3}, I_{p+3,5}, I_{p+5,7}, …, I_{p+2j−3,2j−1}.

C.2.2 Even p = 2g As our starting value, we require

I_{p−1,1} = I_{2g−1,1} = ∫_{0}^{w} \frac{sin^{2g−1} x}{cos x} dx .  (174)

Then substituting c = cos x we get

I_{p−1,1} = ∫_{0}^{w} \frac{sin^{2(g−1)} x (sin x dx)}{cos x}
= ∫_{0}^{w} \frac{(1 − cos^{2} x)^{g−1} (sin x dx)}{cos x}
= ∫_{cos w}^{1}(1 − c^{2})^{g−1} \frac{dc}{c} = ∫_{√\overline{1−μ}}^{1} ∑_{r=0}^{g−1} \binom{g − 1}{r}(−c^{2})^{r} \frac{dc}{c}
= ∑_{r=0}^{g−1} \binom{g − 1}{r}(−1)^{r} ∫_{√\overline{1−μ}}^{1} c^{2r−1} dc
= ∫_{√\overline{1−μ}}^{1} \frac{dc}{c} + ∑_{r=1}^{g−1} \binom{g − 1}{r}(−1)^{r} ∫_{√\overline{1−μ}}^{1} c^{2r−1} dc
= − log \sqrt{1 − μ} + ∑_{r=1}^{g−1} \binom{g − 1}{r}(−1)^{r} \frac{1 − (\sqrt{1 − μ})^{2r}}{2r}
= \frac{1}{2} [log (\frac{1}{1 − μ}) + ∑_{r=1}^{g−1} \binom{g − 1}{r}(−1)^{r} \frac{1 − (1 − μ)^{r}}{r} ] . (175)

Having computed I_{p−1,1}, we then apply the reduction formula (155) j −1 times, to compute successively I_{p+1,3}, I_{p+3,5}, I_{p+5,7}, …, I_{p+2j−3,2j−1}.