ODE Lab III: Advanced Topics

ODE Lab III: Advanced Topics#

The Double Pendulum Problem#

The double pendulum is a well known example of a non-linear, chaotic system in classical mechanics. It consists of a pendulum with another pendulum attached to its end, resulting in a system with two degrees of freedom. This configuration leads to highly complex and sensitive-dependent dynamics, making the double pendulum an excellent subject for studying chaos theory and non-linear dynamics. Because it is not possible to construct analytical solutions, it is also a great example to numerical integrators.

Double pendulum

The equation of motion for double pendulum is pretty ugly. To set up the equations, we assume:

  • The two arms of the pendulums have the same length \(l\).

  • The mass of each arm is \(m\).

  • The angle between the first and second pendulums, with respect to the vertical axis, are denoted by \(\theta_1\) and \(\theta_2\).

Newton’s second law suggests that we will need to solve a system of two second-order ordinary differential equations (ODEs). Using the methods we learn in the lecture, we can cast the problem into a system of four first-order ODEs.

(433)#\[\begin{align} \frac{d\theta_1}{dt} &= \frac{6}{m l^2}\frac{2 p_1 - 3 \cos(\theta_1 - \theta_2) p_2}{16 - 9 \cos^2(\theta_1 - \theta_2)}\\ \frac{d\theta_2}{dt} &= \frac{6}{m l^2}\frac{8 p_2 - 3 \cos(\theta_1 - \theta_2) p_1}{16 - 9 \cos^2(\theta_1 - \theta_2)}\\ \frac{dp_1}{dt} &= -\frac{1}{2} m l^2 \left(\frac{d\theta_1}{dt} \frac{d\theta_2}{dt}\sin(\theta_1 - \theta_2) + 3\frac{g}{l}\sin\theta_1\right)\\ \frac{dp_2}{dt} &= -\frac{1}{2} m l^2 \left(-\frac{d\theta_1}{dt} \frac{d\theta_2}{dt}\sin(\theta_1 - \theta_2) + \frac{g}{l}\sin\theta_2\right) \end{align}\]

where \(p_1\) and \(p_2\) are called the generalized momenta. (There might be typos in the equation. Please double check.)

It would be impossible to implement this in a few minutes as a hands-on exercise!

Thankfully, we have just learned computational Lagrangian mechanics in class using autodiff and ODE integrator. The Lagrangian of the double pendulum system is:

(434)#\[\begin{align} L &= \frac{1}{2} m\left(v_1^2 + v_2^2\right) + \frac{1}{2}I\left(\dot\theta_1^2 + \dot\theta_2^2\right) - mg (y_1 + y_2) \\ &= \frac{1}{6} m l^2 \left(4\dot\theta_1^2 + \dot\theta_2^2 + 3\dot\theta_1\dot\theta_2\cos(\theta_1 - \theta_2)\right) + \frac{1}{2} mgl \left(3\cos\theta_1 + \cos\theta_2\right) \end{align}\]

We will use this to implement a solver in this hands-on lab.

! pip install jax

Requirement already satisfied: jax in /opt/hostedtoolcache/Python/3.12.10/x64/lib/python3.12/site-packages (0.6.0)
Requirement already satisfied: jaxlib<=0.6.0,>=0.6.0 in /opt/hostedtoolcache/Python/3.12.10/x64/lib/python3.12/site-packages (from jax) (0.6.0)
Requirement already satisfied: ml_dtypes>=0.5.0 in /opt/hostedtoolcache/Python/3.12.10/x64/lib/python3.12/site-packages (from jax) (0.5.1)
Requirement already satisfied: numpy>=1.25 in /opt/hostedtoolcache/Python/3.12.10/x64/lib/python3.12/site-packages (from jax) (2.2.5)
Requirement already satisfied: opt_einsum in /opt/hostedtoolcache/Python/3.12.10/x64/lib/python3.12/site-packages (from jax) (3.4.0)
Requirement already satisfied: scipy>=1.11.1 in /opt/hostedtoolcache/Python/3.12.10/x64/lib/python3.12/site-packages (from jax) (1.15.2)

import jax
jax.config.update("jax_enable_x64", True)

from jax import numpy as np

# Step 1. Copy an ODE Integrator from the class note

def odeint(f, x0, t0, t1, dt=0.1, atol=1e-6, rtol=1e-6):
    ...

# Step 2. Copy ELrhs() from the class note

def ELrhs(L):
    ...

# Step 3. Combine ELrhs() with the ODE Integrator.  See, e.g., path() from class

def path(L, xv0, t0, t1, dt=0.1, atol=1e-6, rtol=1e-6):
    return odeint(ELrhs(L), xv0, t0, t1, dt, atol, rtol) # <- make sure it is compatible with the ODE integrator you choose

# Step 4. Implement the Lagrangian of the double pendulum

def L(x, v):
    x1, x2 = x
    v1, v2 = v
    ...

# Step 5. Use path() to solve the problem

xv0   = np.array([[np.pi/2,np.pi/2], [0.0,0.0]])
T, XV = path(L, xv0, 0.0, 10.0)

X = XV[:,0,:]
V = XV[:,1,:]

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[8], line 4
      1 # Step 5. Use path() to solve the problem
      3 xv0   = np.array([[np.pi/2,np.pi/2], [0.0,0.0]])
----> 4 T, XV = path(L, xv0, 0.0, 10.0)
      6 X = XV[:,0,:]
      7 V = XV[:,1,:]

TypeError: cannot unpack non-iterable NoneType object

# Step 6. Plot the result

...

# Step 7. Animate the result

from matplotlib import animation
from IPython.display import HTML

fig = plt.figure(figsize=(8,8))
ax  = plt.axes(xlim=(-2.5, 2.5), ylim=(-2.5, 2.5))
ax.set_aspect('equal')

line, = ax.plot([], [], 'o-', lw=2)

def init():
    line.set_data([], [])
    return line,

def animate(i):
    th1 = X[i,0]
    th2 = X[i,1]

    x1 =   np.sin(th1)
    y1 = - np.cos(th1)

    x2 =   np.sin(th2)
    y2 = - np.cos(th2)

    line.set_data([0, x1, x1+x2], [0, y1, y1+y2])
    return line,

anim = animation.FuncAnimation(fig, animate, init_func=init, frames=len(T), interval=20, blit=True)

plt.close()

HTML(anim.to_html5_video())
Contents