Modeling a chaotic system
The pendulum is a classical physics system that we’re all acquainted with. Be it a grandfather clock or a toddler on a swing, now we have 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 a necessary tool in anyone’s arsenal.
In this text we’ll:
- Find out about 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’s 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 could be damped (decreasing in amplitude attributable to drag forces) or driven (increasing in amplitude attributable to outside force being added), but we’ll start with the best case of indefinite harmonic motion with no outside forces acting on it (undamped motion). That is sort of motion is approximation for modeling a single pendulum that swings at a small angle/low amplitude. On this case we will model the motion with equation 1 below.
We are able to easily put this function into code and simulate an easy 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 certain the pendulum is facing down
Full Pendulum Motion with Lagrangian Mechanics
A straightforward small angle pendulum is start, but we wish to transcend this and model the motion of a full pendulum. Since we will now not use small angle approximations it’s best to model the pendulum using Lagrangian mechanics. That is a necessary tool in physics that switches us from the forces in a system to 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’s one last math step to undergo to rework this into something that we will construct a simulation on. We’d like to bridge back to the dynamic/force oriented reference from the energy reference using the Euler-Lagrange equation. Using this equation we will use the Lagrangian to get the angular acceleration of our pendulum.
After going through the mathematics, now we have angular acceleration which we will use to get angular velocity and angle itself. This can require some numerical integration that will probably be specified by our full pendulum simulation. Even for a single pendulum, the non-linear dynamics means there is no such thing as a analytical solution for solving for theta, thus the necessity for a numerical solution. The mixing is kind 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 realm 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 now have 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 might 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 only pendulum, it just isn’t 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’ll use the identical Lagrangian approach as before (see full derivation).
We will 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 now have successfully modeled a double pendulum, but now it’s time to look at some chaos. Our final simulation will probably 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 the best way to model it. We began from the best harmonic motion model and built as much as the complex and chaotic double pendulum. Along the way in which we learned in regards to the Lagrangian, chaos, and numerical integration.
The double pendulum is the best example of a chaotic system. These systems exist in all places in our world from population dynamics, climate, and even billiards. We are able to take the teachings now we have learned from the double pendulum and apply them every time we encounter a chaotic systems.
Key Take Aways
- Chaotic systems are very sensitive to initial conditions and can evolve in vastly other 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 take 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 make certain 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