At last, we will tackle the infamous double-pendulum. We have a mass $m_1$ connected to a massles, rigid rod of length $l_1$, and displaced by some angle $\theta_1$ from the vertical. Attached to $m_1$ is another massless rigid rod of length $l_2$ with a mass $m_2$ attached to its other end, that is displaced from the vertical by some angle $\theta_2$. We will use Lagrangian mechanics to solve for the system's equations of motion.
It is best to start in cartesian coordinates, then convert to polar coordinates. We begin by looking at $m_1$. We write out $m_1$'s position and velocity.
$$
x_1 = l_1\sin\theta_1
$$
$$
\dot{x}_1 = l_1\dot{\theta}_1\sin\theta_1
$$
$$
y_1 = -l_1\cos\theta_1
$$
$$
\dot{y}_1 = l_1\dot{\theta}_1\sin\theta_1
$$
Now we can find the kinetic energy of $m_1$.
$$
T_1 = \frac{1}{2}m_1\left( \dot{x}_1^2 + \dot{y}_1^2 \right)
$$
$$
\Rightarrow T_1 = \frac{1}{2}m_1 l_1^2\dot{\theta}_1^2\cos^2\theta_1 + \frac{1}{2}m_1 l_1^2\dot{\theta}_1^2\sin^2\theta_1
$$
$$
\Rightarrow T_1 = \frac{1}{2}m_1 l_1^2\dot{\theta}_1^2
$$
Where we have used the trig identity
\begin{equation}
\cos^2\alpha + \sin^2\alpha = 1
\end{equation}
Now we find the potential energy of $m_1$, which only has graviational potential energy.
$$
U_1 = m_1 g y_1 = -m_1gl_1\cos\theta_1
$$
We need to do the same for $m_2$ now. First, we find its position and velocity.
$$
x_2 = l_1\sin\theta_1 + l_2\sin\theta_2
$$
$$
\dot{x}_2 = l_1\dot{\theta}_1\cos\theta_1 + l_2\dot{\theta}_2\cos\theta_2
$$
$$
y_2 = -l_1\cos\theta_1 - l_2\cos\theta_2
$$
$$
\dot{y}_2 = l_1\dot{\theta}_1\sin\theta_1 + l_2\dot{\theta}_2\sin\theta_2
$$
So, $m_2$ has kinetic energy
$$
T_2 = \frac{1}{2}m_2\left( \dot{x}_2^2 + \dot{y}_2^2 \right)
$$
$$
\Rightarrow T_2 = \frac{1}{2}m_2\left[ \left(l_1\dot{\theta}_1\cos\theta_1 + l_2\dot{\theta}_2\cos\theta_2 \right)^2 + \left( l_1\dot{\theta}_1 + l_2\dot{\theta}_2\sin\theta_2 \right)^2\right]
$$
Squaring these binomials, using identity (1), and the following identity
\begin{equation}
\cos(\alpha - \beta) = \cos\alpha\cos\beta + \sin\alpha\sin\beta
\end{equation}
yields our final expression for the kinetic energy of $m_2$:
$$
T_2 = \frac{1}{2}m_2\left[l_1^2\dot{\theta}_1^2 + l_2^2\dot{\theta}_2^2 + 2l_1l_2\dot{\theta}_1\dot{\theta}_2\cos(\theta_1-\theta_2) \right]
$$
Now we find the potential energy of $m_2$, which again only has gravitational potential energy.
$$
U_2=m_2g(y_1+y_2) = -m_2g(l_1\cos\theta_1 + l_2\cos\theta_2)
$$
Finally, we might write out an expression for this system's Lagrangian.
$$
\mathcal{L}=T-U = (T_1 + T_2) - (U_1 + U_2)
$$
\begin{multline}
\Rightarrow \mathcal{L} = \frac{1}{2}m_1l_1^2\dot{\theta}_1^2 + \frac{1}{2}m_2\left[l_1^2\dot{\theta}_1^2 + l_2^2\dot{\theta}_2^2 + 2l_1l_2\dot{\theta}_1\dot{\theta}_2\cos(\theta_1-\theta_2) \right]\\
+m_1gl_1\cos\theta_1 + m_2g(l_1\cos\theta_1 + l_2\cos\theta_2)
\end{multline}
We might now invoke the Euler-Lagrange equation, using $q$ as a generalized coordinate.
\begin{equation}
\frac{\partial \mathcal{L}}{\partial q} - \frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial\dot{q}}\right)=0
\end{equation}
Now we start taking derivatives -- if you are uncomfortable with the chain rule, this is a fantastic problem to work through. We start with $\theta_1$.
\begin{multline}
\frac{\partial \mathcal{L}}{\partial \theta_1}=m_2l_1l_2\dot{\theta}_1\dot{\theta}_2\sin(\theta_1-\theta_2)-m_1gl_1\sin\theta_1-m_2gl_1\sin\theta_1
\end{multline}
\begin{multline*}
\frac{\partial \mathcal{L}}{\partial \dot{\theta}_1}=m_1l_1^2\dot{\theta}_1 + m_2l_1^2\dot{\theta}_1 +l_1l_2\dot{\theta}_2\cos(\theta_1-\theta_2)
\end{multline*}
\begin{multline}
\frac{d}{dt}\left( \frac{\partial \mathcal{L}}{\partial \dot{\theta}_1} \right)=m_1l_1^2\ddot{\theta}_1 + m_2l_1^2\ddot{\theta}_1 + m_2l_1l_2\ddot{\theta}_2\cos(\theta_1-\theta_2)\\
-m_2l_1l_2\dot{\theta}_2(\dot{\theta}_1 - \dot{\theta}_2)\sin(\theta_1-\theta_2)
\end{multline}
Combining equations (4), (5), and (6) yields our first equation of motion.
\begin{multline}
\ddot{\theta}_1\left(m_1l_1^2 + m_2l_1^2 \right) = m_2l_1l_2\dot{\theta}_2(\dot{\theta}_1-\dot{\theta}_2)\sin(\theta_1-\theta_2)\\
-m_2l_1l_2\ddot{\theta}_2\cos(\theta_1-\theta_2) - m_2l_1l_2\dot{\theta}\dot{\theta}_2\sin(\theta_1-\theta_2)\\
-m_1gl_1\sin\theta_1-m_2gl_1\sin\theta_1
\end{multline}
Well... that is kind of gross to look at. Now we turn our attention to $\theta_2$.
\begin{multline}
\frac{\partial\mathcal{L}}{\partial\theta_2} = m_2l_1l_2\dot{\theta}_1\dot{\theta}_2\sin(\theta_1-\theta_2)-m_2gl_2\sin\theta_2
\end{multline}
\begin{multline*}
\frac{\partial\mathcal{L}}{\partial\dot{\theta}_2}=m_2l_2^2\ddot{\theta}_2+m_2l_1l_2\dot{\theta}_1\cos(\theta_1-\theta_2)
\end{multline*}
\begin{multline}
\frac{d}{dt}\left(\frac{\partial\mathcal{L}}{\partial\dot{\theta}_2}\right)=m_2l_2^2\ddot{\theta}_2+m_2l_1l_2\ddot{\theta}_1\cos(\theta_1-\theta_2)
-m_2l_1l_2\dot{\theta}_1(\dot{\theta}_1-\dot{\theta}_2)\sin(\theta_1-\theta_2)
\end{multline}
Combining equations (4), (8), and (9) yields our second equation of motion.
\begin{multline}
\ddot{\theta}_2m_2l_2^2 = m_2l_1l_2\dot{\theta}_1(\dot{\theta}_1-\dot{\theta}_2)\sin(\theta_1-\theta_2)-m_2gl_2\sin\theta_2 - m_2l_1l_2\ddot{\theta}_1\cos(\theta_1-\theta_2)
\end{multline}
Equations (7) and (10) are our equations of motion. However, we have an issue. Equation (7) depends on both $\ddot{\theta}_1$ and $\ddot{\theta}_2$, as does equation (10). Luckily we can just simultaneously solve these equations to fix this issue, unluckliy, this is an algebraic nightmare. I used Mathematica to do this part for me, as it is much better at book-keeping than myslef. I also used Mathematica to check my work as I went. You may download my Mathematica code here:
Mathematica notebook
Now we use the usual trick to turn these two second order differential equations into 4 coupled first order equations. We let
$$
\dot{\theta}_1 = \omega_1 \hspace{1.6cm} \text{and} \hspace{1.6cm} \dot{\theta}_2 = \omega_2
$$
such that
\begin{multline}
\dot{\omega}_1= [-m_1g\sin\theta_1 - m_2g\cos\theta_2\sin(\theta_1 - \theta_2) - m_2l_1\omega_1^2\cos(\theta_1-\theta_2)\sin(\theta_1-\theta_2)\\
-m_2l_2\omega_2^2\sin(\theta_1-\theta_2)]\left[\frac{1}{l_1m_1 + l_1m_2\sin^2(\theta_1-\theta_2)}\right]
\end{multline}
\begin{multline}
\dot{\omega}_2 = [m_1g\cos\theta_1\sin(\theta_1-\theta_2) + m_1l_1\omega_1^2\sin(\theta_1-\theta_2) + m_2g\cos\theta_1\sin(\theta_1-\theta_2) \\
+ m_2l_2\omega_2^2\cos(\theta_1-\theta_2)\sin(\theta_1-\theta_2)\\
+ m_2l_1\omega_1^2\sin(\theta_1-\theta_2)] \left[\frac{1}{l_2m_1 + l_2m_2\sin^2(\theta_1-\theta_2)}\right]
\end{multline}
Wonderful. I solved these equations using the 4th order Runge-Kutta method. See the code block below.
I think coding is gross!
# solve the simple pendulum using fourth order Runge-Kutta method
# Import packages needed
import matplotlib.animation as animation
from numpy import sqrt, sin, cos, arange, pi, append, array, floor
from pylab import plot, xlabel, ylabel, title, show, axhline, savefig, subplots_adjust,\
figure, xlim, rcParams, rc, rc_context, subplot, tight_layout, axvline, xlim, ylim, scatter
# define a function to calculate our slopes at a given position
# eta = (theta1, omega1, theta2, omega2)
def f(eta):
theta1 = eta[0]
omega1 = eta[1]
theta2 = eta[2]
omega2 = eta[3]
f_theta1 = omega1
f_theta2 = omega2
f_omega1 = -(g*sin(theta1)*m1 + sin(theta1-theta2)*m2*(g*cos(theta2) \
+ cos(theta1-theta2)*l1*omega1**2 + l2*omega2**2)) \
/ (l1*m1 + l1*m2*sin(theta1-theta2)**2)
f_omega2 = (sin(theta1-theta2)*(m1*(g*cos(theta1)+l1*omega1**2)\
+ m2*(g*cos(theta1) + l1*omega1**2 + cos(theta1-theta2)*l2*omega2**2))) \
/ (l2*m1 + l2*m2*sin(theta1-theta2)**2)
return(array([f_theta1,f_omega1,f_theta2,f_omega2],float))
# define a function that takes initial angles in as a parameter
# and outputs array of angles and velocities of both masses
def doublePendulumGuy(theta1_initial_deg, theta2_initial_deg):
omega1_initial_deg = 0.0 # initial angular speed
omega2_initial_deg = 0.0 # initial angular speed
theta1_initial = theta1_initial_deg*pi/180 # convert initial angle into degrees
omega1_initial = omega1_initial_deg*pi/180 # convert initial anglular speed into degrees
theta2_initial = theta2_initial_deg*pi/180 # convert initial angle into degrees
omega2_initial = omega2_initial_deg*pi/180 # convert initial anglular speed into degrees
# set up a domain (time interval of interest)
a = 0.0 # interval start
b = 100.0 # interval end
dt = 0.002 # timestep
t_points = arange(a,b,dt) # array of times
# initial conditions eta = (theta1, omega1, theta2, omega2)
eta = array([theta1_initial, omega1_initial, theta2_initial, omega2_initial],float)
# create empty sets to update with values of interest, then invoke Runge-Kutta
theta1_points = []
omega1_points = []
theta2_points = []
omega2_points = []
for t in t_points:
theta1_points.append(eta[0])
omega1_points.append(eta[1])
theta2_points.append(eta[2])
omega2_points.append(eta[3])
k1 = dt*f(eta)
k2 = dt*f(eta + 0.5*k1)
k3 = dt*f(eta + 0.5*k2)
k4 = dt*f(eta + k3)
eta += (k1 + 2*k2 + 2*k3 + k4)/6
return(theta1_points, omega1_points, theta2_points, omega2_points, t_points)
# set up the parameters of our situation
m1 = 5.00 # mass of bob 1
m2 = 3.50 # mass of bob 2
l1 = 0.85 # length of rod 1
l2 = 1.20 # length of rod 2
g = 9.81 # acceleration due to gravity
# call the function and store the arrays of data
theta1_points, omega1_points, theta2_points, omega2_points,\
t_points = doublePendulumGuy(179.00000000,179.000000000)
theta1_points0, omega1_points0, theta2_points0, omega2_points0, \
t_points = doublePendulumGuy(179.00000001,179.000000001)
# create lists of x and y values for bob2
x_points = []
y_points = []
x_points0 = []
y_points0 = []
for i in range(len(t_points)):
x_new = l1*sin(theta1_points[i]) + l2*sin(theta2_points[i])
y_new = -l1*cos(theta1_points[i]) - l2*cos(theta2_points[i])
x_points.append(x_new)
y_points.append(y_new)
x_new0 = l1*sin(theta1_points0[i]) + l2*sin(theta2_points0[i])
y_new0 = -l1*cos(theta1_points0[i]) - l2*cos(theta2_points0[i])
x_points0.append(x_new0)
y_points0.append(y_new0)
# animate our results
# start with styling options
rcParams.update({'font.size': 18})
rc('axes', linewidth=2)
with rc_context({'axes.edgecolor':'white', 'xtick.color':'white', \
'ytick.color':'white', 'figure.facecolor':'darkslategrey',\
'axes.facecolor':'darkslategrey','axes.labelcolor':'white',\
'axes.titlecolor':'white'}):
# set up a figure that will be ultimately be 2x2
fig = figure(figsize=(10,10))
# subplot for animation of pendulum
ax_pend = subplot(1,1,1, aspect='equal')
# get rid of axis ticks
ax_pend.tick_params(axis='both', colors="darkslategrey")
xlim(-1.1*(l1+l2),1.1*(l1+l2))
ylim(-1.1*(l1+l2),1.1*(l1+l2))
### finally we animate ###
# create a list to input images in for each time step
ims = []
index = 0
# only show the first 80seconds or so in the gif
while index <= len(t_points)-1:
ln1, = ax_pend.plot([0,l1*sin(theta1_points[index])],\
[0,-l1*cos(theta1_points[index])],\
color='white',lw=3,zorder=99)
bob1, = ax_pend.plot(l1*sin(theta1_points[index]),\
-l1*cos(theta1_points[index]),'o',\
markersize=22,color="mediumvioletred",zorder=100)
ln2, = ax_pend.plot([l1*sin(theta1_points[index]),l1*sin(theta1_points[index])+\
l2*sin(theta2_points[index])],\
[-l1*cos(theta1_points[index]),-l1*cos(theta1_points[index])\
-l2*cos(theta2_points[index])], color='white',lw=3,zorder=99)
bob2, = ax_pend.plot(l1*sin(theta1_points[index])+l2*sin(theta2_points[index]),\
-l1*cos(theta1_points[index])-l2*cos(theta2_points[index]),'o',\
markersize=22,color="coral",zorder=100)
ln10, = ax_pend.plot([0,l1*sin(theta1_points0[index])],\
[0,-l1*cos(theta1_points0[index])],\
color='k',lw=3,zorder=98)
bob10, = ax_pend.plot(l1*sin(theta1_points0[index]),\
-l1*cos(theta1_points0[index]),'o',\
markersize=22,color="purple",zorder=99)
ln20, = ax_pend.plot([l1*sin(theta1_points0[index]),l1*sin(theta1_points0[index])+\
l2*sin(theta2_points0[index])],\
[-l1*cos(theta1_points0[index]),-l1*cos(theta1_points0[index])\
-l2*cos(theta2_points0[index])], color='k',lw=3,zorder=98)
bob20, = ax_pend.plot(l1*sin(theta1_points0[index])+l2*sin(theta2_points0[index]),\
-l1*cos(theta1_points0[index])-l2*cos(theta2_points0[index]),'o',\
markersize=22,color="green",zorder=99)
if index > 1000:
trail, = ax_pend.plot(x_points[index-990:index],y_points[index-990:index], \
color="cyan",lw=0.8,zorder=20)
trail0, = ax_pend.plot(x_points0[index-990:index],y_points0[index-990:index], \
color="magenta",lw=0.8)
else:
trail, = ax_pend.plot(x_points[:index],y_points[:index], \
color="cyan",lw=0.8,zorder=20)
trail0, = ax_pend.plot(x_points0[:index],y_points0[:index], \
color="magenta",lw=0.8)
# add pictures to ims list
ims.append([ln1, bob1, ln2, bob2, trail,\
ln10, bob10, ln20, bob20, trail0])
# only show every 6 frames
index += 6
# save animations
ani = animation.ArtistAnimation(fig, ims, interval=100)
writervideo = animation.FFMpegWriter(fps=60)
ani.save('./179_179_doublePend.mp4', writer=writervideo)
How hypnotizing these things are to watch. This is a choatic system. Note that some of the driven pendulum examples I showed seemed to also be chaotic in behavior. This means that the system never reaches periodicity. A defining feature of chaos is that small changes in initial conditions will cause vast differences in the outcome of the system's motion. The code above is actually set up to produce the video below. It overlays two double pendulums atop one another. One has initial conditions $\theta_1=\theta_2=179.000000000$ degrees, while the other has intitial conditions $\theta_1=\theta_2=179.000000001$ degrees. A difference of a billionth of a degree. Let us see how long it takes for their motions to differ from one another, and by how much they differ.
You should check out the elastic double pendulum next.