Tyler Jones
University of Wisconsin-Madison
DoublePendulumPerturbationAnimation.m
DoublePendulumFractalEvolution.m
This is the two dimensional restricted (yet classical) double pendulum problem. The masses, sometimes modeled as evenly distributed along its length, will be simplified as two point masses \(m_1=m_2=m \) connected by two light limbs/rods of lengths \(l_1=l_2=l \) with angles \(\theta_1 \) and \(\theta_2 \) from the vertical. If we wish to model the system similar to the inspired paper below, the center of mass will reside at the center of each limb/rod which will result in a moment of inertia: \(I=\frac{1}{12}ml^2 \) where \(l \) is the rigid length of each limb/rod; however, I will treat each pendulum as a point mass residing in the center of each ball. The goal of this project is to apply my knowledge of Lagrangian and Hamiltonian formalisms (Physics 311) as well as applied dynamical systems (Math 415/519) to predict when the second pendulum will 'flip'. Unlike the simple pendulum one might analyze in an introductory physics class, the double pendulum is chaotic. For our purpose, the definition of chaos can be defined as follows: A system that is sensitive and highly dependent on initial conditions, hence small perturbations in these initial conditions will result in vastly different behavior in finite time. After exploring the beauty of chaos, I will employ numerical methods via MATLAB in an attempt to make order from chaos. The goal is to be able to predict initial conditions that DO NOT allow the second pendulum to 'flip', but if it does, how long will it take? To gain a better understanding of my inspiration, follow the link: The Double Pendulum Fractal by Jeremy S. Heyl
In order to define the system, I will use Lagrangian and Hamiltonian formalisms as previously mentioned. To start, I will derive the Lagrangian in standard
Cartesian coordinates as follows:
$$
\begin{cases}
x_1 = l\sin\theta_1 \\
y_1 = -l\cos\theta_1 \\
x_2 = x_1 + l\sin\theta_2 = l\sin\theta_1 + l\sin\theta_2 \\
y_2 = y_1 - l\cos\theta_2 = -l\cos\theta_1 - l\cos\theta_2 \\
\end{cases}
$$
$$
\begin{cases}
\dot{x_1} = l\dot{\theta_1}\cos\theta_1 \\
\dot{y_1} = l\dot{\theta_1}\sin\theta_1 \\
\dot{x_2} = \dot{x_1} + l\dot{\theta_2}\cos\theta_2 = l\dot{\theta_1}\cos\theta_1 + l\dot{\theta_2}\cos\theta_2 \\
\dot{y_2} = \dot{y_1} + l\dot{\theta_2}\sin\theta_2 = l\dot{\theta_1}\sin\theta_1 + l\dot{\theta_2}\sin\theta_2 \\
\end{cases}
$$
$$
T := \frac{1}{2}m|v|^2= \frac{1}{2}m(\dot{x_1}^2 + \dot{y_1}^2 + \dot{x_2}^2 + \dot{y_2}^2)
=\frac{1}{2}ml^2(2\dot{\theta_1}^2 + \dot{\theta_2}^2 + 2\dot{\theta_1}\dot{\theta_2}\cos(\theta_1-\theta_2))
$$
$$
V := -\int F_g(y)dy = mgy_1 + mgy_2 = -mgl(2\cos\theta_1 + \cos\theta_2)
$$
$$ \boxed{L := T-V=\frac{1}{2}ml^2(2\dot{\theta_1}^2 + \dot{\theta_2}^2 + 2\dot{\theta_1}\dot{\theta_2}\cos(\theta_1-\theta_2)) + mgl(2\cos\theta_1 + \cos\theta_2)} $$
Now that we have derived the Lagrangian of our system, we can use the Euler-Lagrange equations in order to obtain the equations
of motion (EOM). Note that (\(q,\dot{q} \)) are the generalized coordinates.
$$
\frac{d}{dt}\Big(\frac{\partial L}{\partial \dot{q}} \Big) - \frac{\partial L}{\partial q}=0
$$
For \(\theta_1\):
$$
\frac{d}{dt}\Big(\frac{\partial L}{\partial \dot{\theta_1}} \Big) - \frac{\partial L}{\partial \theta_1}=0
$$
$$
\frac{\partial L}{\partial \dot{\theta_1}} = ml^2(2\dot{\theta_1} + \dot{\theta_2}\cos(\theta_1 - \theta_2))
$$
$$
\frac{\partial L}{\partial \theta_1} = -mgl(2\sin\theta_1) - ml^2(\dot{\theta_1}\dot{\theta_2}\sin(\theta_1 - \theta_2))
$$
$$
\frac{d}{dt}\Big(\frac{\partial L}{\partial \dot{\theta_1}} \Big) = ml^2(2\ddot{\theta_1} + \ddot{\theta_2}\cos(\theta_1 - \theta_2) - \dot{\theta_2}(\dot{\theta_1} + \dot{\theta_2})\sin(\theta_1 - \theta_2))
$$
$$
ml^2(2\ddot{\theta_1} + \ddot{\theta_2}\cos(\theta_1 - \theta_2) - \dot{\theta_2}(\dot{\theta_1} + \dot{\theta_2})\sin(\theta_1 - \theta_2)) + mgl(2\sin\theta_1) + ml^2(\dot{\theta_1}\dot{\theta_2}\sin(\theta_1 - \theta_2)) = 0
$$
$$
2\ddot{\theta_1} + \ddot{\theta_2}\cos(\theta_1 - \theta_2) + \frac{g}{l}(2\sin\theta_1) - \dot{\theta_2}(\dot{\theta_1} + \dot{\theta_2})\sin(\theta_1 - \theta_2) = 0
$$
For \(\theta_2\):
$$
\frac{d}{dt}\Big(\frac{\partial L}{\partial \dot{\theta_2}} \Big) - \frac{\partial L}{\partial \theta_2}=0
$$
$$
\frac{\partial L}{\partial \dot{\theta_2}} = ml^2(\dot{\theta_2} + \dot{\theta_1}\cos(\theta_1 - \theta_2))
$$
$$
\frac{\partial L}{\partial \theta_2} = -mgl\sin\theta_2 - ml^2(\dot{\theta_1}\dot{\theta_2}\sin(\theta_1 - \theta_2))
$$
$$
\frac{d}{dt}\Big(\frac{\partial L}{\partial \dot{\theta_2}} \Big) = ml^2(\ddot{\theta_2} + \ddot{\theta_1}\cos(\theta_1 - \theta_2) - \dot{\theta_1}(\dot{\theta_1} + \dot{\theta_2})\sin(\theta_1 - \theta_2))
$$
$$
ml^2(\ddot{\theta_2} + \ddot{\theta_1}\cos(\theta_1 - \theta_2) - \dot{\theta_1}(\dot{\theta_1} + \dot{\theta_2})\sin(\theta_1 - \theta_2)) + mgl\sin\theta_2 + ml^2(\dot{\theta_1}\dot{\theta_2}\sin(\theta_1 - \theta_2)) = 0
$$
$$
\ddot{\theta_2} + \ddot{\theta_1}\cos(\theta_1 - \theta_2) + \frac{g}{l}\sin\theta_2 - \dot{\theta_1}(\dot{\theta_1} + \dot{\theta_2})\sin(\theta_1 - \theta_2) = 0
$$
To solve for \( \ddot{\theta_1} \) and \( \ddot{\theta_2} \) in terms of \( \theta_1 \), \( \dot{\theta_1} \), \( \theta_2 \), and \( \dot{\theta_2} \),
we use the final expressions above. This can be done by hand or by using an equation solver.
The equations of motion for \( \ddot{\theta_1} \) and \( \ddot{\theta_2} \) are then given by:
$$
\begin{cases}
\begin{aligned}
\ddot{\theta_1} &= \frac{-2g\sin\theta_1 - g\sin(\theta_1 - 2\theta_2) - 2\sin(\theta_1 - \theta_2)(\dot{\theta_2}^2 l + \dot{\theta_1}^2l\cos(\theta_1 - \theta_2))}{2l(2 - \cos(2\theta_1 - 2\theta_2))} \\
\ddot{\theta_2} &= \frac{2\sin(\theta_1 - \theta_2)(2\dot{\theta_1}^2l + 2g\cos\theta_1 + \dot{\theta_2}^2l\cos(\theta_1 - \theta_2))}{2l(2 - \cos(2\theta_1 - 2\theta_2))}
\end{aligned}
\end{cases}
$$
Note that we are note quite finished. Since we will be using a higher (4th) order scheme to numerically
integrate the equations of motion, we need to rewrite the system with four first-order equations as follows.
To be clear, \(\omega_i \) will be the dummy variable representing \(\dot{\theta_i} \).
$$
\boxed{
\begin{aligned}
\omega_1 &= \dot{\theta_1} \\
\omega_2 &= \dot{\theta_2} \\
\dot{\omega_1} &= \frac{-2g\sin\theta_1 - g\sin(\theta_1 - 2\theta_2) - 2\sin(\theta_1 - \theta_2)(\omega_2^2 l + \omega_1^2l\cos(\theta_1 - \theta_2))}{2l(2 - \cos(2\theta_1 - 2\theta_2))} \\
\dot{\omega_2} &= \frac{2\sin(\theta_1 - \theta_2)(2\omega_1^2l + 2g\cos\theta_1 + \omega_2^2l\cos(\theta_1 - \theta_2))}{2l(2 - \cos(2\theta_1 - 2\theta_2))}
\end{aligned}
}
$$
Before employing numerical methods (Runge-Kutta 4 with adaptive stepsizes), I would like to use Hamiltonian formalisms to predict regions where
it is energetically impossible for the second pendulum to 'flip'. Recall that the conjugate momenta of our system was derived above, but here it is again.
$$
\begin{aligned}
p_{\theta_1} &= \frac{\partial L}{\partial \dot{\theta_1}} = ml^2(2\dot{\theta_1} + \dot{\theta_2}\cos(\theta_1 - \theta_2)) \\
p_{\theta_2} &= \frac{\partial L}{\partial \dot{\theta_2}} = ml^2(\dot{\theta_2} + \dot{\theta_1}\cos(\theta_1 - \theta_2))
\end{aligned}
$$
If we were to follow the exact method used in the inspired research paper, we would also use the following equations of motion that we have already found above.
$$
\begin{aligned}
\dot{p_{\theta_1}} &= \frac{\partial L}{\partial \theta_1} = -mgl(2\sin\theta_1) - ml^2(\dot{\theta_1}\dot{\theta_2}\sin(\theta_1 - \theta_2)) \\
\dot{p_{\theta_2}} &= \frac{\partial L}{\partial \theta_2} = -mgl\sin\theta_2 - ml^2(\dot{\theta_1}\dot{\theta_2}\sin(\theta_1 - \theta_2))
\end{aligned}
$$
Since we are limiting the limbs/rods to be motionless at \(t=0 \), the initial Hamiltonian will be as follows.
Note that the Lagrangian is not explicitly dependent on time which means the Hamiltonian is thus converved.
$$ H_0:=T+V=0+V=V \Rightarrow \boxed{H_0 = -mgl(2\cos\theta_1 + \cos\theta_2)}$$
To maximize our \(H_0 \) we let \(\theta_1=\theta_2=\pi \):
$$
H_0 = -mgl(2\cos\theta_1 + \cos\theta_2)> -mgl
$$
Rearranging yields the curve in which it is energetically impossible for either pendulum to flip:
$$
2\cos\theta_1 + \cos\theta_2 = 1\\
\boxed{\theta_2 = \arccos(1-2\cos\theta_1)}
$$
$$ $$
The figure above (what I will call the 'energetics map') displays a large island where neither pendulum has enough energy to 'flip' inside the region. This fact
holds true for late time which will later be confirmed with a numerical analysis. In order to estimate when the second pendulum will 'flip', I will create
a fractal map. For now, we will first familiarize ourselves with the systems behavior.
To start my analysis, I have numerically solved the equations of motion for a simple initial condition: \(\theta_1=\theta_2=\frac{\pi}{2} \).
We note that this particulat initial condition resides outside of the energetically stable island; hence, we can expect the pendulum to flip in finite
time.
$$ $$
$$ $$
Clearly, we can see that this system is chaotic; however, one might question if the pendulum's trajectory would evolve noticeably different if we
give the initial conditions a small perturbation? To prove that this system is chaotic and sensitive to initial conditions, I will re-run the
simulation, but with another double pendulum slightly perturbed. I will keep the same angle \(\theta_1 \); however, I will add the perturbation \(\epsilon = 0.01 \) radians
to \( \theta_2 \).
Hence, \(\theta_2 = \frac{\pi}{2}+\epsilon \). The results are as follows:
$$
$$
$$
$$
As predicted, the pendulums begin on similar trajectories but then quickly (in finite time) deviate paths.
THIS IS CHAOS IN ONE OF ITS PUREST FORMS. As an engineering exercise, I have decided to
'validate' (more of a sanity check) this numerical model in MATLAB by using Autodesk Inventor's
'Dynamic Simulation' environment. For this brief exercise, please note that I have changed the model to match the
inspiration article. To start, I designed each component of the final assembly: two arms and two pins. Then I brought each
component into an assembly file to complete the double pendulum. Finally, I brought the assembly into the dynamic simulation
environment, applied rigid body constraints/forces, set the initial conditions, and then let the simulation run. Here are the
results:
$$$$
$$$$
As we can see, Autodesk Inventor's model behaves similar to our numerical model despite their respective moments of
inertia being different. (**side note: I am currently building this model**). Now that we have established that the double pendulum
is a chaotic system, I will attempt to make order of this chaos using fractal maps. The goal here is to determine if the second
pendulum will flip within some time domain: let's say \( t\in [0,25s] \). If the second pendulum does flip within the domain,
its initial data will get saved. Clearly, this is just a matter of initial conditions; hence, our map will be \(\theta_1 \)
vs \(\theta_2 \). The numerical appraoch is as follows: create two arrays of \(\theta_1\) and \(\theta_2\) values, numerically
simulate the double pendulum for each pair of initial conditions (\(\theta_1\) and \(\theta_2)\), create a contour plot where 'time'
is the domain we are contouring.
The figure above is my fractal map of the system (~10 million grid points using parallel computing). I claim that this map accurately captures the chaotic nature of our system,
but displays well-defined features/patterns that give more valuable insight into predicting results. As we can see, there is a
large island with a pair of buds. This regime (marked in black) is where the pendulum does not flip. The other colors (blue - red)
represent the time at which the second pendulum does flip. Since the buds branching off of the island were not predicted by
the energetics map, I will investigate.
