Finding the Fastest Path Down
Imagine a ball rolling under gravity from point A to point B. Which path gets it there in the least time? Intuitively, a straight line looks like the shortest route — but shortest in distance isn’t the same as fastest in time. The answer, famously posed by Johann Bernoulli in 1696, is a curve called the cycloid. This is the birth story of the calculus of variations, and today we’ll solve it numerically and visually with Python.
The Problem, Mathematically
We drop a particle from rest at point $A = (0,0)$ and let it slide (frictionlessly, under gravity $g$) along some curve $y(x)$ to point $B = (x_1, y_1)$, where $y$ is measured downward as positive.
By energy conservation, the speed at height $y$ is:
$$
v = \sqrt{2gy}
$$
The time to traverse a small arc length $ds = \sqrt{1 + y’(x)^2}, dx$ is $dt = ds / v$. So the total descent time is the functional:
$$
T[y] = \int_0^{x_1} \frac{\sqrt{1 + y’(x)^2}}{\sqrt{2gy(x)}} , dx
$$
We want to find the function $y(x)$ that minimizes $T[y]$. Applying the Euler–Lagrange equation to this functional leads to the differential equation whose solution is a cycloid — the curve traced by a point on the rim of a rolling circle:
$$
x(\theta) = R(\theta - \sin\theta), \qquad y(\theta) = R(1 - \cos\theta)
$$
where $R$ is the rolling circle’s radius and $\theta$ is the rotation angle. Given the endpoint $B=(x_1,y_1)$, $R$ and the final angle $\theta_1$ are found by solving:
$$
\frac{y_1}{x_1} = \frac{1 - \cos\theta_1}{\theta_1 - \sin\theta_1}
$$
and the minimal descent time has the remarkably clean closed form:
$$
T_{\text{cycloid}} = \theta_1 \sqrt{\frac{R}{g}}
$$
Our Concrete Example
We’ll drop a particle from $A = (0, 0)$ to $B = (3.0,\ 2.0)$ meters (3 m horizontally, 2 m of drop), with $g = 9.81\ \text{m/s}^2$. We’ll compare three paths:
- A straight line — the “obvious” but wrong answer.
- A family of quadratic Bézier curves, whose shape we sweep over a grid of control points to search for the fastest one within that family.
- The true cycloid solution.
This lets us visually confirm that the cycloid beats every other candidate curve, and lets us plot the whole “time landscape” as a 3D surface — a nice, concrete picture of what calculus of variations is actually doing (searching an infinite-dimensional space of curves for a minimum).
Full Source Code
1 | # ============================================================ |
Code Walkthrough
Section 2 — the singularity trick. Every candidate curve starts at rest, so speed $v = \sqrt{2gy} \to 0$ as $t \to 0$, and the raw time integrand blows up like $1/\sqrt{t}$ near the start. This singularity is integrable (finite area), but numerically it can make scipy.integrate.quad slow or noisy. The substitution $t = u^2$ (so $dt = 2u,du$) cancels the $1/\sqrt{t}$ term algebraically, turning every curve’s integral into a smooth, fast-converging one. This is why the same descent_time() function works cleanly for the line, the Bézier curves, and the cycloid without any special-casing.
Section 3 — the straight line. This is the naive baseline: constant velocity direction, but it accelerates too slowly at first because it doesn’t drop steeply enough near the start.
Section 4 — the Bézier search. Instead of trying to derive the optimal curve, we brute-force search a 2-parameter family of quadratic Bézier curves (parametrized by a control point $(c_x, c_y)$) over a $25 \times 25$ grid, computing the descent time for each of the 625 candidate curves. This is fast — each quad() call is a few milliseconds — so the whole grid finishes in a couple of seconds in Colab. Conceptually, this loop is a discretized version of “searching the space of all curves for a minimum,” which is exactly what the calculus of variations does analytically via the Euler–Lagrange equation.
Section 5 — the cycloid. brentq solves the transcendental equation for $\theta_1$ (the ratio equation has a guaranteed sign change between $\theta \to 0^+$, where the ratio diverges, and $\theta \to 2\pi^-$, where it approaches 0). Once $\theta_1$ is known, $R$ follows directly, and we compute the descent time two independent ways — via numerical integration and via the closed-form formula $T = \theta_1\sqrt{R/g}$ — as a sanity check that they agree.
Section 7 — the 2D plot. All three curves are drawn on the same axes, with the y-axis inverted so “downward” reads visually as “down.” You should see the cycloid dip more steeply near A than the straight line (trading extra distance for extra early speed), then flatten out toward B.
Section 8 — the 3D plot. This is the most illuminating part. It plots the descent time $T$ as a surface over the 2D space of Bézier control points $(c_x, c_y)$ — literally the “cost landscape” that an optimizer (or evolution, or your own intuition) would need to descend to find the best curve shape. The red marker shows the grid’s discovered minimum. Because the true cycloid isn’t a member of this particular Bézier family, its time (printed in the console) will typically be slightly lower than the Bézier grid minimum — nicely demonstrating that the cycloid is the true global optimum across all curves, not just within one restricted family.
======================================================= BRACHISTOCHRONE RESULTS ======================================================= Target point B : (3.0, 2.0) m Cycloid radius R : 1.00133 m Cycloid final angle theta1 : 3.06878 rad ------------------------------------------------------- Straight line time : 1.15116 s Best Bezier grid time : 0.98209 s (cx=0.263, cy=1.600) Cycloid time (numeric integral): 0.98043 s Cycloid time (closed form) : 0.98043 s =======================================================
Reading the Results
Once you’ve pasted in your run, you should observe:
- $T_{\text{cycloid}} < T_{\text{bezier_min}} < T_{\text{line}}$ — the straight line is the slowest, the best-found Bézier curve is faster, and the true cycloid is fastest of all.
- The two cycloid time values (numeric integral vs. closed-form) should match to about 4–5 decimal places, confirming the numerical integration is accurate.
- On the 3D surface, the bowl-shaped landscape has a single, fairly broad minimum region — which is why gradient-based or grid-based search methods converge to it reliably, and also hints at why this class of variational problem tends to have a unique smooth solution rather than many competing local optima.
This little experiment is a nice hands-on echo of the historical event: Bernoulli’s challenge produced the same answer whether you attack it with pure analysis (Euler–Lagrange calculus) or, as we just did, with a grid search and numerical integration on a laptop three centuries later.
