# Simulation 105: Double Pendulum Modeling with Numerical Integration | by Le Nguyen | Aug, 2023

## Modeling a chaotic system

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.

• 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 systemtheta_0 = np.radians(15) #degrees to radiansg = 9.8 #m/s^2l = 1.0 #momega = np.sqrt(g/l)phi = 0 #for small angletime_span = np.linspace(0,20,300) #simulate for 20s split into 300 time intervalstheta = []for t in time_span:theta.append(simple_pendulum(theta_0, omega, t, phi))#Convert back to cartesian coordinatesx = 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 Integrationtheta_acceleration = -(g/l)*np.sin(theta) #Get accelerationtheta_velocity += time_step*theta_acceleration #Update velocity with accelerationtheta += time_step*theta_velocity #Update angle with angular velocityreturn theta, theta_velocityg = 9.8 #m/s^2l = 1.0 #mtheta = [np.radians(90)] #theta_0theta_velocity = 0 #Start with 0 velocitytime_step = 20/300 #Define a time steptime_span = np.linspace(0,20,300) #simulate for 20s split into 300 time intervalsfor 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)/normalizationreturn theta1_ddot#Get theta2 accelerationdef 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/normalizationreturn theta2_ddot#Update theta1def theta1_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step):#Numerical Integrationtheta1_velocity += time_step*theta1_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g)theta1 += time_step*theta1_velocityreturn theta1, theta1_velocity#Update theta2def theta2_update(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g,time_step):#Numerical Integrationtheta2_velocity += time_step*theta2_acceleration(m1,m2,l1,l2,theta1,theta2,theta1_velocity,theta2_velocity,g)theta2 += time_step*theta2_velocityreturn theta2, theta2_velocity#Run full double pendulumdef 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 xy1 = -l1*np.cos(theta1_list) #Pendulum 1 yx2 = l1*np.sin(theta1_list) + l2*np.sin(theta2_list) #Pendulum 2 xy2 = -l1*np.cos(theta1_list) - l2*np.cos(theta2_list) #Pendulum 2 yreturn x1,y1,x2,y2`
`#Define system parametersg = 9.8 #m/s^2m1 = 1 #kgm2 = 1 #kgl1 = 1 #ml2 = 1 #mtheta1 = np.radians(90)theta2 = np.radians(45)theta1_velocity = 0 #m/stheta2_velocity = 0 #m/stheta1_list = [theta1]theta2_list = [theta2]time_step = 20/300time_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=0images = []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 += 1plt.close()imageio.mimsave(name, images)!rm -r framesdef simple_pendulum(theta_0, omega, t, phi):theta = theta_0*np.cos(omega*t + phi)return theta#parameters of our systemtheta_0 = np.radians(15) #degrees to radiansg = 9.8 #m/s^2l = 1.0 #momega = np.sqrt(g/l)phi = 0 #for small angletime_span = np.linspace(0,20,300) #simulate for 20s split into 300 time intervalstheta = []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_accelerationtheta += time_step*theta_velocityreturn theta, theta_velocityg = 9.8 #m/s^2l = 1.0 #mtheta = [np.radians(90)] #theta_0theta_velocity = 0time_step = 20/300time_span = np.linspace(0,20,300) #simulate for 20s split into 300 time intervalsfor 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 pendulummakeGif(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)/normalizationreturn theta1_ddotdef 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/normalizationreturn theta2_ddotdef 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_velocityreturn theta1, theta1_velocitydef 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_velocityreturn theta2, theta2_velocitydef 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 pendulumg = 9.8 #m/s^2m1 = 1 #kgm2 = 1 #kgl1 = 1 #ml2 = 1 #mtheta1 = np.radians(90)theta2 = np.radians(45)theta1_velocity = 0 #m/stheta2_velocity = 0 #m/stheta1_list = [theta1]theta2_list = [theta2]time_step = 20/300time_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=0images = []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 += 1plt.close()imageio.mimsave("double_pendulum.gif", images)!rm -r frames`