# Homework 1 Problem 2 # This code block installs ffmpeg to your Deepnote Machine # Allows Deepnote to write MP4 !apt update -y !apt install ffmpeg -y # 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 from matplotlib.pyplot import * from numpy import * # Part b: ------------------------------------------------------------ #Simulate this system using Euler’s method (time step = 0.001 # seconds) in Python with initial conditions, x = 0m, dx/dt = 0m/s, # θ = 1 rad, dθ/dt = 0rad/s over the course of 10 seconds. # Parameters g = 9.81 # m/s^2 m1 = 1 # kg L = 1 # m k = 10 # N/m m2 = 2 # kg h = 1.5 # m # Time array dt = 0.001 t_vec = np.arange(0,10,dt) # Initialize a vector of zeros vec_size = len(t_vec) x_vec = np.zeros(vec_size) dx_vec = np.zeros(vec_size) theta_vec = np.zeros (vec_size) dtheta_vec = np.zeros (vec_size) # Initial Conditions x_vec[0] = 0 dx_vec[0] = 0 theta_vec[0] = 1 dtheta_vec[0] = 0 # Euler Simulation: Using Matrix Form (A * x = B) # Initialize A and B: A = np.array([[m2, 0], [0, m1*L**2]]) b = np.array([[0], [0]]) for i in range (1, vec_size): # Replace non constant values b [0] = -k*x_vec[i-1] + k*L*np.sin(theta_vec[i-1]) b [1] = -k*L**2 * np.sin(theta_vec[i-1]) * np.cos(theta_vec[i-1]) + k*x_vec[i-1]*L*np.cos(theta_vec[i-1]) - k*h*L*np.sin(theta_vec[i-1]) + 0.5*k*L**2 * np.sin(2*theta_vec[i-1]) - m1*g*L*np.sin(theta_vec[i-1]) # Solve for ddx, ddtheta [ddx, ddtheta] = np.linalg.solve(A, b) # Euler integration x_vec [i] = x_vec[i-1] + dx_vec[i-1] * dt theta_vec[i] = theta_vec[i - 1] + dtheta_vec[i - 1] * dt dx_vec[i] = dx_vec[i-1] + ddx * dt dtheta_vec[i] = dtheta_vec[i - 1] + ddtheta * dt # x position vs time plotting plt.plot(t_vec, x_vec) plt.title("x position vs time") plt.show() # theta angle vs time plotting plt.plot(t_vec, theta_vec) plt.title("angle theta vs time") plt.show() # Part c: ------------------------------------------------------------ # Animate the motion from Problem 2b (at 20 frames per second). # Setup Figure: fig, ax = plt.subplots() p, = ax.plot([], [], color='royalblue') min_lim = -5 max_lim = 5 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('Pendulum and Block Simulation:') title = "simulation" # Setup Animation Writer: FPS = 20 sample_rate = int(1 / (dt * FPS)) # Real Time Playback dpi = 300 writerObj = FFMpegWriter(fps=FPS) # Initialize Patch: (Mass Block: m2) width = 1 # Width of Cart height = width / 2 # Height of Cart xy_cart = (x_vec[0] - width / 2, - height / 2) # Bottom Left Corner of Cart r = Rectangle(xy_cart, width, height, color='cornflowerblue') # Rectangle Patch ax.add_patch(r) # Add Patch to Plot # Draw the Ground: ground = ax.hlines(-height / 2, min_lim, max_lim, colors='black') height_hatch = 0.25 width_hatch = max_lim - min_lim xy_hatch = (min_lim, -height / 2 - height_hatch) ground_hatch = Rectangle(xy_hatch, width_hatch, height_hatch, facecolor='None', linestyle='None', hatch='/') ax.add_patch(ground_hatch) # Initialize Patch: (Mass Ball: m1) c = Circle((0, 0), radius=0.1, color='cornflowerblue') # Draws a circle of radius 0.1 at (X, Y) location (0, 0) ax.add_patch(c) # Add the patch to the axes # We need to calculate the cartesian location of the pendulum: # Initialize Array: vec_size = len(t_vec) # We use this alot, so lets store it. x_pendulum_arm = np.zeros(vec_size) y_pendulum_arm = np.zeros(vec_size) origin = (0, h) # Pin Joint (X, Y) Location # Find the X and Y trajectories of the pendulum: for i in range(0, vec_size): x_pendulum_arm[i] = L * np.sin(theta_vec[i]) + origin[0] y_pendulum_arm[i] = origin[1] - L * np.cos(theta_vec[i]) # Animate: with writerObj.saving(fig, title + ".mp4", dpi): for i in range(0, vec_size, sample_rate): # Update Pendulum Arm: x_data_points = [origin[0], x_pendulum_arm[i]] y_data_points = [origin[1], y_pendulum_arm[i]] p.set_data(x_data_points, y_data_points) # Update plot with set_data # Update Pendulum Mass: m1 patch_center = x_pendulum_arm[i], y_pendulum_arm[i] c.center = patch_center # Same idea here, instead of drawing a new patch update its location # Update Mass Block Patch: r.set(xy=(x_vec[i] - width / 2, -height / 2)) # Update Drawing: fig.canvas.draw() # Save Frame: writerObj.grab_frame() from IPython.display import Video Video("/work/simulation.mp4", embed=True, width=640, height=480)