Case Study 1 — How Long Is an Orbit? Measuring an Ellipse with Simpson's Rule

DataField.Dev

Case Study 1 — How Long Is an Orbit? Measuring an Ellipse with Simpson's Rule

Field: Astronomy and orbital mechanics Calculus used: Arc length set-up (preview of Chapter 18), trigonometric integrals (§16.2), and Simpson's rule with its error bound (§16.6) The punchline: the perimeter of an ellipse is the integral that cannot be done by any technique in this chapter — and that historical failure is exactly what made numerical integration indispensable.

A question that should be easy

A satellite traces an ellipse around the Earth; a planet traces an ellipse around the Sun. Mission planners routinely need one very ordinary-sounding number: how far does the body travel in one full orbit? That is just the perimeter of an ellipse. The area of an ellipse is the tidy $\pi a b$. Surely the perimeter is just as clean?

It is not. The perimeter of an ellipse is one of the most famous integrals in mathematics that has no elementary antiderivative at all — and the people who discovered that, in the seventeenth and eighteenth centuries, were precisely the astronomers trying to predict where planets would be. This case study walks the whole road: set up the integral, watch every analytic technique from this chapter fail, and then compute the answer to seven correct digits with Simpson's rule (§16.6).

Take a concrete orbit. Model it as the ellipse $$\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1, \qquad a = 5, \quad b = 3$$ (in units of thousands of kilometers, say), with semi-major axis $a = 5$ and semi-minor axis $b = 3$. We want its perimeter $P$.

Setting up the arc-length integral

Chapter 18 will develop arc length in full, but the idea is already within reach. Parametrize the ellipse by the angle $t$: $$x(t) = a\cos t, \qquad y(t) = b\sin t, \qquad 0 \le t \le 2\pi.$$ A tiny step $dt$ moves the point by $dx = -a\sin t\,dt$ and $dy = b\cos t\,dt$, so the length of that tiny step is $$ds = \sqrt{dx^2 + dy^2} = \sqrt{a^2\sin^2 t + b^2\cos^2 t}\;dt.$$ The whole perimeter is the accumulation of all those steps. By symmetry the four quarter-arcs are equal, so we integrate one quarter and multiply by four: $$P = 4\int_0^{\pi/2} \sqrt{a^2\sin^2 t + b^2\cos^2 t}\;dt.$$ With $a = 5$, $b = 3$ this is $$P = 4\int_0^{\pi/2} \sqrt{25\sin^2 t + 9\cos^2 t}\;dt.$$

Now we try to evaluate it the honest way.

Every analytic technique fails — and that is the lesson

Can we simplify the radicand into a single trig power? Use $\sin^2 t = 1 - \cos^2 t$ to collect terms: $$25\sin^2 t + 9\cos^2 t = 25 - 16\cos^2 t = 25\left(1 - \tfrac{16}{25}\cos^2 t\right).$$ So $$P = 20\int_0^{\pi/2}\sqrt{1 - e^2\cos^2 t}\;dt, \qquad e^2 = \frac{16}{25}, \quad e = \frac45,$$ where $e$ is the orbit's eccentricity (how far from circular it is; $e = 0$ is a circle). This is genuinely tidy — but the square root of $1 - e^2\cos^2 t$ is the wall.

Trigonometric integrals (§16.2)? The parity tricks need integer powers of sine and cosine, not a square root wrapped around them. No identity peels off a clean differential.
Trigonometric substitution (§16.4)? That technique attacks $\sqrt{a^2 - x^2}$, $\sqrt{a^2 + x^2}$, or $\sqrt{x^2 - a^2}$ — a radical over an algebraic expression in $x$. Here the radical sits over a trigonometric expression; there is no substitution that collapses it.
Partial fractions (§16.5)? There is no rational function in sight.

The integral $\displaystyle\int_0^{\pi/2}\sqrt{1 - e^2\cos^2 t}\,dt$ is a complete elliptic integral of the second kind, written $E(e)$. In the 1690s Jakob Bernoulli, working on exactly this curve, conjectured it could not be reduced to the known functions; in 1833 Liouville's theory confirmed that no elementary antiderivative exists. The ellipse perimeter is the historical origin of the phrase "elliptic integral," and the reason a whole branch of nineteenth-century mathematics grew up around functions that are defined by integrals rather than formulas. The §16.4 callout flagged this; here we live it.

So we do the mathematically correct thing when the analytic toolkit is exhausted: we compute numerically.

Simpson's rule, by hand-scale arithmetic

We approximate the quarter-arc $$Q = \int_0^{\pi/2}\sqrt{25\sin^2 t + 9\cos^2 t}\;dt,$$ with Simpson's rule (§16.6). Take $n = 4$ subintervals, so $h = \dfrac{\pi/2 - 0}{4} = \dfrac{\pi}{8} \approx 0.3927$, and nodes $t_0 = 0,\ t_1 = \tfrac\pi8,\ t_2 = \tfrac\pi4,\ t_3 = \tfrac{3\pi}8,\ t_4 = \tfrac\pi2$. Let $f(t) = \sqrt{25\sin^2 t + 9\cos^2 t}$. Evaluating:

$t$	$f(t)$
$0$	$3.0000$
$\pi/8$	$3.3341$
$\pi/4$	$4.1231$
$3\pi/8$	$4.7554$
$\pi/2$	$5.0000$

Apply the Simpson weights $1, 4, 2, 4, 1$: $$Q \approx \frac{h}{3}\big[f_0 + 4f_1 + 2f_2 + 4f_3 + f_4\big] = \frac{0.3927}{3}\big[3 + 4(3.3341) + 2(4.1231) + 4(4.7554) + 5\big].$$ The bracket is $3 + 13.3364 + 8.2462 + 19.0216 + 5 = 48.6042$, so $$Q \approx 0.13090 \times 48.6042 \approx 6.3824, \qquad P = 4Q \approx 25.530.$$

How good is that? The true value (computed below from the elliptic integral) is $P = 25.526999$. Even with only four panels, Simpson's rule lands within $0.003$ — about four correct digits — on an integral that no amount of algebra could ever close. That is the whole argument of this chapter in one number.

Pushing the precision in Python

# Perimeter of the ellipse x^2/25 + y^2/9 = 1 via Simpson's rule,
# checked against the exact complete elliptic integral E.
import numpy as np
from scipy.integrate import simpson, quad
from scipy.special import ellipe

a, b = 5.0, 3.0
f = lambda t: np.sqrt(a**2*np.sin(t)**2 + b**2*np.cos(t)**2)

# Exact perimeter: P = 4a * E(e^2), with e^2 = 1 - (b/a)^2
m = 1 - (b/a)**2
P_exact = 4*a*ellipe(m)
print(f"exact perimeter : {P_exact:.7f}")          # 25.5269989

for n in [4, 8, 16]:
    t = np.linspace(0, np.pi/2, n + 1)
    P = 4*simpson(f(t), x=t)
    print(f"Simpson n={n:>2}: P = {P:.7f}, error {abs(P - P_exact):.1e}")

# Output:
# exact perimeter : 25.5269989
# Simpson n= 4: P = 25.5294891, error 2.5e-03
# Simpson n= 8: P = 25.5270021, error 3.3e-06
# Simpson n=16: P = 25.5269990, error 1.4e-07

Doubling $n$ from $4$ to $8$ cut the error from $2.5\times10^{-3}$ to $3.3\times10^{-6}$ — a factor of roughly $760$. The $1/n^4$ error bound of §16.6 predicts a factor of $2^4 = 16$ at minimum; the integrand here is so smooth and periodic that Simpson does dramatically better than the worst-case guarantee. This is the happy face of the error bound: it is a ceiling on the error, and well-behaved integrands sail far below it.

Why the error bound is not just bookkeeping

For a real mission you cannot peek at the exact answer — there is none in closed form. The §16.6 error bound is what lets you certify precision without it: $$|Q - S_n| \le \frac{(b-a)^5}{180\,n^4}\max_{[a,b]}\big|f^{(4)}\big|.$$ Estimate $\max|f^{(4)}|$ over $[0,\pi/2]$ (numerically, it is around $30$ for this integrand), and the bound at $n = 8$ comes out near $10^{-4}$ — comfortably larger than the true $3\times10^{-6}$ error, exactly as a bound should be. A navigation engineer reports "$P = 25.5270$, guaranteed to $\pm 10^{-4}$" and can defend that claim with a single inequality, no exact answer required. That guarantee is the practical payoff of the whole numerical-integration section.

Connections

§16.2 (trigonometric integrals): the same parity-and-identity machinery handles the radicand's cousins — but the square root here is exactly what pushes the problem outside the analytic toolkit.
§16.4 (trigonometric substitution): flagged the ellipse as the integral that resists trig substitution and forces numerical methods; this case study is that flag cashed in.
§16.6 (numerical integration): trapezoidal and Simpson's rules plus their error bounds are the tools that actually deliver the number.
Chapter 18 (applications of integration): arc length is developed there in full generality; this is a preview applied to a curve that has no closed-form length.
Chapter 24 (applications of series): the same $E(e)$ can be expanded as a series in $e^2$, giving Ramanujan's celebrated perimeter approximations — a second route to the same number.

Discussion Questions

Why does the area of an ellipse have a clean formula ($\pi a b$) while the perimeter does not? What is structurally different about the two integrals? (Hint: compare $\int_0^a \sqrt{1 - x^2/a^2}\,dx$ with $\int \sqrt{a^2\sin^2 t + b^2\cos^2 t}\,dt$.)
As $b \to a$ the ellipse becomes a circle. Show that $e \to 0$ and $P \to 2\pi a$, recovering the familiar circumference. What happens to the difficulty of the integral in that limit?
Trapezoidal vs. Simpson here. If you run the trapezoidal rule on this periodic, symmetric integrand, it converges astonishingly fast — sometimes matching Simpson. Why might the trapezoidal rule be unusually accurate for smooth periodic integrands? (This is the basis of "spectral accuracy" in numerical analysis.)
A nearly parabolic comet orbit has $e$ very close to $1$. Why does the integrand become harder to integrate numerically as $e \to 1$, and how would you respond — more nodes, adaptive quadrature (quad), or a substitution?
Mission engineers often want perimeter to $\pm 1$ meter on an orbit of millions of kilometers. Using the Simpson error bound, roughly how large must $n$ be? Is that feasible to evaluate?

Your Turn — Mini-Project

Pick an orbit (Earth: $a \approx 6{,}378$ km altitude scaling, $e \approx 0.0167$ for Earth's orbit around the Sun; or invent one). Then:

Set up the perimeter integral $P = 4a\int_0^{\pi/2}\sqrt{1 - e^2\cos^2 t}\,dt$.
Implement Simpson's rule yourself (do not call scipy.integrate.simpson).
Tabulate the estimate and its error for $n = 2, 4, 8, 16, 32$ against 4*a*scipy.special.ellipe(e**2).
Plot $\log(\text{error})$ versus $\log n$ and confirm the slope is about $-4$.
Compare your best estimate to Ramanujan's approximation $P \approx \pi\big[3(a+b) - \sqrt{(3a+b)(a+3b)}\big]$. How many digits does Ramanujan's formula get for free?

Case Study 1 — How Long Is an Orbit? Measuring an Ellipse with Simpson's Rule

A question that should be easy

Setting up the arc-length integral

Every analytic technique fails — and that is the lesson

Simpson's rule, by hand-scale arithmetic

Pushing the precision in Python

Why the error bound is not just bookkeeping

Connections

Discussion Questions

Your Turn — Mini-Project

Further Reading