Simple non-linear differential equations can have complicated and chaotic solutions. In this notebook we will study the chaotic motion of the double pendulum. Even though the motion is chaotic, phase-space plots reveals a beautiful underlying order in the chaos.

The double pendulum is shown in figure 1. It consists of two simple pendulums, where the mass of the first pendulum is the pivot point for the second. We refer you to our notebook on the simple pendulum for a quick introduction. Let $m_i$ and $L_i$ be the mass and the length of the $i$th simple pendulum, respectively. $\theta_1$ is the angle between the vertical direction and the first pendulum, and $\theta_2$ is the angle between the vertical direction and the second, as shown in the figure. Note that there is only two degrees of freedom, $\theta_1$ and $\theta_2$. However, this is enough to obtain a chaotic motion.

**Figure 1.** The double pendulum with masses $m_1$ and $m_2$, displacement angles $\theta_1$ and $\theta_2$ and lengths $L_1$ and $L_2$.

We start off by setting some plotting constants and importing needed packages.

In [7]:

```
from scipy.integrate import odeint
import matplotlib.pyplot as plt
from numpy import cos, sin, arange, pi
import matplotlib.cm as cm
%matplotlib inline
figsize = 6
dpi = 600
g = 9.81 # [m/s^2]. Gravitational acceleration
```

The equations of motion can be written as (see appendix for the derivation)

\begin{align} \dot\omega_1 &= \frac{1}{L_1\xi}\left[L_1m_2c_{12}s_{12}\omega_1^2 + L_2m_2s_{12}\omega_2^2 - m_2gc_{12}s_2 + (m_1+m_2)gs_1 \right],\\ \dot\omega_2 &= \frac{1}{L_2\xi}\left[L_2m_2c_{12}s_{12}\omega_2^2 + L_1(m_1+m_2)s_{12}\omega_1^2+(m_1+m_2)gs_1c_{12} - (m_1+m_2)gs_2 \right],\\ \omega_1 &= \dot\theta_1,\\ \omega_2 &= \dot\theta_2, \end{align}where we have defined $c_{12} \equiv \cos(\theta_1-\theta_2)$, $s_{12} \equiv \sin(\theta_1-\theta_2)$, $s_i \equiv \sin(\theta_i)$ and $\xi \equiv c_{12}^2m_2-m_1-m_2$ in order to simplify the notation. Note that we have written the equations as four coupled first order ordinary differential equations, which means that we can easily solve the corresponding initial value problem using e.g. the Euler method or the 4th order Runge-Kutta method. In this notebook we will be using the integrator `odeint`

from `scipy`

. This function solves the initial value problem for stiff or non-stiff systems of first order ode-s. Check out [1] to learn how it works. To use `odeint`

we need to create a function that evaluates the right hand side of the equations of motion. It will also be convinient to define a function that transforms $\omega_1$,$\omega_2$, $\theta_1$ and $\theta_2$ into cartesian coordinates $\mathbf{v}_1$, $\mathbf{v}_2$, $\mathbf{x}_1$ and $\mathbf{x}_2$.

In [8]:

```
def RHS(z, t, L1, L2, m1, m2, g):
""" Return the right hand side of the
ordinary differential equation describing
the double pendulum.
"""
theta1, w1, theta2, w2 = z
cos12 = cos(theta1 - theta2)
sin12 = sin(theta1 - theta2)
sin1 = sin(theta1)
sin2 = sin(theta2)
xi = cos12**2*m2 - m1 - m2
w1dot = ( L1*m2*cos12*sin12*w1**2 + L2*m2*sin12*w2**2
- m2*g*cos12*sin2 + (m1 + m2)*g*sin1)/(L1*xi)
w2dot = -( L2*m2*cos12*sin12*w2**2 + L1*(m1 + m2)*sin12*w1**2
+ (m1 + m2)*g*sin1*cos12 - (m1 + m2)*g*sin2 )/(L2*xi)
return w1, w1dot, w2, w2dot
def to_cartesian(theta1, w1, theta2, w2, L1, L2):
""" Transforms theta and omega to cartesian coordinates
and velocities x1, y1, x2, y2, vx1, vy1, vx2, vy2
"""
x1 = L1 * sin(theta1)
y1 = -L1 * cos(theta1)
x2 = x1 + L2 * sin(theta2)
y2 = y1 - L2 * cos(theta2)
vx1 = L1*cos(theta1)*w1
vy1 = L1*sin(theta1)*w1
vx2 = vx1 + L2*cos(theta2)*w2
vy2 = vy1 + L2*sin(theta2)*w2
return x1, y1, x2, y2, vx1, vy1, vx2, vy2
```

We now need to define the parameters in the setup. Play around with initial conditions and different parameters!

Here, we make the (arbitrary) choice $2L_1=L_2$ and $m_1=3m_2$, and we will let the pendulum swing for 50 seconds.

In [9]:

```
L1, L2 = 1., 2.
m1, m2 = 3., 1.
z0 = [pi/2, 0, pi/2, 0]
tmax, dt = 50, 0.01
t = arange(0, tmax+dt, dt)
```

We are now ready to perform the simulation and transform the results to cartesian coordinates.

In [10]:

```
# Perform simulation
z = odeint(RHS, z0, t, args=(L1, L2, m1, m2, g))
# Extract result
theta1, w1, theta2, w2 = z[:,0], z[:,1], z[:,2], z[:,3]
x1, y1, x2, y2, vx1, vy1, vx2, vy2 = to_cartesian(theta1, w1, theta2, w2, L1, L2)
```

We can now plot the results and create an animation! To this end, we created the functions `plot_position()`

and `create_animation()`

, which are shown in the appendix.

In [11]:

```
plot_position(x1, y1, x2, y2, theta1, theta2, t)
```