import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
def forwardEuler(f, t_i, y_i, dt):
y_next = y_i + dt*f(t_i, y_i)
return y_next
def rungeKutta4(f, t_i, y_i, dt):
f_1 = f(t_i, y_i)
f_2 = f(t_i + dt/2, y_i + (1/2)*dt*f_1)
f_3 = f(t_i + dt/2, y_i + (1/2)*dt*f_2)
f_4 = f(t_i + dt, y_i + dt*f_3)
y_next = y_i + (1/6)*dt*(f_1 + 2*f_2 + 2*f_3 + f_4)
return y_next
Double Pendulum Equations
'''
Attributes
----------
- ts (list) : saves the time steps
- g (float) : value for gravitational constant in m/s
- L1 (float) : length of the 1st rod in m
- L2 (float) : length of the 2nd rod in m
- b (float) : value for dampening term
- step (float) : value for the step size (in seconds)
- theta_1s (list) : list of angular positions for 1st pendulum (angles are in radians)
- angVel_1s (list) : list of angular velocities for 1st pendulum (in radians per second)
- theta_2s (list) : list of angular positions for 2nd pendulum (angles are in radians)
- angVel_2s (list) : list of angular velocities for 2nd pendulum (in radians per second)
- pos1 (list) : x and y coordinates of 1st object (x1, y1)
- pos2 (list) : x and y coordinates of 2nd object (x2, y2)
Methods
-------
- update : uses RK4 to update the variables
- derivatives : returns the derivatives of the system variables
- displayGraphs : display graphs of time vs position, time vs velocity, and position vs velocity
*The animation function was provided in another cell
'''
class DoublePendulum:
'''
This function initializes the class. The attributes are described above and
the following are not optional: theta1, theta2, angVel1, angVel2, stepSize.
'''
def __init__(self, theta1, theta2, angVel1, angVel2, stepSize, L1 = 1, L2 = 1,
m1 = 2, m2 = 2, dampening = 0):
self.theta_1s = [theta1]
self.angVel_1s = [angVel1]
self.theta_2s = [theta2]
self.angVel_2s = [angVel2]
self.step = stepSize
self.L1 = L1
self.L2 = L2
self.m1 = m1
self.m2 = m2
self.ts = [0]
self.b = dampening
self.g = 9.81
x1 = self.L1*np.sin(theta1)
y1 = -self.L1*np.cos(theta1)
x2 = x1 + self.L2*np.sin(theta2)
y2 = y1 -self.L2*np.cos(theta2)
self.pos1 = [(x1,y1)]
self.pos2 = [(x2,y2)]
self.energy = []
def derivative(self, time, state):
'''
return the derivatives of our state variables
'''
theta1, theta2, angVel1, angVel2 = state
m1, m2, L1, L2, g = self.m1, self.m2, self.L1, self.L2, self.g
theta1_prime = angVel1
theta2_prime = angVel2
p1 = -g*(2*m1 + m2)*np.sin(theta1) - m2*g*np.sin(theta1 - 2*theta2) - \
2*np.sin(theta1 - theta2)*m2*((angVel2**2)*L2 + \
(angVel1**2)*L1*np.cos(theta1 - theta2))
p2 = 2*np.sin(theta1 - theta2) * ((angVel1**2)*L1*(m1 + m2) + \
g*(m1 + m2)*np.cos(theta1) + \
angVel2**2 * L2 * m2 * np.cos(theta1 - theta2))
q1 = L1*(2*m1 + m2 - m2*np.cos(2*theta1 - 2*theta2))
q2 = L2 * (2*m1 + m2 - m2 * np.cos(2 * theta1 - 2 * theta2))
angVel1_prime = p1 / q1
angVel2_prime = p2 / q2
return np.array((theta1_prime, theta2_prime, angVel1_prime, angVel2_prime))
def update(self):
'''
Updates the variables using RK4. Remember not only to add the new position and
velocity to their respective lists but also to update the parameter self.pos in
case you want to use the animation provided
'''
state = (self.theta_1s[-1], self.theta_2s[-1], self.angVel_1s[-1], self.angVel_2s[-1])
next_state = rungeKutta4(self.derivative, self.ts[-1], state, self.step)
theta1, theta2, angVel1, angVel2 = next_state
self.theta_1s.append(theta1)
self.theta_2s.append(theta2)
self.angVel_1s.append(angVel1)
self.angVel_2s.append(angVel2)
self.get_position()
self.get_energy()
t_next = self.ts[-1] + self.step
self.ts.append(t_next)
def get_position(self):
theta1, theta2 = self.theta_1s[-1], self.theta_2s[-1]
x1 = self.L1*np.sin(theta1)
y1 = -self.L1*np.cos(theta1)
x2 = x1 + self.L2*np.sin(theta2)
y2 = y1 -self.L2*np.cos(theta2)
self.pos1.append((x1,y1))
self.pos2.append((x2, y2))
def get_energy(self):
self.P1 = self.m1 * self.pos1[-1][1] * self.g # p1 = m1 g y1
self.P2 = self.m2 * self.pos2[-1][1] * self.g # p2 = m2 g y2
# KE = 0.5 * m * (L omega)**2
self.K1 = 0.5 * self.m1 * (self.L1 * self.angVel_1s[-1])**2
self.K2 = 0.5 * self.m2 * (self.L2 * self.angVel_2s[-1])**2
total_energy = self.P1 + self.P2 + self.K1 + self.K2
self.energy.append(total_energy)
def displayGraphs(self):
'''
This function should plot the following graphs:
- x and y coordinate for both objects
- angular position vs time for both objects
- angular velocity vs time for both objects
'''
%matplotlib inline
plt.figure(figsize=(20,6))
plt.subplot(131)
pos1, pos2 = np.array(self.pos1), np.array(self.pos2)
plt.plot(pos1[:, 0], pos1[:, 1], label = 'Bob 1')
plt.plot(pos2[:, 0], pos2[:, 1], label = 'Bob 2')
plt.legend()
plt.xlabel('x-coordinate')
plt.ylabel('y-coordinate')
plt.subplot(132)
plt.plot(self.ts, self.theta_1s, label = '$\\theta_1$')
plt.plot(self.ts, self.theta_2s, label = '$\\theta_2$')
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Angular Position (rad)')
plt.subplot(133)
plt.plot(self.ts, self.angVel_1s, label = '$\omega_1$')
plt.plot(self.ts, self.angVel_2s, label = '$\omega_2$')
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Angular Velocity (rad/s)')
plt.show()
def display(self, animation = False):
lines = [None, None, None, None]
fig = plt.figure(figsize =(8,8))
L = self.L1 + self.L2
ax = plt.subplot(xlim = (-L*1.5, L*1.5), ylim = (-L*1.5, L*1.5))
# initialize lines to certain values
lines[0], = ax.plot([self.pos1[-1][0]], [self.pos1[-1][1]], ms = 20, color = 'red', marker = 'o')
lines[1], = ax.plot([0, self.pos1[-1][0]], [0, self.pos1[-1][1]], color = 'blue')
lines[2], = ax.plot([self.pos2[-1][0]], [self.pos2[-1][1]], ms = 20, color = 'red', marker = 'o')
lines[3], = ax.plot([self.pos1[-1][0], self.pos2[-1][0]], [self.pos1[-1][1], self.pos2[-1][1]], color = 'green')
ax.set_xticks([])
ax.set_yticks([])
ax.hlines(0, -2, 2, ls = 'dotted')
plt.show()
Animation Code
from IPython.display import HTML
theta1 = np.pi/2
theta2 = np.pi/2
angVel1 = 0
angVel2 = 0
k = 0.02
b = 0
A = 0
omega = 0
L1 = 1
L2 = 1
m1 = 2
m2 = 2
sys = DoublePendulum(theta1, theta2, angVel1, angVel2, k, L1, L2, m1, m2)
#############
# ANIMATION #
#############
#%matplotlib notebook
# create figure, axis, and lines for drawing
lines = [None, None, None, None]
fig = plt.figure(figsize =(8,8))
L = L1 + L2
ax = plt.subplot(xlim = (-L*1.2, L*1.2), ylim = (-L*1.2, L*1.2))
# initialize lines to certain values
lines[0], = ax.plot([sys.pos1[-1][0]], [sys.pos1[-1][1]], ms = 20, color = 'red', marker = 'o')
lines[1], = ax.plot([0, sys.pos1[-1][0]], [0, sys.pos1[-1][1]], color = 'blue')
lines[2], = ax.plot([sys.pos2[-1][0]], [sys.pos2[-1][1]], ms = 20, color = 'red', marker = 'o')
lines[3], = ax.plot([sys.pos1[-1][0], sys.pos2[-1][0]], [sys.pos1[-1][1], sys.pos2[-1][1]], color = 'green')
#ax.set_xticks([])
#ax.set_yticks([])
ax.hlines(0, -2, 2, ls = 'dotted')
def animate(i):
sys.update()
lines[0].set_data(sys.pos1[-1][0], sys.pos1[-1][1])
lines[1].set_data([0, sys.pos1[-1][0]], [0, sys.pos1[-1][1]])
lines[2].set_data(sys.pos2[-1][0], sys.pos2[-1][1])
lines[3].set_data([sys.pos1[-1][0], sys.pos2[-1][0]], [sys.pos1[-1][1], sys.pos2[-1][1]])
return lines
# starts animation
ani = animation.FuncAnimation(fig, animate, np.arange(1, 500), interval = 20, blit = True)
HTML(ani.to_html5_video())
sys.displayGraphs()
sys1 = DoublePendulum(theta1, theta2, angVel1, angVel2, k, L1, L2, m1, m2)
sys2 = DoublePendulum(theta1+0.05, theta2+0.05, angVel1, angVel2, k, L1, L2, m1, m2)
sys3 = DoublePendulum(theta1-0.05, theta2-0.05, angVel1, angVel2, k, L1, L2, m1, m2)
for i in range(2000):
sys1.update()
sys2.update()
sys3.update()
sys1.displayGraphs()
sys2.displayGraphs()
sys3.displayGraphs()
sys1 = DoublePendulum(theta1, theta2, angVel1, angVel2, k, L1, L2, m1, m2)
sys2 = DoublePendulum(theta1+0.005, theta2+0.005, angVel1, angVel2, k, L1, L2, m1, m2)
sys3 = DoublePendulum(theta1-0.005, theta2-0.005, angVel1, angVel2, k, L1, L2, m1, m2)
for i in range(901):
sys1.update()
sys2.update()
sys3.update()
if i % 300 == 0 and i != 0:
plt.figure(figsize=(8,8))
sys1_pos2 = np.array(sys1.pos2)
plt.plot(sys1_pos2[:, 0], sys1_pos2[:, 1], label = 'System 1')
sys2_pos2 = np.array(sys2.pos2)
plt.plot(sys2_pos2[:, 0], sys2_pos2[:, 1], label = 'System 2')
sys3_pos2 = np.array(sys3.pos2)
plt.plot(sys3_pos2[:, 0], sys3_pos2[:, 1], label = 'System 3')
plt.legend()
plt.xlabel('x-coordinate')
plt.ylabel('y-coordinate')
plt.title('Stepsize = '+str(i)+', Time = '+str(np.round(sys1.ts[-1], 2))+' s')
plt.show()
theta1 = np.pi/4
theta2 = np.pi/4
angVel1 = 0
angVel2 = 0
k = 0.01
sys = DoublePendulum(theta1, theta2, angVel1, angVel2, k)
sys.display()
for i in range(2000):
sys.update()
if i%200 == 0:
sys.display()
sys.displayGraphs()
plt.plot(sys.ts[1:], sys.energy)
[<matplotlib.lines.Line2D at 0x7fc1f88cd150>]
from IPython.display import HTML
theta1 = np.pi/2
theta2 = np.pi/2
angVel1 = 0
angVel2 = 0
k = 0.02
sys1 = DoublePendulum(theta1, theta2, angVel1, angVel2, k)
sys2 = DoublePendulum(theta1+0.005, theta2+0.005, angVel1, angVel2, k)
sys3 = DoublePendulum(theta1-0.005, theta2-0.005, angVel1, angVel2, k)
#############
# ANIMATION #
#############
#%matplotlib notebook
# create figure, axis, and lines for drawing
lines = np.full(12, None)
fig = plt.figure(figsize =(8,8))
L = L1 + L2
ax = plt.subplot(xlim = (-L*1.2, L*1.2), ylim = (-L*1.2, L*1.2))
# initialize lines to certain values
lines[0], = ax.plot([sys1.pos1[-1][0]], [sys1.pos1[-1][1]], ms = 4, color = 'red', marker = 'o')
lines[1], = ax.plot([0, sys1.pos1[-1][0]], [0, sys1.pos1[-1][1]], color = 'blue')
lines[2], = ax.plot([sys1.pos2[-1][0]], [sys1.pos2[-1][1]], ms = 4, color = 'blue', marker = 'o')
lines[3], = ax.plot([sys1.pos1[-1][0], sys1.pos2[-1][0]], [sys1.pos1[-1][1], sys1.pos2[-1][1]], color = 'blue')
lines[4], = ax.plot([sys2.pos1[-1][0]], [sys2.pos1[-1][1]], ms = 4, color = 'red', marker = 'o')
lines[5], = ax.plot([0, sys2.pos1[-1][0]], [0, sys2.pos1[-1][1]], color = 'green')
lines[6], = ax.plot([sys2.pos2[-1][0]], [sys2.pos2[-1][1]], ms = 4, color = 'green', marker = 'o')
lines[7], = ax.plot([sys2.pos1[-1][0], sys2.pos2[-1][0]], [sys2.pos1[-1][1], sys2.pos2[-1][1]], color = 'green')
lines[8], = ax.plot([sys3.pos1[-1][0]], [sys3.pos1[-1][1]], ms = 4, color = 'red', marker = 'o')
lines[9], = ax.plot([0, sys3.pos1[-1][0]], [0, sys3.pos1[-1][1]], color = 'red')
lines[10], = ax.plot([sys3.pos2[-1][0]], [sys3.pos2[-1][1]], ms = 4, color = 'red', marker = 'o')
lines[11], = ax.plot([sys3.pos1[-1][0], sys3.pos2[-1][0]], [sys3.pos1[-1][1], sys3.pos2[-1][1]], color = 'red')
#ax.set_xticks([])
#ax.set_yticks([])
ax.hlines(0, -2, 2, ls = 'dotted')
ax.grid(True, alpha = 0.75)
def animate(i):
sys1.update()
sys2.update()
sys3.update()
lines[0].set_data(sys1.pos1[-1][0], sys1.pos1[-1][1])
lines[1].set_data([0, sys1.pos1[-1][0]], [0, sys1.pos1[-1][1]])
lines[2], = ax.plot([sys1.pos2[-1][0]], [sys1.pos2[-1][1]], ms = 2, color = 'blue', marker = 'o')
lines[3].set_data([sys1.pos1[-1][0], sys1.pos2[-1][0]], [sys1.pos1[-1][1], sys1.pos2[-1][1]])
lines[4].set_data(sys2.pos1[-1][0], sys2.pos1[-1][1])
lines[5].set_data([0, sys2.pos1[-1][0]], [0, sys2.pos1[-1][1]])
lines[6], = ax.plot([sys2.pos2[-1][0]], [sys2.pos2[-1][1]], ms = 2, color = 'green', marker = 'o')
lines[7].set_data([sys2.pos1[-1][0], sys2.pos2[-1][0]], [sys2.pos1[-1][1], sys2.pos2[-1][1]])
lines[8].set_data(sys3.pos1[-1][0], sys3.pos1[-1][1])
lines[9].set_data([0, sys3.pos1[-1][0]], [0, sys3.pos1[-1][1]])
lines[10], = ax.plot([sys3.pos2[-1][0]], [sys3.pos2[-1][1]], ms = 2, color = 'red', marker = 'o')
lines[11].set_data([sys3.pos1[-1][0], sys3.pos2[-1][0]], [sys3.pos1[-1][1], sys3.pos2[-1][1]])
return lines
# starts animation
ani = animation.FuncAnimation(fig, animate, np.arange(1, 1000), interval = 20, blit = True)
HTML(ani.to_html5_video())
https://www.myphysicslab.com/pendulum/double-pendulum-en.html
Here, they used, $$ L_1 = L_2 = 1\\ m_1 = m_2 = 2\\ \theta_1 = 0.4\\ \theta_2 = 0\\ \omega_1 = \omega_2 = 0 $$
This is the expected graph of theta1 (green) and theta2 (red) over time:
plt.figure(figsize=(13,6))
plt.subplot(121)
plt.plot(sys.theta_1s[:328], sys.angVel_1s[:328])
plt.subplot(122)
plt.plot(sys.theta_2s[:328], sys.angVel_2s[:328])
[<matplotlib.lines.Line2D at 0x7f8386ee3c90>]
plt.plot(sys.theta_1s, sys.theta_2s)
[<matplotlib.lines.Line2D at 0x7f8386dc1550>]
from IPython.display import HTML
theta = np.pi/2
angVel = 0
k = 0.01
L = 10
b = 1.6
sys = SimplePendulum(theta, angVel, k, L, b)
#############
# ANIMATION #
#############
#%matplotlib notebook
# create figure, axis, and lines for drawing
lines = [None, None]
fig = plt.figure(figsize =(3,3))
ax = plt.subplot(xlim = (-sys.L*1.5, sys.L*1.5), ylim = (-sys.L*1.5, sys.L*1.5))
# initialize lines to certain values
lines[0], = ax.plot([sys.pos[0]], [sys.pos[1]], ms = 20, color = 'red', marker = 'o')
lines[1], = ax.plot([0, sys.pos[0]], [0, sys.pos[1]], color = 'blue')
ax.set_xticks([])
ax.set_yticks([])
ax.hlines(0, -2, 2, ls = 'dotted')
def animate(i):
sys.update()
lines[0].set_data(sys.pos[0], sys.pos[1])
lines[1].set_data([0, sys.pos[0]], [0, sys.pos[1]])
return lines
# starts animation
ani = animation.FuncAnimation(fig, animate, np.arange(1, 2000), interval = 10, blit = True)
HTML(ani.to_html5_video())
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-59-e340f77d5d04> in <module>() 6 L = 10 7 b = 1.6 ----> 8 sys = SimplePendulum(theta, angVel, k, L, b) 9 10 ############# NameError: name 'SimplePendulum' is not defined
sys.displayGraphs()