I recently started getting interested in Julia after people insisting I look into it more. Unrelated to that, I have always enjoyed solving Classical Mechanics problems, so decided to try and combine the two with a very familiar problem: The Chaotic Double Pendulum. This is a nice, tractable problem that lets the ODE solving capabilities of Julia really shine.
Deriving Equations of Motion for Double Pendulum
We will derive the equations of motion using what is probably the most popular way: Lagrangian Mechanics. If you are interested in using Hamiltonian Mechanics, there is an excellent explanation here. The Lagrangian method provides a simpler implementation so we’ll use that.
First, we will define two identical pendulums of length \(\frac{l}{2}\), with the first being connected at one end to a central frictionless pivot, and the second being attached to the other end of the first pendulum, by a similar frictionless pivot. Having them be identical to each other allows us to make our equations of motions simpler to implement, however extending it to different masses of pendulum’s and length’s is readily done at the cost of making the mathematics a bit more tedious to derive. We will use the common short-hand \(\dot{\theta} = \frac{\partial \theta}{\partial t}\), where \(t\) is time, throughout.
First start with $$x_1 = \frac{l}{2} \sin\theta_1$$ $$y_1 = \frac{-l}{2}\cos\theta_1$$ $$x_2 = \frac{l}{2}\big(\sin\theta_1 + \sin\theta_2\big)$$ $$y_2 = \frac{-l}{2}\big(\cos\theta_1 + \cos\theta_2\big) $$
Lagrangian
Next let’s sketch out a generic Lagrangian $$ L = K – V $$ where \(K\) is the kinetic energy of the system and \(V\) is the potential energy. We can set $$ K = \frac{1}{2}m\big( v^2_1 + v^2_2 \big) + \frac{1}{2}I\big( \dot{\theta}^2_1 + \dot{\theta}^2_2 \big) $$ $$ V = mg\big(y_1 + y_2\big)$$ where \(v\) and \(\dot{\theta}\) represent the linear and rotational velocities in the system, and \(I\) the moment of inertia of the system. Modelling our pendulums as two rods of uniform density, we can use the known moment of inertia \(I = \frac{ml^2}{12} \) Which gives us our final Lagrangian
$$ L = \frac{1}{2}m\big( \dot{x}^2_1 + \dot{y}^2_1 + \dot{x}^2_2 + \dot{y}^2_2 \big) + \frac{ml^2}{24}\big( \dot{\theta}^2_1 + \dot{\theta}^2_2 \big) – mg\big(y_1 + y_2\big) $$
The derivatives are straightforward to calculate:
$$\dot{x}_1 = \frac{\dot{\theta}_1 l}{2} \cos \theta_1 $$
$$ \dot{y}_1 = \frac{\dot{\theta}_2 l}{2}\sin\theta_1 $$
$$ \dot{x}_2 = l\big( \dot{\theta}_1 \cos\theta_1 + \frac{\dot{\theta}_2}{2}\cos\theta_2 \big) $$
$$\dot{y}_2 = -2\dot{y}_1 + \frac{\dot{\theta}^2_2}{2}l \sin \theta_2 $$
After expanding out the Lagrangian and simplifying, which is a bit tedious but not difficult(omitted for brevity), we arrive at $$L = \frac{1}{6}ml^2 \bigg( \dot{\theta}^2_2 + 4\dot{\theta}_1^2 + 3\dot{\theta}_1\dot{\theta}_2 \cos\big( \theta_1 – \theta_2 \big)\bigg) + \frac{1}{2}mgl \bigg(3\cos\theta_1 + \cos\theta_2 \bigg) $$
Generalised Momenta
It is convenient to derive the “generalised momenta” of the system before we find the equations of motion for the system using the Euler-Lagrange equations, under the principle of the conservation of energy. These generalised momenta can be stated in our system (parameterised by degrees of freedom \(\theta_1\) and \(\theta_2\):
$$\rho_{\theta_1} = \frac{\partial L}{\partial \dot{\theta}_1} = \frac{1}{6}ml^2 \bigg( 8 \dot{\theta}_1 + 3\dot{\theta}_2 \cos\big(\theta_1 – \theta_2\big) \bigg) $$
$$\rho_{\theta_2} = \frac{\partial L}{\partial \dot{\theta}_2} = \frac{1}{6}ml^2 \bigg( 2 \dot{\theta}_2 + 3\dot{\theta}_1 \cos \big( \theta_1 – \theta_2 \big) \bigg) $$
Euler Lagrange Equations
The Euler Lagrange equations allow us to find the minimising function of a differentiable functional, which in this case is our Lagrangian, being the functional of the kinetic and potential energies (and their derivatives) of the system parameterised by its fundamental degrees of freedom \(\theta_1\) and \(\theta_2\):
$$ \frac{d}{dt} \Bigg( \frac{\partial L}{\partial \dot{\theta}_i} \Bigg) – \frac{\partial L}{\partial \theta_i} = 0 $$
The term on the right, \(\frac{\partial L}{\partial \theta_i} \), we can derive in a straightforward manner:
$$ \frac{\partial L}{\partial \theta_1} = \frac{-ml^2}{2} \dot{\theta}_1 \dot{\theta}_2 \sin \big( \theta_1 – \theta_2 \big) – \frac{3}{2} mgl \sin\theta_1$$
$$ \frac{\partial L}{\partial \theta_2} = \frac{=ml^2}{2} \dot{\theta}_1 \dot{\theta}_2 \sin \big( \theta_1 – \theta_2 \big) – \frac{1}{2} mgl \sin\theta_2$$
Handling the term \(\frac{d}{dt} \Big( \frac{\partial L}{\partial \dot{\theta}_i} \Big)\) is more complicated. Previously we did partial derivatives where we only differentiate the parameter being considered, and hold all other symbols as constants. Here, however, we must consider all symbols as \(f(t)\) and thus the differentiation becomes more tedious. Let’s demonstrate this by taking the first parameter, and substituting \(\rho_{\theta_1}\) into \(\frac{d}{dt} \Big( \frac{\partial L}{\partial \dot{\theta}_i} \Big) = \frac{d}{dt}\Big( \rho_{\theta_1} \Big ) \):
$$ \frac{d}{dt} \Big( \rho_{\theta_1} \Big) = \frac{1}{6} \big( 8 \ddot{\theta}_1 + 3\ddot{\theta}_2\cos\big( \theta_1 – \theta_2 \big) – 3 \dot{\theta}_1 \dot{\theta}_2\sin\big( \theta_1 – \theta_2 \big) + 3\dot{\theta}^2_2 \sin\big( \theta_1 – \theta_2 \big)$$
Repeating for the second parameter, \(\theta_2\) yields:
$$ \frac{d}{dt} \Big( \rho_{\theta_2} \Big) = \frac{1}{6} \big( 2 \ddot{\theta}_2 + 3\ddot{\theta}_1\cos\big( \theta_1 – \theta_2 \big) – 3 \dot{\theta}_1 \dot{\theta}_2\sin\big( \theta_1 – \theta_2 \big) – 3\dot{\theta}^2_1 \sin\big( \theta_1 – \theta_2 \big)$$
Now we can combine our terms together into the final Euler Lagrange equation above, to end up with 2 second order Ordinary Differential Equations. After some simplification and cancelling these are:
$$\ddot{\theta}_1 + \frac{3}{8}\ddot{\theta}_2 \cos\big( \theta_1 – \theta_2 \big) = \frac{-3}{8}\dot{\theta}^2_2 \sin \big( \theta_1 – \theta_2 \big) – \frac{9mg}{8l} \sin \theta_1 $$
$$ \ddot{\theta}_2 + \frac{3}{2} \ddot{\theta}_1 \cos \big( \theta_1 – \theta_2 \big) = \frac{3}{2} \dot{\theta}_1^2 \sin \big( \theta_1 – \theta_2 \big) – \frac{3mg}{2l} \sin\theta_2 $$
At this point we could attempt to directly solve these coupled ODE’s numerically, but they will turn out to be unstable, and this is a general problem with solving higher order differential equations on machines with limited numerical precision.
First Order Reduction
As stated, you can try and solve the equations above numerically but I found that they showed significant instability over longer time epochs. It is a generally known property that \(N\) coupled second order ODE’s can be reduced to \(2N\) coupled first order ODE’s. This will help us resolve the ill conditioning of the coupled second order equations, if we can find a way to facilitate this order reduction, which can be often non-trivial.
Lets denote the first and second equation as \(f_1\) and \(f_2\). Secondly, you will notice that the first of our coupled ODE’s above does not depend on \(\dot{\theta}_1\) and the second equation does not depend on \(\dot{\theta}_2\). We can now write the couple second order ODE’s as
$$ A \begin{bmatrix} \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{bmatrix}$$
$$ A = \begin{bmatrix} 1 & \alpha_1 \\ \alpha_2 & 1 \end{bmatrix} $$
We can use the fortunate simplicity of the inverse of a 2\(\times\)2 matrix to write the inverse of \(A\):
$$ A^{-1} = \frac{1}{det(A)} \begin{bmatrix} 1 & -\alpha_1 \\ -\alpha_2 & 1 \end{bmatrix} $$
$$ det(A) = 1 – \alpha_1 \alpha_2 $$
We can work out \(\alpha_1\) and \(\alpha_2\) by inspection of our above equations:
$$ \alpha_1 = \frac{3}{8}\cos\big( \theta_1 – \theta_2 \big) $$
$$ \alpha_2 = \frac{3}{2}\cos\big( \theta_1 – \theta_2 \big) $$
And we can observe the very useful property:
$$ \frac{1}{1 – \alpha_1\alpha_2} > 0 $$
meaning that \(A\) is indeed invertible as the determinant must always be more than 0.
Putting it all together we end up with:
$$ \begin{bmatrix} \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{bmatrix} = A^{-1} \begin{bmatrix} f_1 \\ f_2 \end{bmatrix} $$
$$ \begin{bmatrix} \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{bmatrix} = \frac{1}{1 – \alpha_1\alpha_2} \begin{bmatrix} 1 & -\alpha_1 \\ -\alpha_2 & 1 \end{bmatrix} \begin{bmatrix} f_1 \\ f_2 \end{bmatrix} $$
$$ \begin{bmatrix} \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{bmatrix} = \frac{1}{1 – \alpha_1\alpha_2} \begin{bmatrix} f_1 -\alpha_1 f_2 \\ -\alpha_2 f_1 + f_2 \end{bmatrix} $$
From above we can describe \(f_1\) and \(f_2\):
$$ f_1 = \frac{-3}{8} \dot{\theta}^2_2 \sin \big( \theta_1 – \theta_2 \big) – \frac{9g}{8} \sin \theta_1 $$
$$ f_2 = \frac{3}{2} \dot{\theta}^2_1 \sin \big( \theta_1 – \theta_2 \big) – \frac{3g}{2} \sin \theta_2 $$
where I have set \(m=l=1\) for simplicity.
Now we can form \(2N\) first order ODE’s which allow us to relate \(\theta_i\), \(\dot{\theta}_i \), and \(\ddot{\theta}_i\):
$$ \frac{d}{dt} \begin{bmatrix} \theta_1 \\ \theta_2 \\ \dot{\theta}_1 \\ \dot{\theta}_2 \end{bmatrix} = \begin{bmatrix} \dot{\theta}_1 \\ \dot{\theta}_2 \\ \chi_1 \\ \chi_2 \end{bmatrix} $$
$$ \chi_1 = \frac{f_1 – \alpha_1 f_2}{1 – \alpha_1 \alpha_2 } $$
$$ \chi_2 = \frac{-\alpha_2 f_1 + f_2}{1 – \alpha_1 \alpha_2 } $$
Implementation in Julia
We will show implementations of both the second order ODE solver and the simplified first order ODE solver. These are both implemented in Julia using the DifferentialEquations package. All plotting is done with the Plots package and is the same for both solvers. As can be seen the implementation of the final mathematical derivation is far quicker and simpler than the mathematical derivation itself! Here I use the Rosenbrock23 solver which is resilient to stiff ODE’s with AutoTsit which will adjust the solver hyper parameters dependent on observed evolution parameters. For the second order equation solver, a standard Runge-Kutta 4 solver would also suffice with added benefits of improved runtime.
Second Order Equation Solver
using DifferentialEquations
using Plots
function doublePendulum!(ddu,du,u,p,t)
ddu[1] = ( 3.0 / 8.0 ) * ( du[2] * du[2] * sin(u[1] - u[2]) - ddu[2] * cos( u[1] - u[2] ) ) - (9.0 * 9.81 / 8.0 ) * sin( u[1] )
ddu[2] = 1.5 * ( du[1] * du[1] * sin(u[1] - u[2]) - ddu[1] * cos( u[1] - u[2] ) ) - 1.5 * 9.81 * sin( u[2] )
end
g = 9.81
## Second Order ODE Solver
du0 = [0.0, 0.0]
u0 = [0.2, 0.2]
prob2 = SecondOrderODEProblem(doublePendulum!,du0, u0, tspan)
sol = solve(prob2,AutoTsit5(Rosenbrock23()), reltol=1e-14, abstol=1e-14,dt=0.001);
Sadly as the below image demonstrates, this system is numerically unstable and results in ever increasing amounts of total energy in the system, which is non physical behaviour.
Simplified First Order Equation Solver
using DifferentialEquations
using Plots
g = 9.81
function doublePendulum2!(du,u,p,t)
du[1] = u[3]
du[2] = u[4]
cosO1O2 = cos( u[1] - u[2] )
sinO1O2 = sin( u[1] - u[2] )
alpha_1 = ( 3.0 / 8.0 ) * cosO1O2
alpha_2 = ( 3.0 / 2.0 ) * cosO1O2
det_recip = 1 - alpha_1 * alpha_2
f1 = - ( 3.0 / 8.0 ) * u[4] * u[4] * sinO1O2 - ( 9.0 * g / 8.0 ) * sin( u[1] )
f2 = ( 3.0 / 2.0 ) * u[3] * u[3] * sinO1O2 - ( 3.0 * g / 2.0 ) * sin( u[2] )
# f_1 - \alpha f_2
du[3] = ( f1 - alpha_1 * f2 ) / det_recip
du[4] = ( f2 - alpha_2 * f1 ) / det_recip
end
# Coupled First Order ODE Solver
u0 = [π/2,π/2,0.,0.]
prob3 = ODEProblem(doublePendulum2!,u0,tspan)
sol=solve(prob3, abs_tol=1e-10, rel_tol=1e-10, dt=0.01,adaptive=false);
Display Code
The display code uses the plots library and some basic coordinate conversions and vector arithmetic to display end results in terms of the base parameter evolution of the system, and as a GIF. To visualise the linked together pendulum’s we can hack together a visualisation using a quiver plot.
using Plots
function radiansToXY(l,r)
x = sin.(r) * l
y = cos.(r) * l
return -x,-y
end
p = plot(sol.t,sol[1,:])
p = plot!(sol.t,sol[2,:])
p = plot!(sol.t,sol[3,:])
p = plot!(sol.t,sol[4,:])
savefig(p,"plotAngles.pdf")
pend1 = radiansToXY.(1,sol[1,:])
pend2 = radiansToXY.(1,sol[2,:])
pend1 = [tup[k] for k in 1:2, tup in pend1]
pend2 = [tup[k] for k in 1:2, tup in pend2]
xPend1 = pend1[1,:]
yPend1 = pend1[2,:]
xPend2 = pend2[1,:]
yPend2 = pend2[2,:]
xPend1Plus2 = xPend1+xPend2
yPend1Plus2 = yPend1+yPend2
p = scatter(xPend1,yPend1,markersize = 2)
p = scatter!(xPend1Plus2,yPend1Plus2,markersize = 2)
ylims!((-3.,3.))
xlims!((-3.,3.))
plot!(size=(400,400))
savefig(p,"plotCartesian.pdf")
# Make GIF of system evolution
a = Animation()
for i in 1:size(sol.t)[1]
p = scatter([xPend1[i]],[yPend1[i]],markersize=3)
vx1 = [-xPend1[i]]
vy1 = [-yPend1[i]]
p = quiver!([xPend1[i]],[yPend1[i]],quiver=(vx1,vy1))
p = scatter!([xPend1Plus2[i]],[yPend1Plus2[i]],markersize=3)
vx2 = [-xPend1Plus2[i] + xPend1[i]]
vy2 = [-yPend1Plus2[i] + yPend1[i]]
p = quiver!([xPend1Plus2[i]],[yPend1Plus2[i]],quiver=(vx2,vy2))
ylims!((-3.,3.))
xlims!((-3.,3.))
plot!(size=(400,400))
frame(a,p)
end
gif(a)
Results
Lastly let’s look at some examples of the evolution in time of a double pendulum system starting from no initial velocity and some displacement from the origin of each pendulum.
In the second and third examples we get a great view of the chaotic behaviour of the system in the classical unconstrained double pendulum. Julia’s ability to solve ODE’s with minimal fuss has completely won me over, to the point that I am thinking of going further and writing a Partial Differential Equation solver for some simplified system such as an oscillating dipole with Maxwell’s Equations. It is even better than Python in terms of the brevity of its syntax, and is far faster due to it running natively through C code.