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’re all aware of. Be it a grandfather clock or a toddler on a swing, we’ve seen the common, periodic movement of the pendulum. A single pendulum is properly outlined in classical physics, however the double pendulum (a pendulum connected to the tip of one other pendulum) is literal chaos. On this article, we’re going to construct on our intuitive understanding of pendulums and mannequin the chaos of the double pendulum. The physics is attention-grabbing and the numerical strategies wanted are an important instrument in anybody’s arsenal.

On this article we are going to:

• Find out about harmonic movement and mannequin the conduct of a single pendulum
• Study the basics of chaos principle
• Mannequin the chaotic conduct of a double pendulum numerically

Easy Harmonic Movement

We describe the periodic oscillating motion of a pendulum as harmonic motion. Harmonic movement happens when there may be motion in a system that’s balanced out by a proportional restoring pressure in the other way of mentioned motion. We see an instance of this in determine 2 the place a mass on a spring is being pulled down attributable to gravity, however this places power into the spring which then recoils and pulls the mass again up. Subsequent to the spring system, we see the peak of the mass going round in a circle known as a phasor diagram which additional illustrates the common movement of the system.

Harmonic movement may be damped (lowering in amplitude attributable to drag forces) or pushed (growing in amplitude attributable to exterior pressure being added), however we are going to begin with the best case of indefinite harmonic movement with no exterior forces appearing on it (undamped movement). That is sort of movement is an effective approximation for modeling a single pendulum that swings at a small angle/low amplitude. On this case we will mannequin the movement with equation 1 under.

We are able to simply put this operate into code and simulate a 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) #levels 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 break up into 300 time intervals
theta = []
for t in time_span:
theta.append(simple_pendulum(theta_0, omega, t, phi))
#Convert again to cartesian coordinates
x = l*np.sin(theta)
y = -l*np.cos(theta) #detrimental to verify the pendulum is going through down

Full Pendulum Movement with Lagrangian Mechanics

A easy small angle pendulum is an effective begin, however we wish to transcend this and mannequin the movement of a full pendulum. Since we will now not use small angle approximations it’s best to mannequin the pendulum utilizing Lagrangian mechanics. That is an important instrument in physics that switches us from trying on the forces in a system to trying on the power in a system. We’re switching our body of reference from driving pressure vs restoring pressure to kinetic vs potential power.

The Lagrangain is the distinction between kinetic and potential power 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 power of our system. There may be one final math step to undergo to remodel this into one thing that we will construct a simulation on. We have to bridge again to the dynamic/pressure oriented reference from the power reference utilizing the Euler-Lagrange equation. Utilizing this equation we will use the Lagrangian to get the angular acceleration of our pendulum.

After going via the mathematics, we’ve angular acceleration which we will use to get angular velocity and angle itself. It will 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 isn’t any analytical answer for fixing for theta, thus the necessity for a numerical answer. The mixing is sort of easy (however highly effective), we use angular acceleration to replace angular velocity and angular velocity to replace theta by including the previous amount to the latter and multiplying this by a while step. This will get us an approximation for the world beneath the acceleration/velocity curve. The smaller the time step, the extra correct 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 #Replace velocity with acceleration
theta += time_step*theta_velocity #Replace angle with angular velocity
return theta, theta_velocity
g = 9.8 #m/s^2
l = 1.0 #m
theta = [np.radians(90)] #theta_0
theta_velocity = 0 #Begin with 0 velocity
time_step = 20/300 #Outline a time step
time_span = np.linspace(0,20,300) #simulate for 20s break up 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 again to cartesian coordinates 
x = l*np.sin(theta)
y = -l*np.cos(theta)

We’ve simulated a full pendulum, however that is nonetheless a properly outlined system. It’s now time to step into the chaos of the double pendulum.

Chaos, within the mathematical sense, refers to techniques which can be extremely delicate to their preliminary circumstances. Even slight adjustments within the system’s begin will result in vastly completely different behaviors because the system evolves. This completely describes the movement of the double pendulum. In contrast to the only pendulum, it’s not a properly behaved system and can evolve in a vastly completely different approach with even slight adjustments in beginning angle.

To mannequin the movement of the double pendulum, we are going to use the identical Lagrangian method as earlier than (see full derivation).

We will even be utilizing the identical numerical integration scheme as earlier than when implementing this equation into code and discovering 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)
interplay = -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 + interplay)/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
#Replace 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
#Replace 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

#Outline system parameters
g = 9.8 #m/s^2
m1 = 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 lastly completed it! We’ve efficiently modeled a double pendulum, however now it’s time to watch some chaos. Our ultimate simulation will probably be of two double pendulums with barely completely different beginning situation. We are going to set one pendulum to have a theta 1 of 90 levels and the opposite to have a theta 1 of 91 levels. Let’s see what occurs.

We are able to see that each pendulums begin off with related trajectories however shortly diverge. That is what we imply once we say chaos, even a 1 diploma distinction in angle cascades into vastly completely different finish conduct.

On this article we realized about pendulum movement and easy methods to mannequin it. We began from the best harmonic movement mannequin and constructed as much as the complicated and chaotic double pendulum. Alongside the way in which we realized in regards to the Lagrangian, chaos, and numerical integration.

The double pendulum is the best instance of a chaotic system. These techniques exist in all places in our world from population dynamics, climate, and even billiards. We are able to take the teachings we’ve realized from the double pendulum and apply them every time we encounter a chaotic techniques.

Key Take Aways

• Chaotic techniques are very delicate to preliminary circumstances and can evolve in vastly other ways with even slight adjustments to their begin.
• When coping with a system, particularly a chaotic system, is there one other body of reference to have a look at it that makes it simpler to work with? (Just like the pressure reference body to the power reference body)
• When techniques get too difficult we have to implement numerical options to resolve them. These options are easy however highly effective and supply good approximations to the precise conduct.

All figures used on this article had been both created by the creator or are from Math Images and full beneath 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,title):
!mkdir frames
counter=0
photos = []
for i in vary(0,len(x)):
plt.determine(figsize = (6,6))
plt.plot([0,x[i]],[0,y[i]], "o-", shade = "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")
photos.append(imageio.imread("frames/" + str(counter)+ ".png"))
counter += 1
plt.shut()
imageio.mimsave(title, photos)
!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) #levels 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 break up 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) #detrimental to verify the pendulum is going through 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_velocity
g = 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 break up 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 again to cartesian coordinates 
x = l*np.sin(theta)
y = -l*np.cos(theta)
#Use similar operate 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)
interplay = -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 + interplay)/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

#Outline system parameters, run double pendulum
g = 9.8 #m/s^2
m1 = 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 frames
counter=0
photos = []
for i in vary(0,len(x1)):
plt.determine(figsize = (6,6))
plt.determine(figsize = (6,6))
plt.plot([0,x1[i]],[0,y1[i]], "o-", shade = "b", markersize = 7, linewidth=.7 )
plt.plot([x1[i],x2[i]],[y1[i],y2[i]], "o-", shade = "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")
photos.append(imageio.imread("frames/" + str(counter)+ ".png"))
counter += 1
plt.shut()
imageio.mimsave("double_pendulum.gif", photos)
!rm -r frames