# MOD HW 1 # Question 2 # Tessany Schou # Install ffmpeg (allows us to create videos) !apt update -y !apt install ffmpeg -y # Normal Python Mathy Code # Import Libraries: import numpy as np import math import matplotlib.pyplot as plt from matplotlib.animation import FFMpegWriter from matplotlib.patches import Circle from matplotlib.patches import Rectangle # Parameters m1 = 1 # kg m2 = 2 # kg g = 9.8 # m/s^2 k = 10 # N/m h = 1.5 # m L = 1 # m # We need an array of time points from 0 to 10 in increments of 0.001 seconds dt = 0.001 t_vec = np.arange(0,10,dt) sim_length = len(t_vec) # Initialize a vector of zeros x_vec = np.zeros(len(t_vec)) dx_vec = np.zeros(len(t_vec)) th_vec = np.zeros(len(t_vec)) dth_vec = np.zeros(len(t_vec)) # Set our initial condition x_vec[0] = 0 # initial x position (m) dx_vec[0] = 0 # initial x velocity (m/s) th_vec[0] = 1 # initial theta position (rad) dth_vec[0] = 0 # initial theta velocity (rad/s) # Loop through time # Euler's Method (approximately integrates the differential equation) for i in range(1, len(t_vec)): x_vec[i] = x_vec[i-1] + dx_vec[i-1]*dt dx_vec[i] = dx_vec[i-1] + ((k*(L - ((h - L*math.cos(th_vec[i-1]))**2 + (x_vec[i-1] - L*math.sin(th_vec[i-1]))**2)**(1/2))*(2*x_vec[i-1] - 2*L*math.sin(th_vec[i-1])))/(2*((h - L*math.cos(th_vec[i-1]))**2 + (x_vec[i-1] - L*math.sin(th_vec[i-1]))**2)**(1/2)))/m2*dt th_vec[i] = th_vec[i-1] + dth_vec[i-1]*dt dth_vec[i] = dth_vec[i-1] + (- (k*(2*L*math.cos(th_vec[i-1])*(x_vec[i-1] - L*math.sin(th_vec[i-1])) - 2*L*math.sin(th_vec[i-1])*(h - L*math.cos(th_vec[i-1])))*(L - ((h - L*math.cos(th_vec[i-1]))**2 + (x_vec[i-1] - L*math.sin(th_vec[i-1]))**2)**(1/2)))/(2*((h - L*math.cos(th_vec[i-1]))**2 + (x_vec[i-1] - L*math.sin(th_vec[i-1]))**2)**(1/2)) - L*g*h*m1*math.sin(th_vec[i-1]))/(L**2*m1)*dt # Plotting x positon plt.xlabel('t') plt.ylabel('x') plt.title('x(t)') plt.plot(t_vec,x_vec) plt.show() # Plotting y position plt.xlabel('t') plt.ylabel('θ') plt.title('θ(t)') plt.plot(t_vec,th_vec) plt.show() # Setup Figure: Initialize Figure / Axe Handles fig, ax = plt.subplots() p1, = ax.plot([], [], color='cornflowerblue') p2, = ax.plot([], [], color='cornflowerblue') ax.axis('equal') ax.set_xlim([-2.5, 2.5]) ax.set_ylim([-0.5, 2]) ax.set_xlabel('X') ax.set_ylabel('Y') ax.set_title('Question 2 Simulation:') video_title = "Q2simulation" # Now we want to plot stuff on those axes c = Circle((L*math.sin(th_vec[0]),h-L*math.cos(th_vec[0])), radius=0.1, color='cornflowerblue') ax.add_patch(c) # Setup Animation Writer: FPS = 20 sample_rate = int(1 / (dt * FPS)) dpi = 300 writerObj = FFMpegWriter(fps=FPS) # Draw Pin Joint: pin_joint1 = Circle((0, h), radius=0.025, color='black', zorder=10) ax.add_patch(pin_joint1) pin_joint2 = Circle((x_vec[0],0), radius=0.025, color='black', zorder=10) ax.add_patch(pin_joint2) # Plot and Create Animation: with writerObj.saving(fig, video_title+".mp4", dpi): for i in range(0, sim_length, sample_rate): # Update Pendulum Arm: xp_data_points = [0, L*math.sin(th_vec[i])] yp_data_points = [h, h-L*math.cos(th_vec[i])] xb_data_points = [L*math.sin(th_vec[i]), x_vec[i]] yb_data_points = [h-L*math.cos(th_vec[i]), 0] p1.set_data(xp_data_points, yp_data_points) p2.set_data(xb_data_points, yb_data_points) # Update Pendulum Patch: pendulum_center = L*math.sin(th_vec[i]), h-L*math.cos(th_vec[i]) c.center = pendulum_center block_center = x_vec[i],0 pin_joint2.center = block_center # Update Drawing: fig.canvas.draw() # Update the figure with the new changes # Grab and Save Frame: writerObj.grab_frame() from IPython.display import Video Video("/work/"+video_title+".mp4", embed=True, width=640, height=480)