[ad_1]
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.
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
Source link
