The Brachistochrone problem

Posted on 30 March 2017

The Brachistochrone problem asks the question "what is the shape of the curve down which a bead sliding from rest and accelerated by gravity will slip, without friction, from one fixed point, $P_1$ to another $P_2$ in the least time?"

Its solution is well-known and given in various forms in many places on the internet, but is given in outline at the bottom of this post.

Paths between $P_1 = (0,0)$ and $P_2 = (1, 0.65)$:

enter image description here

Note that the path giving the shortest time (the cycloid) is not the shortest distance between the points and may have a minimum (i.e. the bead slides down and then up before reaching $P_2$, as shown by the following example for $P_2 = (1, 0.3)$.

enter image description here

The following Python program plots the Brachistochrone curve (an arc of a cycloid) and calculates the time of travel. It is also compared with three other paths: (i) a straight line, (ii) a circular path with a vertical tangent at $P_1$ and (iii) a parabola with vertical tangent at $P_1$.

import numpy as np
from scipy.optimize import newton
from scipy.integrate import quad
import matplotlib.pyplot as plt

# Acceleration due to gravity (m.s-2); final position of bead (m).
g = 9.81
x2, y2 = 1, 0.65

def cycloid(x2, y2, N=100):
    """Return the path of Brachistochrone curve from (0,0) to (x2, y2).

    The Brachistochrone curve is the path down which a bead will fall without
    friction between two points in the least time (an arc of a cycloid).
    It is returned as an array of N values of (x,y) between (0,0) and (x2,y2).

    """

    # First find theta2 from (x2, y2) numerically (by Newton-Rapheson).
    def f(theta):
        return y2/x2 - (1-np.cos(theta))/(theta-np.sin(theta))
    theta2 = newton(f, np.pi/2)

    # The radius of the circle generating the cycloid.
    R = y2 / (1 - np.cos(theta2))

    theta = np.linspace(0, theta2, N)
    x = R * (theta - np.sin(theta))
    y = R * (1 - np.cos(theta))

    # The time of travel
    T = theta2 * np.sqrt(R / g)
    print('T(cycloid) = {:.3f}'.format(T))
    return x, y, T

def linear(x2, y2, N=100):
    """Return the path of a straight line from (0,0) to (x2, y2)."""

    m = y2 / x2
    x = np.linspace(0, x2, N)
    y = m*x

    # The time of travel
    T = np.sqrt(2*(1+m**2)/g/m * x2)
    print('T(linear) = {:.3f}'.format(T))
    return x, y, T

def func(x, f, fp):
    """The integrand of the time integral to be minimized for a path f(x)."""

    return np.sqrt((1+fp(x)**2) / (2 * g * f(x)))

def circle(x2, y2, N=100):
    """Return the path of a circular arc between (0,0) to (x2, y2).

    The circle used is the one with a vertical tangent at (0,0).

    """

    # Circle radius
    r = (x2**2 + y2**2)/2/x2

    def f(x):
        return np.sqrt(2*r*x - x**2)
    def fp(x):
        return (r-x)/f(x)

    x = np.linspace(0, x2, N)
    y = f(x)

    # Calcualte the time of travel by numerical integration.
    T = quad(func, 0, x2, args=(f, fp))[0]
    print('T(circle) = {:.3f}'.format(T))
    return x, y, T

def parabola(x2, y2, N=100):
    """Return the path of a parabolic arc between (0,0) to (x2, y2).

    The parabola used is the one with a vertical tangent at (0,0).

    """

    c = (y2/x2)**2

    def f(x):
        return np.sqrt(c*x)
    def fp(x):
        return c/2/f(x)

    x = np.linspace(0, x2, N)
    y = f(x)

    # Calcualte the time of travel by numerical integration.
    T = quad(func, 0, x2, args=(f, fp))[0]
    print('T(parabola) = {:.3f}'.format(T))
    return x, y, T

# Plot a figure comparing the four paths.
fig, ax = plt.subplots()

for curve in ('cycloid', 'circle', 'parabola', 'linear'):
    x, y, T = globals()[curve](x2, y2)
    ax.plot(x, y, lw=4, alpha=0.5, label='{}: {:.3f} s'.format(curve, T))
ax.legend()

ax.set_xlabel('$x$')
ax.set_xlabel('$y$')
ax.set_xlim(0, 1)
ax.set_ylim(0.8, 0)
plt.savefig('brachistochrone.png')
plt.show()

Appendix A: Solution the the Brachistochrone problem

We choose $P_1 = (0,0)$ and set $P_2 = (x_2, y_2)$, with the $y-$axis pointing in the same direction as the acceleration due to gravity, $\boldsymbol{g}$. Then we seek to minimize $$ T = \int_{P_1}^{P_2}\frac{\mathrm{d}s}{v}, $$ where the arc length element along the curve we seek is $\mathrm{d}s = \sqrt{\mathrm{d}x^2 + \mathrm{d}y^2}$ and the speed along the curve, $v$ satisfies $\frac{1}{2}mv^2 = mgy$ by conservation of energy. Then the functional to be minimized is $$ T = \int_{P_1}^{P_2}\sqrt{\frac{1+y'^2}{2gy}}\,\mathrm{d}x, $$ where $y' = \mathrm{d}y/\mathrm{d}x$. The function to be varied in this minimization, $$ F(x,y,y') = \sqrt{\frac{1+y'^2}{2gy}} $$ does not explicitly depend on $x$ so the Beltrami identity applies: $$ F - y'\frac{\partial F}{\partial y'} = C = \mathrm{constant}. $$ Now, $$ \frac{\partial F}{\partial y'} = \frac{y'}{\sqrt{2gy}\sqrt{1+y'^2}} $$ and substitution into the Beltrami identity gives $$ y(1+y'^2) = \frac{1}{2gC^2} = 2R. $$ This differential equation has a parametric solution in the form of a cycloid: \begin{align*} x(\theta) &= R(\theta - \sin\theta),\\ y(\theta) &= R(1 - \cos\theta). \end{align*} The parameter $\theta$ varies from $0$ at $P_1 = (0,0)$ to $\theta_2$ at $P_2 = (x_2, y_2)$ which may be found numerically as the solution to $$ \frac{y_2}{x_2} = \frac{1-\cos\theta_2}{\theta_2 - \sin\theta_2} $$ Substitution back into the second parametric equation then gives $$ R = \frac{y_2}{1-\cos\theta_2}. $$

Appendix B: Paths and times for other curves

Straight line

The straight line between $P_1 = (0,0)$ and $P_2 = (x_2,y_2)$ has the equation $y = mx$ with $m=y_2/x_2$. The time taken is $$ T = \int_{P_1}^{P_2}\sqrt{\frac{1+y'^2}{2gy}}\,\mathrm{d}x = \int_0^{x_2} \sqrt{\frac{1+m^2}{2gmx}}\,\mathrm{d}x = \sqrt{\frac{(1+m^2)x_2}{gm}}. $$

Circular arc

The circle through $P_1 = (0,0)$ and $P_2 = (x_2,y_2)$ with vertical tangent at $P_1$ has the formula $y^2 = r^2 - (x-r)^2$ with radius $$ r = \frac{x_2^2 + y_2^2}{2x_2}. $$ The integral giving $T$ may be evaluated numerically using $y = \sqrt{2xr-x^2}$ and $y' = (r-x)/y$.

Parabolic arc

The parabola through $P_1 = (0,0)$ and $P_2 = (x_2,y_2)$ with vertical tangent at $P_1$ has the formula $y = \sqrt{cx}$ with $c = (y_2/x_2)^2$. The integral giving $T$ may be evaluated numerically using this and $y' = c/(2y)$.