!apt update -y
!apt install ffmpeg -y
Hit:1 http://deb.debian.org/debian buster InRelease
Hit:2 http://deb.debian.org/debian buster-updates InRelease
Hit:3 http://security.debian.org/debian-security buster/updates InRelease
2 packages can be upgraded. Run 'apt list --upgradable' to see them.
ffmpeg is already the newest version (7:4.1.9-0+deb10u1).
0 upgraded, 0 newly installed, 0 to remove and 2 not upgraded.
# Import Libraries:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FFMpegWriter
from matplotlib.patches import Circle
from matplotlib.patches import Rectangle
# Parameters
g = 9.8
k = 300
a = 0.5
L1 = 1
L2 = 0.5
m1 = 150
m2 = 50
b = a/0.8
# We need an array of time points from 0 to 10 in increments of 0.01 seconds
# Simulation Parameters
dt = 0.001
sim_time = 10
t_vec = np.arange(0,sim_time,dt)
vec_size = len(t_vec)
F = np.arange(0,sim_time,dt)
#F = 0.10*np.sin(t_vec)
for i in range(0,vec_size):
F[i] = 0.10*np.sin(t_vec[i])
#print(F)
# Initialize a vector of zeros
th1_vec = np.zeros(len(t_vec))
dth1_vec = np.zeros(len(t_vec))
th2_vec = np.zeros(len(t_vec))
dth2_vec = np.zeros(len(t_vec))
# Set our initial condition
th1_vec[0] = -np.pi/2 # Initial th1
dth1_vec[0] = 0 # Initial th1 velocity
th2_vec[0] = 0 # Initial th2
dth2_vec[0] = 0 # Initial th2 velocity
# Loop through time
# Euler's Method (approximately integrates the differential equation)
for i in range(1, len(t_vec)):
th1_vec[i] = th1_vec[i-1] + dth1_vec[i-1]*dt
dth1_vec[i] = dth1_vec[i-1] + (-(g/L1)*np.sin(th1_vec[i-1]) + ((k*a**2)/(m1*L1**2))*(th2_vec[i-1]-th1_vec[i-1]) + (b/m1)*np.sin(th1_vec[i-1])*np.cos(th1_vec[i-1]))*dt
th2_vec[i] = th2_vec[i-1] + dth2_vec[i-1]*dt
dth2_vec[i] = dth2_vec[i-1] + (-(g/L2)*np.sin(th2_vec[i-1]) - ((k*a**2)/(m2*L2**2))*(th2_vec[i-1]-th1_vec[i-1]) + (b/m2)*np.sin(th2_vec[i-1])*np.cos(th2_vec[i-1]))*dt
th1_vec = th1_vec[::-1]
dth1_vec = dth1_vec[::-1]
th2_vec = th2_vec[::-1]
dth2_vec = dth2_vec[::-1]
plt.plot(t_vec,th1_vec)
plt.title("th1 vs Time")
plt.show()
plt.plot(t_vec,dth1_vec)
plt.title("dth1 vs Time")
plt.show()
plt.plot(t_vec,th2_vec)
plt.title("th2 vs Time")
plt.show()
plt.plot(t_vec,dth2_vec)
plt.title("dth vs Time")
plt.show()
class Pendulum_Calculations:
def __init__(self): # Constructor
self.x_pendulum_arm = None
self.y_pendulum_arm = None
self.origin = None
self.x_data_points = None
self.y_data_points = None
self.patch_center = None
# Making each link in system an Object
Pend1 = Pendulum_Calculations()
Pend2 = Pendulum_Calculations()
Spring = Pendulum_Calculations()
# Create animation
# Setup figure:
fig, ax = plt.subplots()
p1, = ax.plot([],[],color='dimgray')
p2, = ax.plot([],[],color='dimgray')
spring, = ax.plot([],[],color='red')
min_lim = -3
max_lim = 3
ax.axis('equal')
ax.set_xlim([min_lim, max_lim])
ax.set_ylim([min_lim, max_lim])
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_title('Double Pendulum Spring System')
title = "Double Pendulum Spring System"
# Setup animation writer
FPS = 20
sample_rate = int(1/(dt*FPS)) # Real time playback
dpi = 300
writerObj = FFMpegWriter(fps=FPS)
# We need to calculate the cartesian location of the pendulum:
# Initialize Array:
simulation_size = len(t_vec) # We use this alot, so lets store it.
Pend1.origin = (-1/2, 0) # Pin Joint (X, Y) Location of pendulum 1
Pend2.origin = ( 1/2, 0) # Pin Joint (X, Y) Location of pendulum 2
Pend1.x_pendulum_arm = np.zeros(simulation_size)
Pend1.y_pendulum_arm = np.zeros(simulation_size)
Pend2.x_pendulum_arm = np.zeros(simulation_size)
Pend2.y_pendulum_arm = np.zeros(simulation_size)
Spring.x_pendulum_arm = np.zeros(simulation_size)
Spring.y_pendulum_arm = np.zeros(simulation_size)
# Find the X and Y trajectories of the pendulums:
for i in range(0, simulation_size):
Pend1.x_pendulum_arm[i] = L1 * np.sin(th1_vec[i]) + Pend1.origin[0]
Pend1.y_pendulum_arm[i] = Pend1.origin[1] - L1 * np.cos(th1_vec[i])
Pend2.x_pendulum_arm[i] = L2 * np.sin(th2_vec[i]) + Pend2.origin[0]
Pend2.y_pendulum_arm[i] = Pend2.origin[1] - L2 * np.cos(th2_vec[i])
Spring.x_pendulum_arm[i] = L1/2 * np.sin(th1_vec[i]) + Pend1.origin[0]
Spring.y_pendulum_arm[i] = Pend1.y_pendulum_arm[i]/2
# Initialize Patch: (Pendulum 1)
c1 = Circle(Pend1.origin, radius=0.1, color='cornflowerblue') # Draws a circle of radius 0.1 at (X, Y) location (0, 0)
ax.add_patch(c1) # Add the patch to the axes
hinge1 = Circle((Pend1.origin), radius=0.05, color='black')
ax.add_patch(hinge1)
# Initialize Patch: (Pendulum 2)
c2 = Circle(Pend2.origin, radius=0.1, color='cornflowerblue') # Draws a circle of radius 0.1 at (X, Y) location (0, 0)
ax.add_patch(c2) # Add the patch to the axes
hinge2 = Circle((Pend2.origin), radius=0.05, color='black')
ax.add_patch(hinge2)
# Plot and Create Animation
with writerObj.saving(fig, title+".mp4", dpi):
for i in range(0, vec_size, sample_rate):
# Update Pendulum Arm 1:
Pend1.x_data_points = [Pend1.origin[0], Pend1.x_pendulum_arm[i]]
Pend1.y_data_points = [Pend1.origin[1], Pend1.y_pendulum_arm[i]]
p1.set_data(Pend1.x_data_points, Pend1.y_data_points)
# Update Pendulum Patch:
Pend1.patch_center = Pend1.x_pendulum_arm[i], Pend1.y_pendulum_arm[i]
c1.center = Pend1.patch_center
# Update Pendulum Arm 2:
Pend2.x_data_points = [Pend2.origin[0], Pend2.x_pendulum_arm[i]]
Pend2.y_data_points = [Pend2.origin[1], Pend2.y_pendulum_arm[i]]
p2.set_data(Pend2.x_data_points, Pend2.y_data_points)
# Update Pendulum Patch:
Pend2.patch_center = Pend2.x_pendulum_arm[i], Pend2.y_pendulum_arm[i]
c2.center = Pend2.patch_center
# Update Spring:
Spring.x_data_points = [Spring.x_pendulum_arm[i], Pend2.x_pendulum_arm[i]]
Spring.y_data_points = [Spring.y_pendulum_arm[i], Pend2.y_pendulum_arm[i]]
spring.set_data(Spring.x_data_points, Spring.y_data_points)
# Update drawing, grab and save frame
fig.canvas.draw()
writerObj.grab_frame()
from IPython.display import Video
Video("/work/Double Pendulum Spring System.mp4", embed=True, width=640, height=480)