#Parabolic Constrained Harmonic Motion #install ffmpeg !apt update -y !apt install ffmpeg -y # Import Libraries: import numpy as np import matplotlib.pyplot as plt import sympy as sp from matplotlib.animation import FFMpegWriter from matplotlib.patches import Circle # Pendulum Parameters g = -9.81 m1 = 1 m2 = 2 k = 10 l1 = 1 lk = 1.5 # We need an array of time points from 0 to 10 in increments of 0.01 seconds dt = 0.001 t_vec = np.arange(0,10,dt) # Initialize a vector of zeros x_vec = np.zeros(len(t_vec)) dx_vec = np.zeros(len(t_vec)) ddx_vec = np.zeros(len(t_vec)) theta_vec = np.zeros(len(t_vec)) dtheta_vec = np.zeros(len(t_vec)) ddtheta_vec = np.zeros(len(t_vec)) y_loc = np.zeros(len(t_vec)) # Set our initial condition x_vec[0] = 0 # initial X location dx_vec[0] = 0 # initial horizontal velocity theta_vec[0] = 1 # radian dtheta_vec[0] = 0 # radians per second # Loop through time # Euler's Method (approximately integrates the differential equation (Equations Solved via Matlab)) for i in range(1, len(t_vec)): x_vec[i] = x_vec[i-1] + dx_vec[i-1]*dt ddx_vec[i-1] = -(5*(np.sin(theta_vec[i-1]) - x_vec[i-1])*(2*(x_vec[i-1]**2 - 3*np.cos(theta_vec[i-1]) - 2*np.sin(theta_vec[i-1])*x_vec[i-1] + np.cos(theta_vec[i-1])**2 + np.sin(theta_vec[i-1])**2 + 9/4)**(1/2) - 3))/(2*(x_vec[i-1]**2 - 3*np.cos(theta_vec[i-1]) - 2*np.sin(theta_vec[i-1])*x_vec[i-1] + np.cos(theta_vec[i-1])**2 + np.sin(theta_vec[i-1])**2 + 9/4)**(1/2)) dx_vec[i] = dx_vec[i-1] + ddx_vec[i-1]*dt theta_vec[i] = theta_vec[i-1] + dtheta_vec[i-1]*dt ddtheta_vec[i-1] = -(2250*np.sin(theta_vec[i-1]) - 519*np.sin(theta_vec[i-1])*(x_vec[i-1]**2 - 3*np.cos(theta_vec[i-1]) - 2*np.sin(theta_vec[i-1])*x_vec[i-1] + np.cos(theta_vec[i-1])**2 + np.sin(theta_vec[i-1])**2 + 9/4)**(1/2) - 1500*np.cos(theta_vec[i-1])*x_vec[i-1] + 1000*np.cos(theta_vec[i-1])*x_vec[i-1]*(x_vec[i-1]**2 - 3*np.cos(theta_vec[i-1]) - 2*np.sin(theta_vec[i-1])*x_vec[i-1] + np.cos(theta_vec[i-1])**2 + np.sin(theta_vec[i-1])**2 + 9/4)**(1/2))/(100*(x_vec[i-1]**2 - 3*np.cos(theta_vec[i-1]) - 2*np.sin(theta_vec[i-1])*x_vec[i-1] + np.cos(theta_vec[i-1])**2 + np.sin(theta_vec[i-1])**2 + 9/4)**(1/2)) dtheta_vec[i] = dtheta_vec[i-1] + ddtheta_vec[i-1]*dt y_loc[i] = -1.5 fig1= plt.plot(t_vec, x_vec) plt.title('X Positon of m2 vs. Time') plt.xlabel('Time') plt.ylabel('X Position of m2 (m)') plt.show() fig2 = plt.plot(t_vec, theta_vec) plt.title('Theta vs. Time') plt.xlabel('Time') plt.ylabel('Theta (rad)') plt.show() # Setup Figure: Initialize Figure / Axe Handles fig, ax = plt.subplots() # Initialize Figure and Axes p, = ax.plot([], [], color='royalblue') # Initialize Empty Plot q, = ax.plot([], [], color='slategray') ax.axis('equal') ax.set_xlim([-3, 3]) # X Lim ax.set_ylim([-3, 3]) # Y Lim ax.set_xlabel('X') # X Label ax.set_ylabel('Y') # Y Label ax.set_title('Mass-Spring-Pendulum Simulation:') # Plot Title video_title = "spmapen" # name of the mp4 # Initialize Patches: c = Circle((1, -0.65), radius=0.1, color='red') # Draws a circle of radius 0.1 at (X, Y) location (0, 0) ax.add_patch(c) # Add the patch to the axes b = Circle((0, -1.5), radius=0.15, color='lime') ax.add_patch(b) # Add the patch to the axes # Setup Animation Writer: FPS = 20 # Frames per second for the video sample_rate = int(1 / (dt * FPS)) # Rate at which we sample the simulation to visualize dpi = 300 # Quality of the Video writerObj = FFMpegWriter(fps=FPS) # Video Object # 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. x_pendulum_arm = np.zeros(simulation_size) y_pendulum_arm = np.zeros(simulation_size) origin = (0,0) # pin joint at 0,0 for i in range(0, simulation_size): x_pendulum_arm[i] = l1 * np.sin(theta_vec[i]) + origin[0] y_pendulum_arm[i] = origin[1] - l1 * np.cos(theta_vec[i]) # Plot and Create Animation: with writerObj.saving(fig, video_title+".mp4", dpi): # We want to create video that represents the simulation # So we need to sample only a few of the frames: for i in range(0, simulation_size, sample_rate): # Update Pendulum Arm: x_data_points = [origin[0], x_pendulum_arm[i]] y_data_points = [origin[1], y_pendulum_arm[i]] kx_data_points = [x_vec[i],x_pendulum_arm[i]] ky_data_points = [y_loc[i],y_pendulum_arm[i]] p.set_data(x_data_points, y_data_points) # Update plot with set_data q.set_data(kx_data_points, ky_data_points) circ_center = x_pendulum_arm[i], y_pendulum_arm[i] c.center = circ_center m_center = x_vec[i], y_loc[i] b.center = m_center # Update Drawing: fig.canvas.draw() # Update the figure with the new changes # Grab and Save Frame: writerObj.grab_frame() # we can import the video and run it on the notebook from IPython.display import Video Video("/work/spmapen.mp4",embed=True, width=640,height=480)