The pendulum is a classical physics system that we’re all acquainted with. Be it a grandfather clock or a baby on a swing, we’ve seen the regular, periodic motion of the pendulum. A single pendulum is well defined in classical physics, however the double pendulum (a pendulum attached to the tip of one other pendulum) is literal chaos. In this text, we’re going to construct 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 important tool in anyone’s arsenal.
In this text we are going to:
- Study harmonic motion and model the behavior of a single pendulum
- Learn the basics of chaos theory
- Model the chaotic behavior of a double pendulum numerically
Easy Harmonic Motion
We describe the periodic oscillating movement of a pendulum as harmonic motion. Harmonic motion occurs when there may be movement in a system that’s balanced out by a proportional restoring force in the wrong way of said movement. We see an example of this in figure 2 where a mass on a spring is being pulled down attributable 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 peak of the mass going around in a circle called a phasor diagram which further illustrates the regular motion of the system.
Harmonic motion may be damped (decreasing in amplitude attributable to drag forces) or driven (increasing in amplitude attributable to outside force being added), but we are going to start with the only case of indefinite harmonic motion with no outside forces acting on it (undamped motion). That is type of motion is a superb approximation for modeling a single pendulum that swings at a small angle/low amplitude. On this case we are able to model the motion with equation 1 below.
We are able to easily put this function into code and simulate a straightforward 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 ensure the pendulum is facing down
Full Pendulum Motion with Lagrangian Mechanics
An easy small angle pendulum is a superb start, but we would like to transcend this and model the motion of a full pendulum. Since we are able to not use small angle approximations it’s best to model the pendulum using Lagrangian mechanics. That is an important tool in physics that switches us from taking a look at the forces in a system to taking a look at the energy in a system. We’re 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 within 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 may be one last math step to undergo to rework this into something that we are able to construct a simulation on. We want to bridge back to the dynamic/force oriented reference from the energy reference using the Euler-Lagrange equation. Using this equation we are able to use the Lagrangian to get the angular acceleration of our pendulum.
After going through the mathematics, we’ve angular acceleration which we are able to use to get angular velocity and angle itself. This can require some numerical integration that shall be specified by our full pendulum simulation. Even for a single pendulum, the non-linear dynamics means there isn’t any analytical solution for solving for theta, thus the necessity for a numerical solution. The combination is sort of easy (but powerful), we use angular acceleration to update angular velocity and angular velocity to update theta by adding the previous quantity to the latter and multiplying this by a while step. This gets us an approximation for the world 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 now simulated a full pendulum, but this remains to be a well defined system. It’s now time to step into the chaos of the double pendulum.
Chaos, within the mathematical sense, refers to systems which can be highly sensitive to their initial conditions. Even slight changes within the system’s start will result in vastly different behaviors because the system evolves. This perfectly describes the motion of the double pendulum. Unlike the one pendulum, it shouldn’t be a well behaved system and can evolve in a vastly different way with even slight changes in starting angle.
To model the motion of the double pendulum, we are going to use the identical Lagrangian approach as before (see full derivation).
We can even be using the identical 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 now successfully modeled a double pendulum, but now it’s time to watch some chaos. Our final simulation shall be of two double pendulums with barely different starting condition. We are going to set one pendulum to have a theta 1 of 90 degrees and the opposite to have a theta 1 of 91 degrees. Let’s see what happens.
We are able to see that each pendulums start off with similar trajectories but quickly diverge. That is what we mean after we say chaos, even a 1 degree difference in angle cascades into vastly different end behavior.
In this text we learned about pendulum motion and model it. We began from the only harmonic motion model and built as much as the complex and chaotic double pendulum. Along the way in which we learned concerning the Lagrangian, chaos, and numerical integration.
The double pendulum is the only example of a chaotic system. These systems exist in every single place in our world from population dynamics, climate, and even billiards. We are able to take the teachings we’ve learned from the double pendulum and apply them each time we encounter a chaotic systems.
Key Take Aways
- Chaotic systems are very sensitive to initial conditions and can evolve in vastly alternative ways with even slight changes to their start.
- When coping with a system, especially a chaotic system, is there one other frame of reference to have a look at it that makes it easier to work with? (Just like the force reference frame to the energy reference frame)
- When systems get too complicated we want to implement numerical solutions to unravel them. These solutions are easy but powerful and offer good approximations to the actual behavior.
All figures utilized in this text were either created by the creator 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
Easy 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 ensure 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 easy 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