Posted by: christian on 30 Mar 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)$:
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)$.
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)$.
Share on Twitter Share on Facebook
Comments
Comments are pre-moderated. Please be patient and your comment will appear soon.
Matt 5 years, 9 months ago
Wonderful. Thanks for sharing this. I'm very happy to have found your website.
Link | ReplyIf y2 > x2 is a brachistochrone still a tautochrone as Vsauce says here, https://youtu.be/skvnj67YGmw?t=1090 ?
Keep up the good work!
christian 5 years, 9 months ago
I'm glad you enjoyed it. As far as I know, a brachistochrone is still a tautocrone for y2 > x2. Wikipedia has some proofs here:
Link | Replyhttps://en.wikipedia.org/wiki/Tautochrone_curve
and taking e.g. Abel's solution, I don't see why this condition should invalidate the proof.
Matt 5 years, 9 months ago
Thanks for the reply.
Link | ReplyIt occurred to me that when y2 >> x2 say, y2 = 1 and x2 = 0.01, the curve approaches a vertical drop. And a vertical drop is not a tautochrone: objects dropped from different heights reach the ground (end point) at different times. However, a not-quite-a-vertical-drop could still be described by the equation to a brachistochrone (one with a large cycloid radius), but presumably not fulfill the definition of a tautochrone. I'm curious to know the parameters whereby the brachistochrone ceases to be a tautochrone.
Trying to do this with Python, I hit a wall about here,
x2, y2 = 0.26173,1
...where slightly smaller values for x2 give,
scipy\optimize\zeros.py:173: RuntimeWarning: Tolerance of 4.76539196989e+27 reached
In my web-search for an answer, this excellent post appears to completely contradict my earlier statement, “So it is easy to break the Brachistochrone condition while keeping the isochrone [tautochrone] condition satisfied.” - https://physics.stackexchange.com/questions/16819/is-there-an-intuitive-reason-the-brachistochrone-and-the-tautochrone-are-the-sam
Do you think the brachistochrone is a general solution to the tautochrone or vice versa or are they perhaps mutually exclusive under certain circumstances?
christian 5 years, 9 months ago
I think the brachistochrone has to include the cusp of the cycloid (where the acceleration is greatest) and the tautochrone must end at its minimum. There's more on this here: https://math.stackexchange.com/a/2853385/126411
Link | ReplyAnathnath Ghosh 3 years, 5 months ago
Fantastic works.
Link | Replyhats off to you sir
Joaquín Herreros 2 years, 6 months ago
I have been playing with this curves for unitary radius (for aplying it to skateboard slides), and I have found a pretty good cartesian approximation (0<x<pi/2):
Link | Replyy(x) = cos^2(sqrt(5/2)+pi^2/2*{[log(2)]^(1/8)-[log(pi/x)]^(1/8)})
Is not too complicated and its fits very well. I hope it helps you too.
Best regards.
New Comment