An ellipsoid is the three-dimensional figure bounded by the surface described by the equation
$$
\frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1,
$$
where $a$, $b$ and $c$ are the semi-principal axes. If $a=b=c$, the ellipsoid is a sphere. The volume of an ellipsoid has a simple form,
$$
V = \frac{4}{3}\pi abc.
$$
There is no closed formula for the surface area of a general ellipsoid, but it may be expressed in terms of incomplete elliptic integrals of the first and second kinds, $K(\phi, k)$ and $E(\phi, k)$:
$$
S = 2\pi c^2 + \frac{2\pi ab}{\sin\phi} \left( K(\phi, k^2)\cos^2\phi + E(\phi, k^2)\sin^2\phi \right),
$$
where
$$
\cos\phi = \frac{c}{a}, \quad k = \frac{a\sqrt{b^2-c^2}}{b\sqrt{a^2-c^2}}
$$
and the co-ordinate system has been chosen such that $a \ge b \ge c$.
Define a function, ellipsoid_surface to calculate the surface area of a general ellipsoid, and compare the results for different-shaped ellipsoids with the following approximate formula:
$$
\begin{align*}
& S \approx 2\pi c^2 + 2\pi abr\left( 1 - \frac{b^2-c^2}{6b^2}r^2 \left( 1 - \frac{3b^2 + 10c^2}{56b^2}r^2\right) \right),\\
\mathrm{where} \quad& r = \frac{\phi}{\sin\phi}.
\end{align*}
$$
Solution P8.1.4
Here is one solution.
import numpy as np
from scipy.special import ellipkinc, ellipeinc
def ellipsoid_surface(a, b, c):
# We must have a >= b >= c, so order the arguments
a, b, c = reversed(sorted([a, b, c]))
if a == c:
# The ellipsoid is a sphere
return 4 * np.pi * a**2
cosphi = c / a
phi = np.arccos(cosphi)
cosphi2 = cosphi**2
sinphi2 = 1 - cosphi2
k2 = a**2 * (b**2 - c**2) / b**2 / (a**2 - c**2)
# Round to 8 decimal places to avoid an unfortunate numerical instability
k2 = round(k2, 8)
return 2 * np.pi * c**2 + 2 * np.pi * a * b / np.sin(phi) * (
ellipkinc(phi, k2) * cosphi2 + ellipeinc(phi, k2) * sinphi2
)
def ellipsoid_surface_approx(a, b, c):
# We must have a >= b >= c, so order the arguments
a, b, c = reversed(sorted([a, b, c]))
cosphi = c / a
phi = np.arccos(cosphi)
r = phi / np.sin(phi)
return 2 * np.pi * c**2 + 2 * np.pi * a * b * r * (
1
- (b**2 - c**2)
/ 6
/ b**2
* r**2
* (1 - (3 * b**2 + 10 * c**2) / 56 / b**2 * r**2)
)
# Fix a = 1 and construct an array of values for b = 0.02, 0.04, ... 0.98
a = 1
db = 0.02
nb = int((1 - 2 * db) / db + 1)
bgrid = np.linspace(db, 1 - db, nb)
# de will hold the relative error for different combinations of (b, c)
# de is symmetric, so only calculate for b >= c
de = np.zeros((nb, nb))
for i, b in enumerate(bgrid):
for j in range(0, i + 1):
c = bgrid[j]
Sexact = ellipsoid_surface(a, b, c)
Sapprox = ellipsoid_surface_approx(a, b, c)
de[i, j] = np.abs(Sexact - Sapprox) / Sexact
# Fill in the symmetric upper triangle of the de array
de = de + de.T - np.diag(de.diagonal())
maxde = np.max(np.abs(de))
print(f"Maximum percentage error in approximation: {maxde * 100:.2f}")
import matplotlib.pyplot as plt
plt.pcolor(bgrid, bgrid, de)
plt.xlabel(r"$b$")
plt.ylabel(r"$c$")
plt.colorbar()
plt.show()
A couple of points about the above code should be explained. Firstly, the algorithm requires $a, b, c$ to satisfy $a \ge b \ge c$ so we order them that way and check that the ellipsoid isn't a sphere (the definition of k2 will run into problems if it is since its denominator will be zero).
There turns out to be a curious numerical instability in SciPy's implementation of the incomplete elliptic integrals for some particular exact combinations of $\phi$ and $k^2$, so we round k2 to 8 decimal places (more than enough) to skirt it.
To analyse the approximate ellipsoid surface area formula we need to investigate the results produce for different combinations of $a$, $b$ and $c$. There are really only two degrees of freedom here so we can pick one of them, $a$, to be anything we like and we choose $a=1$. Then we consider values of $b$ satisfying $0 < b < 1$, at some convenient resolution (steps of 0.02) and values of $c$ at the same resolution, with $b >= c$ (the volume clearly doesn't depent on which axis is called $b$ and which is called $c$). We can fill in the symmetric triangle of the difference array by adding it to the transpose of itself and subtracting the diagonal. The relative error is plotted below.
Error in the approximate formula for an ellipsoid's surface area.
