# maths with Python 5: Double Compound Pendulum Chaotic Map

A few weeks ago 6 months ago, instead of holidays, I was spending my time in Cambridge doing measurements there with a good friend, part of the MetMags team. While I’m there, I always enjoy quite a lot taking a look at the students posters in the physics building (Cavendish Laboratory, we work on the Mott building). Because from time to time you find out something new or interesting. And this time I was not disappointed.

I found a poster about the double compound pendulum. And it’s chaotic behaviour. And this time, an interesting plot. This plot In this plot, we put in the x axis the initial value of $\theta_1$ in the range -3,3 and in the y axis, the initial value of $\theta_2$ in the range -3,3. And the colors indicate the amount of elapsed time before the pendulum “flips over”. The colour of each pixel indicates whether either pendulum flips within $10\sqrt(l/g)$ (green), within $100\sqrt(l/g)$ (red), within $1000\sqrt(l/g)$ (purple), within $10000\sqrt(l/g)$ (blue) or in white for other cases. (This description and the values are taken from Wikipedia’s Double Pendulum article).

Ok… nice… can I reproduce that plot? I like chaotic plots a lot, and this is a good challenge. Let’s do it!

First let’s write down the position for the centres of mass of the two bars in the pendulum. Using the first image to define variables, for the first centre of mass we have $\left\{ \begin{matrix} x_1=\frac{l}{2} \sin\theta_1\\ y_1=-\frac{l}{2} \cos\theta_1\end{matrix} \right.$

And for the second one. $\left\{ \begin{matrix} x_2=l(\sin\theta_1+\frac{1}{2} \sin\theta_2) \\ y_2=-l(\cos\theta_1+\frac{1}{2} \cos\theta_2) \end{matrix} \right.$

And now let’s go for the Lagrangian of the system. (For those who didn’t have a classical mechanics course, the Lagrangian is basically where all the information of the system is stored in terms of energy). Since we have to sticks with uniform mass, we have to take into account the rotational energy around the centre of mass. $\begin{matrix} L=Kinetic Energy - Potential Energy \\=\frac{1}{2}m(v_1^2+v_2^2)+\frac{1}{2}I(\dot{\theta}_1^2+\dot{\theta}_2^2)-mg(y_1+y_2) \\=\frac{1}{2}m(\dot{x}_1^2+\dot{y}_1^2+\dot{x}_2^2+\dot{y}_2^2)+\frac{1}{2}I(\dot{\theta}_1^2+\dot{\theta}_2^2)-mg(y_1+y_2) \end{matrix}$

The main feature here is the usage of the moment of inertia and how it is going to affect the rest of the equations. We assume two longitudinal sticks whit uniform mass. $I=\frac{1}{12}ml^2$

It will be interesting to use other geometries and see what properties arise, but solving the problem for a general case will be crazy, so we will do it for this particular case. It’s very important to notice here that the article in the Wikipedia tries to solve the Compound Pendulum, but it makes a mistake using the moment of inertia of the simple pendulum, so the analysis from here is going to be different.

Ok, so now we substitute the previous expression for $x$ and $y$ in the Lagrangian. $L=\frac{1}{2}ml^2[\frac{1}{3}\dot{\theta}_2^2+\frac{4}{3}\dot{\theta}_1^2+\dot{\theta}_1\dot{\theta}_2cos(\theta_1-\theta_2)]+\frac{1}{2}mgl(3\cos\theta_1+\cos\theta_2)$

Nice! And now let’s go for the Euler-Lagrange equations $\frac{d}{dt} \left( \frac{\partial L}{\partial\dot{\theta}_i} \right) -\frac{\partial L}{\partial \theta_i }=0\:\:\:i=1,2$

Wohooo. Hold a second. Solving equations? Me? As my friend say, I’m too senior for that now. So, let’s ask somebody else to do it. Of course looking for the solution is too easy, and nobody is going to do it for free… so, let’s do it in Maple! (Remember that we already use it to show the Rössler system? ) It took me a while to prepare everything, but this is what I made.

First we load the necessary libraries. We load the parameters of this problem.        And the Lagrangian we found before. Because Maple didn´t like too much
subindexes, we substitute theta_1 for theta and theta_2 for mu.    Now let´s define the momenta as    This equations can be inverted to obtain       And now, using the momenta in the Euler-Lagrange equations we get    In this way we have turned the initial second degree equations
into first degree equations.
For solving these equations numerically we use the rkf45 method
and impose an absolute error limit of 10^-16. The time range is from zero to 30 time units.   It is now possible to show the time evolution of both angles.  Another thing we can do, is use the solution to make an animation of the pendulum.
First we need to put the solution into arrays.Remember that our solution is inside
the error bounds only for a range of 30 time units (Here, to limit the size of the
animation, we are going to plot only 10 time units).     And now plot the animation.        Wow Nice!

Ok, so de animation is amazing (a little bit rough, but it is ok, the simulation is correct, it’s just the visualization that can be better), but what is useful here? The reduction of the degree on the equations. These are the first degree equations for angles and moments:    We can write them back in terms of $\theta_1 , \theta_2$ $\frac{d}{dt}\theta_1(t)=\frac{-6(2p_1(t)-3p_2(t)\cos(-\theta_1(t)+\theta_2(t)))}{-16+9\cos(-\theta_1(t)+\theta_2(t))^2}$ $\frac{d}{dt}\theta_2(t)=\frac{6(-8p_2(t)+3p_1(t)\cos(-8\theta_1(t)+\theta_2(t)))}{-16+9\cos(-\theta_1(t)+\theta_2(t))^2}$ $\frac{d}{dt}p_1(t)=-\frac{1}{2}(\frac{d}{dt}\theta_1(t))(\frac{d}{dt}\theta_2(t))\sin(\theta_1-\theta_2)-15\sin(\theta_1(t))$ $\frac{d}{dt}p_2(t)=-\frac{1}{2}(\frac{d}{dt}\theta_1(t))(\frac{d}{dt}\theta_2(t))\sin(-\theta_1+\theta_2)-5\sin(\theta_2(t))$

Great, now we have expressions for the evolution of the pendulum. Since it’s not possible (doesn’t looks like it) to compute when the pendulum is going to flip over, we need to turn these expressions into something more suitable to do some numeric solving.

So, it is now, Python time! Remeber my other posts about Python?

# Raspberry Pi 001: Setup and Run first Arduino-Python project.

It is now time for some heavy stuff.

First attempt. Using Scipy Integrate.

My first attempt at the problem was trying to integrate the equations using the scipy.integrate.odeint. This is my code.

<br /># -*- coding: utf-8 -*-<br />"""<br />Created on Mon Dec 30 21:55:07 2013<br /><br />@author: Hector<br />"""<br />#Load necessary libraries.<br />from numpy import sin, cos, pi, array, arange, zeros_like<br />import matplotlib.pyplot as plt<br />import scipy.integrate as integrate<br />import matplotlib.animation as animation<br /><br />#Here we define the differential equations.<br />#The input is the actual state of the system and<br />#an array with the time points were we want to calculate<br /># the solutions of the system<br />def derivs(state, t):<br /><br />    dydx = zeros_like(state)<br />    dydx =(6*(-8*state+3*state*cos(-state+state)))/(-16+9*(cos(-state+state))*(cos(-state+state)))<br />    dydx =-(6*(2*state-3*state*cos(-state+state)))/(-16+9*(cos(-state+state))*(cos(-state+state)))<br />    dydx =-1/2*dydx*dydx*sin(state-state)-15*sin(state)<br />    dydx =-1/2*dydx*dydx*sin(-state+state)-5*sin(state)<br /><br />    return dydx<br /><br />#This is going to be the initial conditions of our problem<br />mu = -1/4*pi<br />theta = 1.6<br />p1 = 0.0<br />p2 = 0.0<br /><br />#And this the time points where we want the solution.<br />dt = 0.005<br />t = arange(0.0, 30, dt)<br /><br /># initial state<br />state = array([mu, theta, p1, p2])<br /><br /># integrate your ODE using scipy.integrate.<br />y = integrate.odeint(derivs, state, t)<br /><br />#Now the solution are inside array y as columns,<br />#in the same order as we define the state and the diff equations.<br />#We need now to reconstruct the positions of the 2 sticks that made the pendulum.<br />#The first stick will have one end at (0,0) and the other end at:<br />x1 = sin(y[:,1])<br />y1 = -cos(y[:,1])<br />#The second stick will start where the first one ends, and its other part will be at:<br />x2 = sin(y[:,0]) + x1<br />y2 = -cos(y[:,0]) + y1<br /><br />#Now it's time to plot and animate the pendulum.<br />fig = plt.figure()<br />ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2))<br />ax.grid()<br /><br />#This definitions are used to plot the pendulum<br />line, = ax.plot([], [], 'o-', lw=2)<br />punto, = ax.plot([], [], '-', lw=1,markersize=1)<br />time_template = 'time = %.1fs'<br />time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)<br />def init():<br />    punto.set_data([], [])<br />    line.set_data([], [])<br />    time_text.set_text('')<br />    return punto,line, time_text<br />def animate(i):<br />    thisx = [0, x1[i], x2[i]]<br />    thisy = [0, y1[i], y2[i]]<br />    punto.set_data(x2[0:i],y2[0:i])<br />    line.set_data(thisx, thisy)<br />    time_text.set_text(time_template%(i*dt))<br />    return punto,line, time_text<br /><br />#This function animates the plot. It creates a frame for every point in the simulation.<br />#Initialize the ploting functions, used the "animate" to change then, and uses fig as the output.<br />ani = animation.FuncAnimation(fig, animate, arange(1, len(y)),<br />    interval=10, blit=True, init_func=init)<br />ani.save('double_pendulum.mp4',writer='ffmpeg',fps=15)<br />plt.show()<br />

And the output of the program once saved must look like this.

Extra things you will need to run this code.

First of all, close all the command windows and Python interpreters.

Install FFmpeg. To do it, you need to extract the downloaded files, put them on a folder called FFmpeg and add that folder to the system variables path. How to do it is explained here.

Ok, but this attempt failed tremendously when I tried to compute the chaotic map for the initial conditions. Simply it takes too much time to do the integration for all the possible initial conditions.

Second attempt. Try a more direct approach and reduce the computation steps.

I wrote another code, similar to the previous one. Basically, It has a matrix with all the possible initial states. It calculates the next state from that one (next time step), and if the new state is “over the top” that point is eliminated from the matrix and assigned to a map with the colour indicating when it passed over the top.  So, the change in this approach is that now we are not calculating anything extra, just what we need. Here it is the code.

# -*- coding: utf-8 -*-
"""
Created on Sat Jan 04 17:50:11 2014

@author: Hector
"""
from numpy import sin, cos, pi, array, ones, zeros_like,floor,zeros,double,meshgrid,arange
import matplotlib.pyplot as plt
from scipy import integrate
import pickle

#import time
#start_time = time.time() This must be uncommented in case we want to stimate the amount of time needed.

###################################################
#We define the equations as before.
def derivs(t,state):

dydx = zeros_like(state)
dydx =(6*(-8*state+3*state*cos(-state+state)))/(-16+9*(cos(-state+state))*(cos(-state+state)))
dydx =-(6*(2*state-3*state*cos(-state+state)))/(-16+9*(cos(-state+state))*(cos(-state+state)))
dydx =-1/2*dydx*dydx*sin(state-state)-15*sin(state)
dydx =-1/2*dydx*dydx*sin(-state+state)-5*sin(state)

return dydx

######################################################

#These are the parameters of this simulation. It's going to be 400x400 states within the desired range.
nx=4
ny=4
#The simulation is going to run for 96 units of time, and the time step is going to be one.
t_start = 0.0
t_final = 96.0
delta_t = 1
num_steps = floor((t_final - t_start)/delta_t) + 1

#These 4 matrices have the actual state of the system. 2 angles and 2 moments.
mus = [None for x in range(nx*ny)]
thetas = [None for x in range(nx*ny)]
p1s = [None for x in range(nx*ny)]
p2s = [None for x in range(nx*ny)]
#This one will keep the record of the time passed for each estate.
times = t_start*ones(nx*ny)
#The next matrix is going to have the measurement of when the pendulum passed over the top.
p =zeros([nx,ny])
#This is the trick, this list indicates which states need to be computed and which ones are already finished.

compute=range(nx*ny)
compute2=compute
#Load equations and start some variables.
r=integrate.ode(derivs).set_integrator('dop853', method='bdf')
ii=0
jj=0
iii=0
jjj=0
state = array([0, 0, 0, 0])
angle=0

m=len(compute)-1
a=0

#Fill the state matrix's with the initial states.
for kk in range(nx*ny):
mus[kk] = -3+double(iii)/nx*6
thetas[kk] = -3+double(jjj)/ny*6
p1s[kk]=0
p2s[kk]=0
if jjj<(ny-1):
jjj=jjj+1
else:
jjj=0
iii=iii+1

#This is the main loop. It will run until all the states reach the time specified, or until all the states passed over the top.
while max(times)<t_final and m>0:
print(t_final-max(times))
for k in reversed(compute):
state = array([mus[k], thetas[k], p1s[k], p2s[k]])
r.set_initial_value(state, times[k])
r.integrate(r.t + delta_t)
angle=abs(1*(r.y))

#Replace the old state with the new one
mus[k]=r.y
thetas[k]=r.y
p1s[k]=r.y
p2s[k]=r.y

if angle>0.7*pi: #If we want to be rigurous, we can increase the value and make it closer to pi
#This state is now to be included in the color map.
ii=int(k/nx)
jj=k%ny
p[ii,jj]=times[k]
compute2.pop(m) #And eliminated from the computation list.
#Since we are reading values from this list and we ant to eliminate entries, we need 2 lists.

else:
times[k]=times[k]+delta_t
m -= 1
compute=compute2
m=len(compute)-1

#This commend will save the data of the matrix p into the data file.
with open('data_file_name.pickle', 'w') as f:
pickle.dump(p, f)

#This is the command to load some data from a data file.
#with open('data_file_name.pickle') as f:
#with open('data_file_name.pickle') as f:

#elapsed_time = time.time() - start_time
#print(elapsed_time)

fig = plt.figure()
ax1.set_xlabel('Theta')
ax1.set_ylabel('Mu')
ax1.axis([-3,3,-3,3])
x=-3+6.0/(nx-1)*arange(nx)
y=-3+6.0/(nx-1)*arange(ny)
xx,yy=meshgrid(x,y)
plt.pcolormesh(xx,yy,p,cmap='jet')
plt.show()


And after 2 days running the code, this is the result…. Voilá!!! I think for now this code is finish. Am I happy?

No.

The problem is that this map is not yet as good as the one in wikipedia. What to do now?

Our next step is to go for CUDA and parallelization of calculus. But that will be our next post.

Hope you enjoyed this one.