The pendulum is a classical physics system that we are all familiar with. Be it a grandfather clock or a child on a swing, we have seen the regular, periodic motion of the pendulum. A single pendulum is well defined in classical physics, but the double pendulum (a pendulum attached to the end of another pendulum) is literal chaos. In this article, we are going to build on our intuitive understanding of pendulums and model the chaos of the double pendulum. The physics is interesting and the numerical methods needed are an essential tool in anyone’s arsenal.

In this article we will:

- Learn about harmonic motion and model the behavior of a single pendulum
- Learn the fundamentals of chaos theory
- Model the chaotic behavior of a double pendulum numerically

## Simple Harmonic Motion

We describe the periodic oscillating movement of a pendulum as harmonic motion. Harmonic motion occurs when there is movement in a system that is balanced out by a proportional restoring force in the opposite direction of said movement. We see an example of this in figure 2 where a mass on a spring is being pulled down due to gravity, but this puts energy into the spring which then recoils and pulls the mass back up. Next to the spring system, we see the height of the mass going around in a circle called a phasor diagram which further illustrates the regular motion of the system.

Harmonic motion can be damped (decreasing in amplitude due to drag forces) or driven (increasing in amplitude due to outside force being added), but we will start with the simplest case of indefinite harmonic motion with no outside forces acting on it (undamped motion). This is kind of motion is a good approximation for modeling a single pendulum that swings at a small angle/low amplitude. In this case we can model the motion with equation 1 below.

We can easily put this function into code and simulate a simple pendulum over time.

`def simple_pendulum(theta_0, omega, t, phi):`

theta = theta_0*np.cos(omega*t + phi)

return theta#parameters of our system

theta_0 = np.radians(15) #degrees to radians

g = 9.8 #m/s^2

l = 1.0 #m

omega = np.sqrt(g/l)

phi = 0 #for small angle

time_span = np.linspace(0,20,300) #simulate for 20s split into 300 time intervals

theta = []

for t in time_span:

theta.append(simple_pendulum(theta_0, omega, t, phi))

#Convert back to cartesian coordinates

x = l*np.sin(theta)

y = -l*np.cos(theta) #negative to make sure the pendulum is facing down

## Full Pendulum Motion with Lagrangian Mechanics

A simple small angle pendulum is a good start, but we want to go beyond this and model the motion of a full pendulum. Since we can no longer use small angle approximations it is best to model the pendulum using Lagrangian mechanics. This is an essential tool in physics that switches us from looking at the forces in a system to looking at the energy in a system. We are switching our frame of reference from driving force vs restoring force to kinetic vs potential energy.

The Lagrangain is the difference between kinetic and potential energy given in equation 2.

Substituting in the Kinetic and Potential of a pendulum given in equation 3 yields the Lagrangain for a pendulum seen is equation 4

With the Lagrangian for a pendulum we now describe the energy of our system. There is one last math step to go through to transform this into something that we can build a simulation on. We need to bridge back to the dynamic/force oriented reference from the energy reference using the Euler-Lagrange equation. Using this equation we can use the Lagrangian to get the angular acceleration of our pendulum.

After going through the math, we have angular acceleration which we can use to get angular velocity and angle itself. This will require some numerical integration that will be laid out in our full pendulum simulation. Even for a single pendulum, the non-linear dynamics means there is no analytical solution for solving for *theta, *thus the need for a numerical solution. The integration is quite simple (but powerful), we use angular acceleration to update angular velocity and angular velocity to update *theta *by adding the former quantity to the latter and multiplying this by some time step. This gets us an approximation for the area under the acceleration/velocity curve. The smaller the time step, the more accurate the approximation.

`def full_pendulum(g,l,theta,theta_velocity, time_step):`

#Numerical Integration

theta_acceleration = -(g/l)*np.sin(theta) #Get acceleration

theta_velocity += time_step*theta_acceleration #Update velocity with acceleration

theta += time_step*theta_velocity #Update angle with angular velocity

return theta, theta_velocityg = 9.8 #m/s^2

l = 1.0 #m

theta = [np.radians(90)] #theta_0

theta_velocity = 0 #Start with 0 velocity

time_step = 20/300 #Define a time step

time_span = np.linspace(0,20,300) #simulate for 20s split into 300 time intervals

for t in time_span:

theta_new, theta_velocity = full_pendulum(g,l,theta[-1], theta_velocity, time_step)

theta.append(theta_new)

#Convert back to cartesian coordinates

x = l*np.sin(theta)

y = -l*np.cos(theta)

We have simulated a full pendulum, but this is still a well defined system. It is now time to step into the chaos of the double pendulum.

Chaos, in the mathematical sense, refers to systems that are highly sensitive to their initial conditions. Even slight changes in the system’s start will lead to vastly different behaviors as the system evolves. This perfectly describes the motion of the double pendulum. Unlike the single pendulum, it is not a well behaved system and will evolve in a vastly different way with even slight changes in starting angle.

To model the motion of the double pendulum, we will use the same Lagrangian approach as before (see full derivation).

We will also be using the same numerical integration scheme as before when implementing this equation into code and finding theta.

`#Get theta1 acceleration `

def theta1_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g):

mass1 = -g*(2*m1 + m2)*np.sin(theta1)

mass2 = -m2*g*np.sin(theta1 - 2*theta2)

interaction = -2*np.sin(theta1 - theta2)*m2*np.cos(theta2_velocity**2*l2 + theta1_velocity**2*l1*np.cos(theta1 - theta2))

normalization = l1*(2*m1 + m2 - m2*np.cos(2*theta1 - 2*theta2))theta1_ddot = (mass1 + mass2 + interaction)/normalization

return theta1_ddot

#Get theta2 acceleration

def theta2_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g):

system = 2*np.sin(theta1 - theta2)*(theta1_velocity**2*l1*(m1 + m2) + g*(m1 + m2)*np.cos(theta1) + theta2_velocity**2*l2*m2*np.cos(theta1 - theta2))

normalization = l1*(2*m1 + m2 - m2*np.cos(2*theta1 - 2*theta2))

theta2_ddot = system/normalization

return theta2_ddot

#Update theta1

def theta1_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step):

#Numerical Integration

theta1_velocity += time_step*theta1_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g)

theta1 += time_step*theta1_velocity

return theta1, theta1_velocity

#Update theta2

def theta2_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step):

#Numerical Integration

theta2_velocity += time_step*theta2_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g)

theta2 += time_step*theta2_velocity

return theta2, theta2_velocity

#Run full double pendulum

def double_pendulum(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step,time_span):

theta1_list = [theta1]

theta2_list = [theta2]

for t in time_span:

theta1, theta1_velocity = theta1_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step)

theta2, theta2_velocity = theta2_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step)

theta1_list.append(theta1)

theta2_list.append(theta2)

x1 = l1*np.sin(theta1_list) #Pendulum 1 x

y1 = -l1*np.cos(theta1_list) #Pendulum 1 y

x2 = l1*np.sin(theta1_list) + l2*np.sin(theta2_list) #Pendulum 2 x

y2 = -l1*np.cos(theta1_list) - l2*np.cos(theta2_list) #Pendulum 2 y

return x1,y1,x2,y2

`#Define system parameters`

g = 9.8 #m/s^2m1 = 1 #kg

m2 = 1 #kg

l1 = 1 #m

l2 = 1 #m

theta1 = np.radians(90)

theta2 = np.radians(45)

theta1_velocity = 0 #m/s

theta2_velocity = 0 #m/s

theta1_list = [theta1]

theta2_list = [theta2]

time_step = 20/300

time_span = np.linspace(0,20,300)

x1,y1,x2,y2 = double_pendulum(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step,time_span)

We’ve finally done it! We have successfully modeled a double pendulum, but now it’s time to observe some chaos. Our final simulation will be of two double pendulums with slightly different starting condition. We will set one pendulum to have a *theta 1 *of 90 degrees and the other to have a *theta 1* of 91 degrees. Let’s see what happens.

We can see that both pendulums start off with similar trajectories but quickly diverge. This is what we mean when we say chaos, even a 1 degree difference in angle cascades into vastly different end behavior.

In this article we learned about pendulum motion and how to model it. We started from the simplest harmonic motion model and built up to the complex and chaotic double pendulum. Along the way we learned about the Lagrangian, chaos, and numerical integration.

The double pendulum is the simplest example of a chaotic system. These systems exist everywhere in our world from population dynamics, climate, and even billiards. We can take the lessons we have learned from the double pendulum and apply them whenever we encounter a chaotic systems.

## Key Take Aways

- Chaotic systems are very sensitive to initial conditions and will evolve in vastly different ways with even slight changes to their start.
- When dealing with a system, especially a chaotic system, is there another frame of reference to look at it that makes it easier to work with? (Like the force reference frame to the energy reference frame)
- When systems get too complicated we need to implement numerical solutions to solve them. These solutions are simple but powerful and offer good approximations to the actual behavior.

All figures used in this article were either created by the author or are from Math Images and full under the GNU Free Documentation License 1.2

Classical Mechanics, John Taylor https://neuroself.files.wordpress.com/2020/09/taylor-2005-classical-mechanics.pdf

## Simple Pendulum

`def makeGif(x,y,name):`

!mkdir framescounter=0

images = []

for i in range(0,len(x)):

plt.figure(figsize = (6,6))

plt.plot([0,x[i]],[0,y[i]], "o-", color = "b", markersize = 7, linewidth=.7 )

plt.title("Pendulum")

plt.xlim(-1.1,1.1)

plt.ylim(-1.1,1.1)

plt.savefig("frames/" + str(counter)+ ".png")

images.append(imageio.imread("frames/" + str(counter)+ ".png"))

counter += 1

plt.close()

imageio.mimsave(name, images)

!rm -r frames

def simple_pendulum(theta_0, omega, t, phi):

theta = theta_0*np.cos(omega*t + phi)

return theta

#parameters of our system

theta_0 = np.radians(15) #degrees to radians

g = 9.8 #m/s^2

l = 1.0 #m

omega = np.sqrt(g/l)

phi = 0 #for small angle

time_span = np.linspace(0,20,300) #simulate for 20s split into 300 time intervals

theta = []

for t in time_span:

theta.append(simple_pendulum(theta_0, omega, t, phi))

x = l*np.sin(theta)

y = -l*np.cos(theta) #negative to make sure the pendulum is facing down

## Pendulum

`def full_pendulum(g,l,theta,theta_velocity, time_step):`

theta_acceleration = -(g/l)*np.sin(theta)

theta_velocity += time_step*theta_acceleration

theta += time_step*theta_velocity

return theta, theta_velocityg = 9.8 #m/s^2

l = 1.0 #m

theta = [np.radians(90)] #theta_0

theta_velocity = 0

time_step = 20/300

time_span = np.linspace(0,20,300) #simulate for 20s split into 300 time intervals

for t in time_span:

theta_new, theta_velocity = full_pendulum(g,l,theta[-1], theta_velocity, time_step)

theta.append(theta_new)

#Convert back to cartesian coordinates

x = l*np.sin(theta)

y = -l*np.cos(theta)

#Use same function from simple pendulum

makeGif(x,y,"pendulum.gif")

## Double Pendulum

`def theta1_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g):`

mass1 = -g*(2*m1 + m2)*np.sin(theta1)

mass2 = -m2*g*np.sin(theta1 - 2*theta2)

interaction = -2*np.sin(theta1 - theta2)*m2*np.cos(theta2_velocity**2*l2 + theta1_velocity**2*l1*np.cos(theta1 - theta2))

normalization = l1*(2*m1 + m2 - m2*np.cos(2*theta1 - 2*theta2))theta1_ddot = (mass1 + mass2 + interaction)/normalization

return theta1_ddot

def theta2_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g):

system = 2*np.sin(theta1 - theta2)*(theta1_velocity**2*l1*(m1 + m2) + g*(m1 + m2)*np.cos(theta1) + theta2_velocity**2*l2*m2*np.cos(theta1 - theta2))

normalization = l1*(2*m1 + m2 - m2*np.cos(2*theta1 - 2*theta2))

theta2_ddot = system/normalization

return theta2_ddot

def theta1_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step):

theta1_velocity += time_step*theta1_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g)

theta1 += time_step*theta1_velocity

return theta1, theta1_velocity

def theta2_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step):

theta2_velocity += time_step*theta2_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g)

theta2 += time_step*theta2_velocity

return theta2, theta2_velocity

def double_pendulum(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step,time_span):

theta1_list = [theta1]

theta2_list = [theta2]

for t in time_span:

theta1, theta1_velocity = theta1_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step)

theta2, theta2_velocity = theta2_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step)

theta1_list.append(theta1)

theta2_list.append(theta2)

x1 = l1*np.sin(theta1_list)

y1 = -l1*np.cos(theta1_list)

x2 = l1*np.sin(theta1_list) + l2*np.sin(theta2_list)

y2 = -l1*np.cos(theta1_list) - l2*np.cos(theta2_list)

return x1,y1,x2,y2

`#Define system parameters, run double pendulum`

g = 9.8 #m/s^2m1 = 1 #kg

m2 = 1 #kg

l1 = 1 #m

l2 = 1 #m

theta1 = np.radians(90)

theta2 = np.radians(45)

theta1_velocity = 0 #m/s

theta2_velocity = 0 #m/s

theta1_list = [theta1]

theta2_list = [theta2]

time_step = 20/300

time_span = np.linspace(0,20,300)

for t in time_span:

theta1, theta1_velocity = theta1_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step)

theta2, theta2_velocity = theta2_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step)

theta1_list.append(theta1)

theta2_list.append(theta2)

x1 = l1*np.sin(theta1_list)

y1 = -l1*np.cos(theta1_list)

x2 = l1*np.sin(theta1_list) + l2*np.sin(theta2_list)

y2 = -l1*np.cos(theta1_list) - l2*np.cos(theta2_list)

`#Make Gif`

!mkdir framescounter=0

images = []

for i in range(0,len(x1)):

plt.figure(figsize = (6,6))

plt.figure(figsize = (6,6))

plt.plot([0,x1[i]],[0,y1[i]], "o-", color = "b", markersize = 7, linewidth=.7 )

plt.plot([x1[i],x2[i]],[y1[i],y2[i]], "o-", color = "b", markersize = 7, linewidth=.7 )

plt.title("Double Pendulum")

plt.xlim(-2.1,2.1)

plt.ylim(-2.1,2.1)

plt.savefig("frames/" + str(counter)+ ".png")

images.append(imageio.imread("frames/" + str(counter)+ ".png"))

counter += 1

plt.close()

imageio.mimsave("double_pendulum.gif", images)

!rm -r frames