Classic variational problems: brachistochrone, catenary and geodesics

Three problems shaped the calculus of variations, and each still teaches something. The brachistochrone asks for the track of quickest descent and shows that the answer to a variational problem can be surprising. The catenary asks for the shape of a hanging chain and introduces constraints handled by a Lagrange multiplier. Geodesics, the shortest paths on a curved surface, connect variational calculus to geometry and, through general relativity, to gravity itself. This lesson works through all three with the tools of lesson 2 and checks each result numerically.

The brachistochrone

In June 1696 Johann Bernoulli published a challenge in the journal Acta Eruditorum: given two points at different heights, not one directly below the other, find the curve along which a bead, sliding without friction under gravity and starting from rest, travels from the upper point to the lower in the least time. He called it the brachistochrone, from the Greek for shortest time. Solutions came from Johann himself, his brother Jakob, Leibniz, l'Hôpital and an anonymous English author who was quickly recognized as Newton; Johann is said to have remarked that he knew the lion by its claw.

Set up coordinates with the start at the origin and \(y\) measured downward. Energy conservation gives the speed at depth \(y\) as \(v = \sqrt{2gy}\), independent of the shape of the track. The time to travel an element of arc \(ds = \sqrt{1 + y'^2}\,dx\) is \(ds/v\), so

The integrand has no explicit \(x\), so the Beltrami identity applies. With \(F = \sqrt{1 + y'^2}/\sqrt y\) (dropping the constant \(\sqrt{2g}\)), \(\partial F/\partial y' = y'/(\sqrt y\sqrt{1 + y'^2})\), and

for some constant \(2R\). This first-order equation is solved by the substitution \(y = R(1 - \cos\theta)\). Then \(1 + y'^2 = 2R/y = 2/(1 - \cos\theta)\), so \(y'^2 = (1 + \cos\theta)/(1 - \cos\theta) = \cot^2(\theta/2)\). Since \(dy = R\sin\theta\,d\theta\), \(dx = dy/y' = R\sin\theta\tan(\theta/2)\,d\theta = R(1 - \cos\theta)\,d\theta\). Integrating from the origin,

These are the equations of a cycloid, the curve traced by a point on the rim of a wheel of radius \(R\) rolling along the underside of the line \(y = 0\). The constant \(R\) is fixed by requiring the cycloid to pass through the end point. The time along it is particularly simple: since \(ds = 2R\sin(\theta/2)\,d\theta\) and \(v = \sqrt{2gR(1 - \cos\theta)} = 2\sqrt{gR}\sin(\theta/2)\), \(dt = ds/v = \sqrt{R/g}\,d\theta\), and the time to reach parameter \(\theta\) is \(\theta\sqrt{R/g}\).

Why is the cycloid quickest? It drops steeply at the start, where the bead is slow, so that it gains speed early and covers the long middle section quickly. A straight line is the shortest path but wastes the early part of the trip at low speed. The trade-off between length and speed is exactly what the variational principle balances.

The tautochrone

The cycloid has a second remarkable property, discovered by Christiaan Huygens in 1659. A bead released from rest anywhere on an inverted cycloid reaches the bottom in the same time, \(\pi\sqrt{R/g}\), no matter where it starts. The curve is a tautochrone, an equal-time curve. The reason is that along the cycloid the restoring force is exactly proportional to the arc length from the bottom, so the motion is simple harmonic for any amplitude. Huygens used this to design a pendulum clock whose bob was guided by cycloidal cheeks at the top of the string, making its period independent of amplitude; friction in the cheeks made the idea less useful in practice than in theory.

import math
g = 9.81

# 1. Brachistochrone: slide from (0, 0) to (x1, y1) with y measured downward. Compare travel times.
def travel_time(path, n=200000):
    # path(s) -> (x, y) for s in [0, 1]; time = sum ds / v with v = sqrt(2 g y)
    # substitute s = u^2 so that the 1/sqrt(y) singularity at the start becomes harmless
    t = 0.0; x0, y0 = path(0.0)
    for i in range(1, n + 1):
        x1_, y1_ = path((i / n) ** 2)
        ds = math.hypot(x1_ - x0, y1_ - y0); ym = 0.5 * (y0 + y1_)
        t += ds / math.sqrt(2 * g * ym) if ym > 0 else 0.0
        x0, y0 = x1_, y1_
    return t

X, Y = math.pi, 2.0                    # end point; this is the cycloid with R = 1 and theta from 0 to pi
R = 1.0
cycloid = lambda s: (R * (s * math.pi - math.sin(s * math.pi)), R * (1 - math.cos(s * math.pi)))
line = lambda s: (X * s, Y * s)
parab = lambda s: (X * s, Y * (2 * s - s * s))   # steep start, horizontal at the end
print(f"straight line : {travel_time(line):.4f} s (exact {math.sqrt(2 * (X * X + Y * Y) / (g * Y)):.4f})")
print(f"parabola      : {travel_time(parab):.4f} s")
print(f"cycloid       : {travel_time(cycloid):.4f} s (exact pi sqrt(R/g) = {math.pi * math.sqrt(R / g):.4f})")

# Tautochrone: from any starting point on the cycloid the time to the bottom is the same
for th0 in (0.3, 1.0, 2.0, 3.0):
    # parametrize from theta0 to pi; speed from energy, v = sqrt(2 g (y - y0))
    n = 200000; t = 0.0
    y0 = R * (1 - math.cos(th0))
    for i in range(n):
        u = (i + 0.5) / n                        # theta = theta0 + (pi - theta0) u^2 removes the singularity
        a = th0 + (math.pi - th0) * u * u
        ds = 2 * R * math.sin(a / 2) * (math.pi - th0) * 2 * u / n
        t += ds / math.sqrt(2 * g * (R * (1 - math.cos(a)) - y0))
    print(f"start at theta0 = {th0:.1f}: time to the lowest point = {t:.4f} s")

# 2. Catenary: a 12 m chain hung between supports 10 m apart. Solve 2a sinh(d/(2a)) = L for a.
d, L = 10.0, 12.0
lo, hi = 0.5, 100.0
for _ in range(200):
    a = 0.5 * (lo + hi)
    lo, hi = (a, hi) if 2 * a * math.sinh(d / (2 * a)) > L else (lo, a)
sag = a * (math.cosh(d / (2 * a)) - 1)
print(f"catenary parameter a = {a:.4f} m, sag = {sag:.4f} m")
# a parabola with the same span and sag has a different shape and length
para_len = 0.0; n = 100000
for i in range(n):
    x = -d / 2 + d * (i + 0.5) / n
    slope = 8 * sag * x / d ** 2
    para_len += math.sqrt(1 + slope * slope) * d / n
print(f"parabola with the same span and sag: length {para_len:.4f} m")
print(f"height above the lowest point, 2.5 m from the middle: catenary {a * (math.cosh(d / (4 * a)) - 1):.4f} m, "
      f"parabola {sag / 4:.4f} m")

# 3. Geodesics on a sphere: two points at latitude 50 degrees, 90 degrees of longitude apart
Re = 6371.0; lat = math.radians(50); dlon = math.radians(90)
along_parallel = Re * math.cos(lat) * dlon
cosc = math.sin(lat) ** 2 + math.cos(lat) ** 2 * math.cos(dlon)
great = Re * math.acos(cosc)
top = math.degrees(math.atan(math.tan(lat) / math.cos(dlon / 2)))
print(f"along the parallel: {along_parallel:.0f} km; great circle: {great:.0f} km; "
      f"great circle reaches latitude {top:.1f} degrees")

straight line : 1.1890 s (exact 1.1890)
parabola      : 1.0461 s
cycloid       : 1.0030 s (exact pi sqrt(R/g) = 1.0030)
start at theta0 = 0.3: time to the lowest point = 1.0030 s
start at theta0 = 1.0: time to the lowest point = 1.0030 s
start at theta0 = 2.0: time to the lowest point = 1.0030 s
start at theta0 = 3.0: time to the lowest point = 1.0030 s
catenary parameter a = 4.6954 m, sag = 2.9234 m
parabola with the same span and sag: length 11.9529 m
height above the lowest point, 2.5 m from the middle: catenary 0.6814 m, parabola 0.7309 m
along the parallel: 6433 km; great circle: 6012 km; great circle reaches latitude 59.3 degrees

The program integrates \(\int ds/v\) along each track with a change of variable that removes the infinite integrand at the start, where \(v = 0\). The straight line takes 1.189 s, matching the exact formula \(\sqrt{2(X^2 + Y^2)/(gY)}\); the parabola 1.046 s; and the cycloid 1.003 s, which is \(\pi\sqrt{R/g}\) for \(R = 1\ \mathrm m\). Starting the bead from rest at four different points on the cycloid, from barely above the bottom to nearly the top, gives the same 1.0030 s each time: the tautochrone property, checked to four figures.

The catenary

A uniform flexible chain of length \(\ell\) and mass per unit length \(\mu\) hangs between two supports. It settles into the shape of least potential energy, subject to the constraint that its length is \(\ell\). The potential energy is \(\mu g\int y\,ds = \mu g\int y\sqrt{1 + y'^2}\,dx\), and the length is \(\int\sqrt{1 + y'^2}\,dx\). With a Lagrange multiplier \(\lambda\), we make stationary

Write \(u = y - \lambda\); the integrand \(u\sqrt{1 + u'^2}\) has no explicit \(x\), and the Beltrami identity gives \(u\sqrt{1 + u'^2} - u u'^2/\sqrt{1 + u'^2} = u/\sqrt{1 + u'^2} = a\), a constant. So \(u'^2 = u^2/a^2 - 1\), whose solution symmetric about \(x = 0\) is

This is the catenary, from the Latin for chain. The constant \(a\) has a physical meaning: the horizontal component of tension is the same everywhere along the chain and equals \(\mu g a\). For supports at equal heights a distance \(d\) apart, the length condition is \(2a\sinh(d/2a) = \ell\), which must be solved numerically for \(a\).

Galileo had suggested that a hanging chain forms a parabola. It does not, although the two are close when the sag is small. The parabola is the shape of a cable carrying a load spread uniformly along the horizontal, like the deck of a suspension bridge, while the catenary carries a load spread uniformly along its own length. The correct shape was found in 1691 by Jakob Bernoulli's rivals Leibniz, Huygens and Johann Bernoulli, in answer to Jakob's challenge. The program solves \(2a\sinh(5/a) = 12\) for a 12 m chain on supports 10 m apart, giving \(a = 4.695\ \mathrm m\) and a sag of 2.923 m. A parabola with the same span and sag is 11.953 m long, shorter than the chain.

Turned upside down, a catenary is the ideal shape for an arch made of blocks that can carry compression but no bending: every part of the arch is in pure compression along its own length. Robert Hooke stated this in 1675 in an anagram, decoded after his death as: as hangs the flexible line, so but inverted will stand the rigid arch. The Gateway Arch in St. Louis is a weighted catenary, and Antoni Gaudí designed vaults by hanging chains loaded with small bags of shot and photographing them upside down.

Minimal surfaces of revolution

A closely related problem asks for the curve \(y(x)\) between two points which, rotated about the \(x\) axis, sweeps out the smallest area. The area is \(2\pi\int y\sqrt{1 + y'^2}\,dx\), the same integrand as the chain without the multiplier, so the answer is again \(y = a\cosh(x/a)\): the surface is a catenoid. Lesson 2's program checked that the Beltrami expression is constant on such a curve. A soap film stretched between two coaxial rings forms a catenoid; if the rings are pulled too far apart, no catenoid fits, and the film snaps into two flat discs. That is a reminder that extremals need not exist for all boundary conditions.

Geodesics

On a curved surface, the shortest path between two points is called a geodesic. Describe the surface by two coordinates \((u, v)\) and its geometry by the line element \(ds^2 = E\,du^2 + 2F\,du\,dv + G\,dv^2\), the first fundamental form of MATH-430. The length of a curve \(v(u)\) is \(\int\sqrt{E + 2Fv' + Gv'^2}\,du\), and geodesics satisfy its Euler–Lagrange equation.

On a sphere of radius \(R\), with latitude \(\lambda\) and longitude \(\phi\), \(ds^2 = R^2(d\lambda^2 + \cos^2\lambda\,d\phi^2)\). A cleaner method treats the path as parametrized by a parameter \(s\) and uses the integrand \(\tfrac12(\dot\lambda^2 + \cos^2\lambda\,\dot\phi^2)\), whose extremals are geodesics traversed at constant speed. Longitude is cyclic, so \(\cos^2\lambda\,\dot\phi\) is constant; this is Clairaut's relation, the geometric version of angular momentum conservation. The solutions are great circles, the intersections of the sphere with planes through its centre.

Two consequences surprise people on first meeting. A parallel of latitude other than the equator is not a geodesic, so the shortest route between two cities at the same latitude bulges toward the pole. The program computes two points at 50 degrees north and 90 degrees of longitude apart: the route along the parallel is 6433 km, the great circle 6012 km, and the great circle reaches 59.3 degrees north at its midpoint. This is why flights from western Canada or the northern United States to Europe pass over Greenland. Second, geodesics can stop being shortest: on a sphere, a great-circle arc longer than half the circumference is a geodesic but not the shortest path. This is the geometric analogue of the action being a saddle point for long enough times in lesson 1.

Geodesics and mechanics

There is a precise sense in which free motion is geodesic motion. A particle sliding without friction on a smooth surface, with no other forces, has Lagrangian \(L = \tfrac12 m\,(E\dot u^2 + 2F\dot u\dot v + G\dot v^2)\), which is the geodesic integrand. Its speed is constant, and it moves along a geodesic of the surface. Jacobi showed in 1837 that even motion in a potential can be described this way, by choosing a modified metric whose geodesics are the trajectories at a given energy. Einstein's general relativity, the subject of PHYS-450, takes this idea to its conclusion: gravity is not a force at all, and freely falling bodies follow geodesics of curved spacetime.