using PyPlot, Interact, QuadGK, ForwardDiff
This IJulia notebook contains some computational examples and homework exercises to accompany a January 2018 lecture by Prof. Steven G. Johnson on the calculus of variations, motivated by the famous brachistochrone problem. This lecture is part of the MIT IAP course 18.095 (Mathematics Lecture Series).
For the most part, this notebook does not include any derivations, which were presented on the blackboard. Hopefully, a transcription of those notes can be posted soon.
This is the unlikely name of a famous problem in mathematics, with a long history going back to Bernoulli in 1696:
Suppose a point mass is sliding down a curve $y(x)$, without friction, from point $A$ to point $B$. What curve $y(x)$ yields the smallest time?
Before we analyze this problem mathematically, let's have the computer crunch some numbers for us to see what happens for a few different curve shapes:
const g = 9.80655 # the acceleration of gravity, in m/s²
9.80655
Given a curve $y(x)$ from $x=a$ to $x=b$, a little bit of physics and calculus tells us that the transit time is the functional
$$ T[y] = \sqrt{\frac{1}{2g}} \int_a^b \sqrt{\frac{1 + \left(\frac{dy}{dx}\right)^2}{y(a) - y(x)}} dx $$where $g$ is the gravitational acceleration.
This is a tricky integral to do for an arbitrary curve $y(x)$, but the computer can calculate it (to any desired accuracy) with no problem:
# numerical integral of the transit time for y(x) from x=a to x=b, explained elsewhere
T(y, a,b) = quadgk(x -> sqrt((1 + ForwardDiff.derivative(y, x)^2)/(2g*(y(a)-y(x)))), a,b, reltol=1e-4, abstol=1e-3)[1]
T (generic function with 1 method)
x = linspace(0,1,200)
fig = figure()
@manipulate for p = slider(0.05:0.05:4, value=1.0)
withfig(fig) do
y(x) = (1-x)*exp(-(-2+2^p)*x)
Ty = T(y,0,1)
plot(x, y.(x), "r-")
plot(0,1, "ko")
plot(1,0, "ko")
text(0.05,1,"A")
text(0.98,0.05,"B")
xlabel(L"$x$ (meters)")
ylabel(L"$y$ (meters)")
title("total time = $Ty seconds")
end
end
Obviously, it is better if the curve is steeper in the beginning, to get the mass moving more quickly. However, if it is too steep, then the time is longer simply because the distance is larger.
There is an optimum. Some curve that is steep but not too steep, just right to minimize the time. Actually finding this curve is quite tricky, and leads to the beautiful topic of the calculus of variations: a general approach to solving problems of minimizing functionals.
The exact solution (the minimum-time curve) turns out to be a shape called a cycloid. This is most simply written as a parametric curve in terms of $x(\theta)$ and $y(\theta)$ for $\theta \in [0,\phi]$:
$$ x(\theta) = a + k (\theta - \sin\theta) \\ y(\theta) = y(a) - k (1 - \cos\theta) $$where $k (1-\cos\phi) = y(a) - y(b)$ and $\phi$ is determined by solving the transcendental equation:
$$ \frac{\phi - \sin\phi}{1 - \cos\phi} = \frac{b-a}{y(a)-y(b)} $$On the computer, we can solve for $\phi$ by Newton's method:
# find root of r(ϕ) = ϕ - sinϕ - α (1 - cosϕ), where α = (b-a)/(y(a)-y(b)) = 1 in our case
# by Newton iterations ϕ ⟵ ϕ - r(ϕ)/r′(ϕ)
ϕ = 2.5 # initial guess
while true
Δϕ = (ϕ - sin(ϕ) + cos(ϕ) - 1) / (1 - cos(ϕ) - sin(ϕ))
ϕ -= Δϕ
abs(Δϕ) < 1e-12 && break
end
θ = linspace(0,pi,200)
plot(θ, (θ .- sin.(θ))./(1 .- cos.(θ)), "b-")
plot(θ, ones(θ), "k--")
plot(ϕ, (ϕ - sin(ϕ))/(1 - cos(ϕ)), "ro")
xlabel(L"\theta")
ylabel(L"\frac{\theta - \sin\theta}{1 - \cos\theta}")
title("solving the transcendental equation for the endpoint \$\\phi\$")
println("ϕ = $ϕ")
ϕ = 2.4120111439135252
k = 1 / (1 - cos(ϕ))
0.5729170375317504
Now that we know $\phi$ and $k$, let's plot the curve:
θ = linspace(0,ϕ,200)
plot(k*(θ .- sin.(θ)), 1 .- k*(1. - cos.(θ)), "b-")
plot(0,1, "ko")
plot(1,0, "ko")
text(0.05,1,"A")
text(0.98,0.05,"B")
title("Brachistochrone solution: Cycloid curve")
ylabel(L"$y$ (meters)");
In terms of this parameterization, the minimum time is $$T_0 = \sqrt{\frac{1}{2g}} \int_0^{\phi} \sqrt{\frac{\left( \frac{dx}{d\theta} \right)^2 + \left( \frac{dy}{d\theta} \right)^2}{y(a) - y}} d\theta = \sqrt{\frac{k}{2g}} \int_0^{\phi} \sqrt{\frac{(1 - \cos\theta)^2 + \sin^2 \theta}{1 - \cos\theta}} d\theta $$. We can compute this integral numerically:
T₀ = sqrt(k/(2g)) * quadgk(θ -> sqrt(((1-cos(θ))^2+sin(θ)^2) / (1 - cos(θ))), 0, ϕ)[1]
0.5829979871060635
This minimum time $T_0 \approx 0.583$ is very similar to what we observed above in our little numerical experiment for different curves. Let's try it again, but this time show $T[y]/T_0$ to see how the different curves compare to the minimum:
x = linspace(0,1,200)
fig = figure()
@manipulate for p = slider(0.05:0.05:4, value=1.0)
withfig(fig) do
y(x) = (1-x)*exp(-(-2+2^p)*x)
Ty = T(y,0,1) / T₀
plot(x, y.(x), "r-")
plot(k*(θ .- sin.(θ)), 1 .- k*(1. - cos.(θ)), "b-")
legend([L"y(x)", "brachistochrone"])
plot(0,1, "ko")
plot(1,0, "ko")
text(0.05,1,"A")
text(0.98,0.05,"B")
xlabel(L"$x$ (meters)")
ylabel(L"$y$ (meters)")
title("\$T[y] / T\_0\$ = $Ty")
end
end
At best, this "guessed" $y(x)$ is takes about 1% longer than the minimum-time (brachistochrone) curve.
Although the calculus of variations can sometimes lead to very complicated ordinary differential equations (ODEs), the following problems lead to trivial ODEs that you should be able to solve even if you haven't taken 18.03:
The functional $L[y] = \int_a^b \sqrt {1 + (y')^2} dx$ gives the length of the curve $y(x)$ from $x=a$ to $x=b$. Assuming the endpoints are fixed $y(a)=0$ and $y(b)=y_B$, show that the solution of the Euler–Lagrange equations (via the Beltrami identity) give a straight line. (i.e. the shortest curve between two points is a straight line.)
Consider the functional $F[y] = \frac{\int_0^\pi (y')^2 dx} {\int_0^\pi y^2 dx}$. Suppose that we want to minimize this over all real-valued functions $y(x)$ with $y(0) = y(\pi) = 0$. This is not quite in the same form as the functionals covered in class, so you can't just blindly quote the Euler–Lagrange equations!