A double pendulum consists of 2 pendula, one of which hangs off of the second. Below, the angles $\theta_1$ and $\theta_2$ give the position of the red ball ($m_1$) and green ball ($m_2$) respectively. The lengths of the (massless) rod holding the balls to the pivot are $L_1$ and $L_2$ respectively as well.
The cartesian coordinates $x_1, y_1, x_2, y_2$ are given by:
| $x_1 = L_1\sin\theta_1$ | $y_1 = -L_1\cos\theta_1$ | |
| $x_2 = x_1 + L_2\sin\theta_2$ | $y_2 = y_1 - L_2\cos\theta_2$ |
Their derivatives (needed later) are:
| $\dot x_1 = L_1\cos\theta_1$ | $\dot y_1 = L_1\sin\theta_1$ | |
| $\dot x_2 = \dot x_1 + L_2\cos\theta_2$ | $\dot y_2 = \dot y_1 + L_2\sin\theta_2$ |
The total kinetic energy is given by:
$KE = \frac{1}{2} m_1 (\dot x_1^2 + \dot y_1^2) + \frac{1}{2} m_2 (\dot x_2^2 + \dot y_2^2)$
which after some algebra reduces to: $$KE = \frac{1}{2} M L_1^2 \dot\theta_1^2 + \frac{1}{2} m_2 L_2^2 \dot\theta_2^2 + m_2 L_1 L_2 \dot\theta_1\dot\theta_2\cos\Delta\theta\label{KE}$$ where $M\equiv m_1 + m_2$ and $\Delta\theta\equiv \theta_1 - \theta_2$.
The potential energy due to gravity is the usual:
$PE = m_1 g y_1 + m_2 g y_2$
which reduces to
$$PE = -MgL_1\cos\theta_1 - m_2gL_2\cos\theta_2\label{PE}$$
The full Lagrangian is therefore: $$\mathcal{L} = \frac{1}{2} M L_1^2 \dot\theta_1^2 + \frac{1}{2} m_2 L_2^2 \dot\theta_2^2 + m_2 L_1 L_2 \dot\theta_1\dot\theta_2\cos\Delta\theta + MgL_1\cos\theta_1 + m_2gL_2\cos\theta_2\label{Lagrangian}$$ After applying the Euler-Lagrange equations to $(\dot\theta_1,\theta_1)$ and $(\dot\theta_2,\theta_2)$ and simplying a bit we get the following two equations:
$$ML_1\ddot\theta_1 + m_2L_2\ddot\theta_2\cos\Delta\theta + m_2L_2\dot\theta_2^2\sin\Delta\theta + Mg\sin\theta_1 = 0\label{eqth1}$$ $$L_1\ddot\theta_1\cos\Delta\theta + L_2\ddot\theta_2 - L_1\dot\theta_1^2\sin\Delta\theta + g\sin\theta_2 = 0\label{eqth2}$$
These are of course very non-linear coupled systems, not solvable analytically.
It is useful first to take a few limits. For instance, setting $m_2 = 0$, $L_2 = 0$, and $\theta_1 = \theta_2$, both equations become the simple pendulum equation $\ddot\theta_1 + \omega^2\sin\theta_1 = 0$ where $\omega^2 \equiv g/L$.
To solve these equations numerically in a simulation, we first have to rearrange into two equations, each of which have only 1 second derivative in the time. After a lot of algebra you should get the following two equations: $$\ddot\theta_1 = [-\sin\Delta\theta(m_2L_1\dot\theta_1^2\cos\Delta\theta + m_2L_2\dot\theta_2^2) - g(M\sin\theta_1 - m_2\sin\theta_2\cos\Delta\theta)]/L_1\alpha\label{eq1}$$ $$\ddot\theta_2 = [+\sin\Delta\theta(ML_1\dot\theta_1^2 + m_2L_2\dot\theta_2^2\cos\Delta\theta) + g(M\sin\theta_1\cos\Delta\theta - M\sin\theta_2)]/L_2\alpha\label{eq2}$$ and where $\alpha \equiv m_1 + m_2\sin^2\Delta\theta$.
In the simulation below, we use 3 common methods for the numerical integration: Euler's method; the modified Euler-Cromer; and Runge-Kutta (order 2, RK2). All 3 start with the 2 basic equations, ($\ref{eq1}$) and ($\ref{eq2}$), written in the following way:
$\ddot\theta_1 = f_1(\theta_1,\theta_2,\dot\theta_1,\dot\theta_2)$
$\ddot\theta_2 = f_2(\theta_1,\theta_2,\dot\theta_1,\dot\theta_2)$
where $f_1$ is given by equation ($\ref{eq1}$) and $f_2$ by equation ($\ref{eq2}$).
What we want is a way to use our knowledge of $\theta_1(t)$ to get $\theta_1(t+\dt)$ (and the same for $\theta_2, \dot\theta_1, and \dot\theta_2$). To do this, we first consider $\theta_1$ and make the substitution:
$\alpha\equiv\theta_1$
$\beta\equiv\dot\theta_1$
This allows us to write equation ($\ref{eq1}$), a second order differential equation, in terms of 2 first order differential equations:
$\dot\beta = f_1(\theta_1,\theta_2,\dot\theta_1,\dot\theta_2)$
$\dot\alpha = f_1'(\beta)=\beta$
As discussed in the section on numerical integration, the Euler technique uses the above 2 equations to numerically find $\alpha(t+\dt)$ and $\beta(t+\dt)$ based on knowledge of $\alpha(t)$, $\beta(t)$, using the approximation (based on calculus definitions of derivatives): $$\alpha(t+\dt)=\alpha(t)+\dt\cdot f_1'(\beta)\label{eqa1}$$ $$\beta(t+\dt)=\beta(t)+\dt\cdot f_1\label{eqa2}$$ This is basically just using the definition of derivatives before taking the limit:
$\frac{d\alpha}{dt}=f_1'=\lim_{\dt\to 0}\frac{\alpha(t+\dt)-\alpha(t)}{\dt}$
and solving for $\alpha(t+\dt)$ to get equation ($\ref{eqa1}$). And similarly for $\beta(t+\dt)$.
In words, what we are doing here is using the knowledge of $\alpha(t)$, the value at the beginning of the interval, and "swimming" to $\alpha(t+\dt)$ using the value of the derivative at the beginning of the interval. So to be explicit: $$\alpha(t+\dt)=\alpha(t)+\dt\cdot f_1'(\beta(t))$$ $$\beta(t+\dt)=\beta(t)+\dt\cdot f_1(\alpha(t),\dot\alpha(t),\beta(t),\dot\beta(t))$$ The Euler-Cromer technique is just a variation where we use the slope (derivative) not at the beginning, but at the end of the interval: $$\alpha(t+\dt)=\alpha(t)+\dt\cdot f_1'(\beta(t+\dt))$$ $$\beta(t+\dt)=\beta(t)+\dt\cdot f_1(\alpha(t+\dt),\dot\alpha(t+\dt),\beta(t+\dt),\dot\beta(t+\dt))$$ And finally, the Runge-Kutta (2nd order, RK2) uses the value of the slope in the middle of the interval: $$\alpha(t+\dt)=\alpha(t)+\dt\cdot f_1'(\beta_h)$$ $$\beta(t+\dt)=\beta(t)+\dt\cdot f_1(\alpha_h,\dot\alpha_h,\beta_h,\dot\beta_h)$$ where $$\alpha_h\equiv\alpha(t)+\half\dt\cdot \dot\alpha(t)$$ $$\dot\alpha_h\equiv\dot\alpha(t)+\half\dt\cdot f_1'(\beta(t))$$ and etc for $\beta$ and $\dot\beta$. The RK2 really amounts to starting at the beginning of the interval and using Euler to calculate the slope at the midpoint of the interval, $f_1(\alpha_h,\dot\alpha_h,\beta_h,\dot\beta_h)$, to find to get us to the end of the interval. For more, surf to the section on numerical_integration.
In the simulation below, you can choose which 3 numerical integration methods you want to use, change the time interval $\dt$, the value of the two masses $m_1$ and $m_2$, and the value of the two pendulum lengths $L_1$ and $L_2$ by clicking on the arrows (hold the arrows down and the masses and pendlums and it will change continuously). You can even change the masses and lengths while the thing is in motion, just for fun! Note that each of the 3 integration methods will have a different color, so you can run them all and see them on the two plots on the left (which are described below).
To get a feel for how good each method is as an approximation, the left graph labeled "Energy (arbitrary scale..." plots the total energy E = KE + PE as a function of time. Any deviation from a constant value (dashed line) shows how good (or bad) an approximation the method is. As you can see quite readily, Euler's approximation shows the total energy diverging, which makes it a pretty poor approximation for such a non-linear system. The graph below this, labeled "Phase space", plots the time derivative of the angle $\phi_1$ ($\dot\phi_1$) vs the angle $\phi$. The title "phase space" comes out of the idea that for a moving particle, what's important is not just position but also velocity (or momentum), and plotting those two shows the functional relationship. We know that the Euler approximation violates energy conservation quite readily, and this is clear by looking at either plots: the path will eventually be way off scale and you will quickly see this as the signal gets beyond the plot boundaries. The Runge-Kutta approximation, however, shows that this is a much better approximation for this highly non-linear system. What's interesting is to shrink the length $L_2$ so that it's close to zero, run the simulation, and look at the phase space plot. This should look pretty much just like a regular single pendulum, and be highly repetitive. The phase space plot should show this - the curve repeats itself. As you let $L_2$ get bigger, you will see more and more chaotic behavior. But if you click on the "Single" radio button and hit start, the simulation will show a single pendulum only. You can see that the Euler-Cromer is actually pretty good for a single pendulum, and RK2 is more or less perfect over a large amount of time.
| $\dt$: 0.1 $L_1$: 150 $L_2$: 150 $m_1$: 15 $m_2$: 15 | |
| Euler Euler-Cromer Runge-Kutta | |
| Single Double | |
| Energy (arbitrary scale, dashed line is $E_{tot}$)
Phase space $\dot\phi_1$ vs $\phi_1$ | |
The double pendulum has some amazing properties, as you can see. The pendulum itself goes through a path that is ultra sensitive to initial conditions, and the phase space plot has very irregular orbits. All of this leads to the field of chaos theory. For more, try here.
Back to TOP
Copyright © 2019 by Drew Baden
All rights reserved. No part of this publication may be reproduced, distributed, or transmitted in any form or by any means, including photocopying, recording, or other electronic or mechanical methods, without the prior written permission of the publisher, except in the case of brief quotations embodied in critical reviews and certain other noncommercial uses permitted by copyright law.