################## Homework 1 Problem 2 ########################### # 992759-MOD-HW1 !apt update -y !apt install ffmpeg -y import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FFMpegWriter from matplotlib.patches import Circle g = 9.8 # m/s^2 L = 1 # m m_1 = 1 # kg m_2 = 2 # kg k = 10 # N/m h = 1.5 # height from the hinge to the ground (m) # define our time vector for simulation dt = 0.001 t_vec = np.arange(0,10,dt) # define vectors that will hold our state variables over time (initializing vectors) 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 = np.zeros(len(t_vec)) # hiinge location # initial conditions x_vec[0] = 0 dx_vec[0] = 0 theta_vec[0] = 1 dtheta_vec[0] = 0 # loop for theta for i in range(1, len(t_vec)): theta_vec[i] = theta_vec[i-1] + dtheta_vec[i-1]*dt ddtheta_vec[i-1] = -((g*np.sin(theta_vec[i-1]))/(L)) - (k*(np.sqrt((L*np.sin(theta_vec[i-1]) - x_vec[i-1])**2 + (h - L*np.cos(theta_vec[i-1]))**2) - L)*((L*np.sin(theta_vec[i-1]) - x_vec[i-1])*(L*np.cos(theta_vec[i-1])) + (h - L*np.cos(theta_vec[i-1]))*(L*np.sin(theta_vec[i-1]))))/(m_1*L**2*np.sqrt((L*np.sin(theta_vec[i-1]) - x_vec[i-1])**2 + (h - L*np.cos(theta_vec[i-1]))**2)) dtheta_vec[i] = dtheta_vec[i-1] + ddtheta_vec[i-1]*dt # x loop for i in range(1, len(t_vec)): x_vec[i] = x_vec[i-1] + dx_vec[i-1]*dt ddx_vec[i-1] = (k*(np.sqrt((L*np.sin(theta_vec[i-1]) - x_vec[i-1])**2 + (h - L*np.cos(theta_vec[i-1]))**2) - L)*(L*np.sin(theta_vec[i-1]) - x_vec[i-1]))/(m_2*np.sqrt((L*np.sin(theta_vec[i-1]) - x_vec[i-1])**2 + (h - L*np.cos(theta_vec[i-1]))**2)) dx_vec[i] = dx_vec[i-1] + ddx_vec[i-1]*dt y[i] = -1.5 # hinge ########### Plotting ############# # x plot plt.plot(t_vec, x_vec) plt.title('X position') plt.show() # y plot plt.plot(t_vec, theta_vec) plt.title('Theta angle position') plt.show() ####### Animation ########## fig, ax = plt.subplots() p1, = ax.plot([], [], color='cornflowerblue') p2, = ax.plot([], [], color='cornflowerblue') ax.axis('equal') ax.set_xlim([-3, 3]) ax.set_ylim([-3, 3]) ax.set_xlabel('X') ax.set_ylabel('Y') ax.set_title('Problem 2 Simulation') video_title = "Problem 2 Simulation" c = Circle((0, -0.75), radius=0.1, color='cornflowerblue') ax.add_patch(c) # more circle b = Circle((0, -1.5), radius=0.2, color='cornflowerblue') ax.add_patch(b) # Setup Animation Writer: FPS = 20 sample_rate = int(1 / (dt * FPS)) dpi = 300 writerObj = FFMpegWriter(fps=FPS) # calculate the pendulum location in cartesian coords simulation_size = len(t_vec) x_link = np.zeros(simulation_size) # pendulum link y_link = np.zeros(simulation_size) origin = (0, 0) # # loop to get the trajectories for i in range(0, simulation_size): x_link[i] = L * np.sin(theta_vec[i]) + origin[0] y_link[i] = origin[1] - L * np.cos(theta_vec[i]) # Plot and Create Animation: with writerObj.saving(fig, video_title+".mp4", dpi): for i in range(0, simulation_size, sample_rate): x_data = [origin[0], x_link[i]] # updating link y_data = [origin[1], y_link[i]] x2_data = [x_vec[i], x_link[i]] y2_data = [y[i], y_link[i]] p1.set_data(x_data, y_data) # updating plot p2.set_data(x2_data, y2_data) patch_center = x_link[i], y_link[i] # patch update c.center = patch_center m_center = x_vec[i], y[i] b.center = m_center fig.canvas.draw() writerObj.grab_frame() from IPython.display import Video Video("/work/Problem 2 Simulation.mp4", embed=True, width=640, height=480)