# import the libraries !apt update -y !apt install ffmpeg -y from mpl_toolkits.mplot3d import Axes3D import numpy as np import matplotlib.pyplot as plt import math import sympy from matplotlib.animation import FFMpegWriter from matplotlib.patches import Circle from matplotlib.patches import Rectangle from matplotlib import pyplot as plt import matplotlib.animation as manimation

# Question 1 #Derive the equations of motion for the following bead that is constrained to slide without friction #on the parabolic wire below, where m = 1kg and g = 9.8 m/s2. Use the Euler-Lagrange equation #with constraints. # Define the parameters for the simulation g=9.8 # gravity in [m/s^2] m=1 # Length [Kg] ## Part A # Write the EOM using symbolic too box ## Part B # Simulate this system using Euler’s method (time step = 0.001 seconds) in Python with # initial conditions, x = -1m, dx/dt = 0m/s, over the course of 10 seconds. # Define our time vector for simulation sim_t=10 #[sec] dt=0.001 #[sec] t_vec=np.arange(0,sim_t,dt) # our time vector in increments of dt # define vectors that will hold our state variables over time [x,xd;y,yd] and our constraints force x_vec=np.zeros(len(t_vec)) xd_vec=np.zeros(len(t_vec)) y_vec=np.zeros(len(t_vec)) yd_vec=np.zeros(len(t_vec)) Lamda_vec=np.zeros(len(t_vec)) Lam_vec=np.zeros(len(t_vec)) # Assign Intial Values: x_vec[0] = -1 xd_vec[0] = 0 y_vec[0] = 1 yd_vec[0] = 0 # Euler Simulation: Using Matrix Form (A * x = B) # Initialize A and B: A = np.array([[m, 0,0], [0, m ,-1],[0,-1,0]]) B = np.array([0,-m*g,0]) for i in range(1, len(t_vec)): # Only the off diagonal needs to be Updated: A[2, 0] = 2*x_vec[i-1] A[0, 2] = 2*x_vec[i-1] # b must be updated every iteration: B[2] =-2*np.power(xd_vec[i-1], 2) [xdd, ydd,Lamda_vec[i]] = np.linalg.solve(A, B) # Use ddx and ddtheta to solve: x_vec[i] = x_vec[i - 1] + xd_vec[i - 1] * dt+xdd*np.power(dt,2) y_vec[i] = np.power(x_vec[i-1],2) xd_vec[i] = xd_vec[i - 1] + xdd * dt yd_vec[i] = yd_vec[i - 1] + ydd * dt Lam_vec[i]=(-m*xdd)/2*x_vec[i-1] plt.plot(t_vec,x_vec) plt.xlabel('Time[sec]') plt.ylabel('Horizontal Position') plt.title('X(t)') plt.grid(True) plt.show() plt.plot(t_vec,y_vec) plt.xlabel('Time[sec]') plt.ylabel('Vertical Position') plt.title('Y(t) ') plt.grid(True) plt.show() plt.plot(t_vec,Lamda_vec) plt.xlabel('Time[sec]') plt.ylabel('Constrained Force') plt.title('Lamba') plt.grid(True) plt.show()

## Part C # Animate the motion from Problem 1b, showing a bead sliding on a wire (at 20 frames per # second). x = x_vec y = y_vec vec_size=len(t_vec) title = "simulation" # Setup Animation Writer: FPS = 20 sample_rate = int(1 / (dt * FPS)) # Real Time Playback dpi = 300 writerObj = FFMpegWriter(fps=FPS) # Initialize the movie fig = plt.figure() # plot the sine wave line sine_line, = plt.plot(x_vec, y_vec, 'b') red_circle, = plt.plot([], [], 'ro', markersize = 10) plt.xlabel('x') plt.ylabel('sin(x)') # Update the frames for the movie with writerObj.saving(fig, title + ".mp4", dpi): for i in range(0, vec_size, sample_rate): x0 = x_vec[i] y0 = y_vec[i] red_circle.set_data(x0, y0) # 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)

## Question 2 #In the following system, a frictionless pendulum of length L = 1m and bob mass m1 = 1kg. is #attached with a spring of linear stiffness, k = 10N/m, to a frictionless sliding mass, m2 = 2kg. The #pendulum hinge point is 1.5m above the ground, the rest length of the spring is 1m, and g = 9.8 #m/s2. ## Part A #Write the equations of motion which include d2x/dt2, d2y/dt2, and any constraint forces #Report: Attach your mathematical derivation using Lagrangian Mechanics. Derive with #original symbols (e.g. m1, not 1kg) (20 points) # Pendulum Parameters: m_1 = 1 # kg m_2 = 2 # kg L = 1 # m Ls=1 g = 9.8 # m/s^2 K=10 # Create Symbols for Time: t = sympy.Symbol('t') # Creates symbolic variable t # Create Generalized Coordinates as a function of time: q = [th, x] th = sympy.Function('th')(t) x = sympy.Function('x')(t) # Position Equation: r = [x, y] r1 = np.array([L * sympy.sin(th), -L * sympy.cos(th)+1.5]) # Position of first mass = pendulum r2=np.array([x,0]) # Position of block sliding print('This print the position of the pendulum block r1=',r1) print('This print the position the smaller mass r2=',r2) # Derive the change in the spring length Spring_x=(x-r1[0])**2 Spring_y=(r1[1])**2 Spring_total=Spring_x+Spring_y L_spring=Spring_total**0.5 D_spring=L_spring print('This print the change in length',D_spring) # Velocity Equation: d/dt(r) = [dx/dt, dy/dt] v1 = np.array([r1[0].diff(t), r1[1].diff(t)]) # Velocity of first pendulum v2 = np.array([r2[0].diff(t),0]) # Velocity of second pendulum print('This print the velocity of the pendulum block v1=',v1) print('This print the velocity the smaller mass v2=',v2) # Energy Equations: T = 1/2 * m_1 * np.dot(v1, v1) + 1/2 * m_2 * np.dot(v2, v2) # Kinetic Energy V_2G= m_1 * g * r1[1] # Potential Energy due to gravity V_2S= 0.5*K*(D_spring**2) # potential energy due to the spring V=V_2G+V_2S L = T - V # Lagrangian print('The lagrangian term L=',L) # Lagrange Terms: dL_dth = L.diff(th) dL_dth_dt = L.diff(th.diff(t)).diff(t) dL_dx = L.diff(x) dL_dx_dt = L.diff(x.diff(t)).diff(t) # Euler-Lagrange Equations: dL/dq - d/dt(dL/ddq) = 0 th_eqn = dL_dth - dL_dth_dt x_eqn = dL_dx - dL_dx_dt print('X equation of motion x_eqn',x_eqn) print('th equation of motion x_eqn',th_eqn) # Replace Time Derivatives and Functions with Symbolic Variables: replacements = [(th.diff(t).diff(t), sympy.Symbol('ddth')), (th.diff(t), sympy.Symbol('dth')), (th, sympy.Symbol('th')), (x.diff(t).diff(t), sympy.Symbol('ddx')), (x.diff(t), sympy.Symbol('dx')), (x, sympy.Symbol('x'))] th_eqn = th_eqn.subs(replacements) x_eqn = x_eqn.subs(replacements) r1 = r1[0].subs(replacements), r1[1].subs(replacements) r2 = r2[0].subs(replacements),r2[1] print('X equation of motion x_eqn',x_eqn) print('th equation of motion x_eqn',th_eqn) # Simplfiy then Force SymPy to Cancel factorization: [Sometimes needed to use .coeff()] th_eqn = sympy.simplify(th_eqn) x_eqn = sympy.simplify(x_eqn) th_eqn = th_eqn.cancel() x_eqn = x_eqn.cancel() print('X equation of motion x_eqn',x_eqn) print('th equation of motion x_eqn',th_eqn) # Solve for Coefficients for A * x = B where x = [ddth1 ddth2] A1 = th_eqn.coeff(sympy.Symbol('ddth')) A2 = th_eqn.coeff(sympy.Symbol('ddx')) A3 = x_eqn.coeff(sympy.Symbol('ddth')) A4 = x_eqn.coeff(sympy.Symbol('ddx')) print('This print from 96') print(A1) print(A2) print(A3) print(A4) # Multiply remaining terms by -1 to switch to other side of equation: A * x - B = 0 -> A * x = B remainder = [(sympy.Symbol('ddth'), 0), (sympy.Symbol('ddx'), 0)] B1 = -1 * th_eqn.subs(remainder) B2 = -1 * x_eqn.subs(remainder) print('This is my B matrix') print(B1) print(B2) ddth_fin=B1/A1 ddx_fin=B2/A4 print('This print from 122 my acceleration ') print('Angular acceleration=',ddth_fin) print('Linear acceleration',ddx_fin) # Generate Lambda Functions for A and B and Position Equations: replacements = (sympy.Symbol('th'), sympy.Symbol('dth'), sympy.Symbol('x'), sympy.Symbol('dx')) ddth_fin=sympy.utilities.lambdify(replacements, ddth_fin, "numpy") ddx_fin=sympy.utilities.lambdify(replacements, ddx_fin, "numpy") A1 = sympy.utilities.lambdify(replacements, A1, "numpy") A2 = sympy.utilities.lambdify(replacements, A2, "numpy") A3 = sympy.utilities.lambdify(replacements, A3, "numpy") A4 = sympy.utilities.lambdify(replacements, A4, "numpy") B1 = sympy.utilities.lambdify(replacements, B1, "numpy") B2 = sympy.utilities.lambdify(replacements, B2, "numpy") r1 = sympy.utilities.lambdify(replacements, r1, "numpy") r2 = sympy.utilities.lambdify(replacements, r2, "numpy") #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. Report: Plot of x(t), θ(t) for 10 seconds. Publish the code and plots together as a #PDF (20 points) # Simulate System: x0 = 1, 0, 0, 0 # th1, dth1, th2, dth2 dt = 0.001 sim_time = 10 Time = np.arange(0, sim_time, dt) sim_length = len(Time) # Initialize Arrays: th_vec = np.zeros(sim_length) dth_vec = np.zeros(sim_length) x_vec = np.zeros(sim_length) dx_vec = np.zeros(sim_length) # Evaluate Initial Conditions: th_vec[0] = x0[0] dth_vec[0] = x0[1] x_vec[0] = x0[2] dx_vec[0] = x0[3] # Initialize A and B: A = np.array([[0, 0], [0, 0]]) B = np.array([0, 0]) # Initialize Arrays: th1_vec=np.zeros(sim_length) dth1_vec=np.zeros(sim_length) x_vec=np.zeros(sim_length) dx_vec=np.zeros(sim_length) x1_vec = np.zeros(sim_length) y1_vec = np.zeros(sim_length) x2_vec = np.zeros(sim_length) y2_vec = np.zeros(sim_length) # Evaluate Initial Conditions: th1_vec[0]=x0[0] dth1_vec[0]=x0[1] x_vec[0]=x0[2] dx_vec[0]=x0[3] x1_vec[0], y1_vec[0] = r1(x0[0], x0[1], x0[2], x0[3]) x2_vec[0], y2_vec[0] = r2(x0[0], x0[1], x0[2], x0[3]) # Initialize A and B: A=np.array([[0,0],[0,0]]) B=np.array([0,0]) # Euler Integration: for i in range(1, sim_length): ddth1=ddth_fin(th1_vec[i-1], dth1_vec[i-1], x_vec[i-1], dx_vec[i-1]) ddx=ddx_fin(th1_vec[i-1], dth1_vec[i-1], x_vec[i-1], dx_vec[i-1]) # Euler Step Integration: th1_vec[i] = th1_vec[i-1] + dth1_vec[i-1] * dt dth1_vec[i] = dth1_vec[i-1] + ddth1 * dt x_vec[i] = x_vec[i-1] + dx_vec[i-1] * dt dx_vec[i] = dx_vec[i-1] + ddx * dt # Animation States: x1_vec[i], y1_vec[i] = r1(th1_vec[i-1], dth1_vec[i-1], x_vec[i-1], dx_vec[i-1]) x2_vec[i], y2_vec[i] = r2(th1_vec[i-1], dth1_vec[i-1], x_vec[i-1], dx_vec[i-1]) plt.xlabel('Time [s]') plt.ylabel('X(t) [m]') plt.title('X(t) [m] vs Time[s]') plt.grid() plt.plot(Time,x_vec) plt.show() plt.xlabel('Time [s]') plt.ylabel('Theta(t) [rad]') plt.title('Theta(t) vs Time') plt.grid() plt.plot(Time,th1_vec) plt.show()

# Animation # Set up our Figure for drawing our Pendulum fig, ax=plt.subplots() # Create a plot on these axes, which is currently empty p, = ax.plot([],[], color='green',linewidth=1) # initializes an empty plot ax.set_xlim([-2 ,2]) # X Lim ax.set_ylim([-2 ,2]) # Y Lim ax.set_xlabel('X') ax.set_ylabel('Y ') ax.set_title('Cart-Pendulum Simulation') video_title = "Cart-Pendulum Simulation" # We want to plot stuff on those axes c2= Circle((0,0), radius=0.1, color='cornflowerblue') ax.add_patch(c2) r=Rectangle((0,0),width=0.25,height=0.25,color='b') ax.add_patch(r) # Define information for our Animation FPS = 20 sample_rate= int(1/(FPS*dt)) dpi = 300 # Quality of the Video writerObj = FFMpegWriter(fps=FPS) sim_size= len(Time) # number of sim time points x_pend_arm=np.zeros(sim_size) y_pend_arm=np.zeros(sim_size) x_bl = np.zeros(sim_size) y_bl =np.zeros(sim_size) for i in range(0,sim_size): x_bl[i]=x2_vec[i] y_bl[i]= y2_vec[i] x_pend_arm[i]=x1_vec[i] y_pend_arm[i]= y1_vec[i] #Plot and Create Animation with writerObj.saving(fig,video_title+".mp4",dpi): # Loop through the frames we want to animate for i in range(0,sim_size,sample_rate): # Pendulum Arm xpt = [0,x_pend_arm[i],x_bl[i]] ypt = [1.5, y_pend_arm[i],y_bl[i]] # Update the data in plotting object p.set_data(xpt, ypt) #Update pendulum end circ_center=x_pend_arm[i],y_pend_arm[i] c2.set_center(circ_center) # Update our cart r_c = x_bl[i],y_bl[i] # ^ note: the commas without brackets create a LIST r.set_xy(r_c) #^ updates the cart #Update Drawing fig.canvas.draw() #Grab and save the frame writerObj.grab_frame() # Import the video you just made and run it on the notebook from IPython.display import Video Video('/work/simulation2.mp4', embed=True,width=640, height=480)

## Grad Students Only. A bead of zero radius and mass of 1kg slides without friction on the inside # surface of a paraboloid of equation, z = x2 + y2. # Part A # Write the equations of motion which include d2x/dt2, d2y/dt2, d2z/dt2 and any constraint forces (10 points) # Part B # Simulate this system using Euler’s method (time step = 0.001 seconds) in Python with # initial conditions, x = 1m, dx/dt = 0m/s, y = 0m, dy/dt = 1m/s over the course of 10 # seconds. Report: Plot of x(t), y(t) for 10 seconds. Publish the code and plots together as a # PDF (10 points) # Define the parameters for the simulation g=9.8 # gravity in [m/s^2] m=1 # Length [Kg] # Define our time vector for simulation sim_t=10 #[sec] dt=0.001 #[sec] t_vec=np.arange(0,sim_t,dt) # our time vector in increments of dt # define vectors that will hold our state variables over time [x,xd;y,yd] and our constraints force x_vec=np.zeros(len(t_vec)) dx_vec=np.zeros(len(t_vec)) y_vec=np.zeros(len(t_vec)) dy_vec=np.zeros(len(t_vec)) z_vec=np.zeros(len(t_vec)) dz_vec=np.zeros(len(t_vec)) lamb_vec=np.zeros(len(t_vec)) # Assign Intial Values: x_vec[0] = 1 dx_vec[0] = 0 y_vec[0] = 0 dy_vec[0] = 1 z_vec[0] = 1 dz_vec[0] = 0 for i in range(1, len(t_vec)): A = np.array([[m, 0, 0, -2*x_vec[i-1]], [0, m, 0, -2*y_vec[i-1]], [0, 0, m, 1], [-2*x_vec[i-1], -2*y_vec[i-1], 1, 0]]) B = np.array([0, 0, -m*g, 2*dx_vec[i-1]**2 + 2*dy_vec[i-1]**2]) [ddx, ddy, ddz, lamb] = np.linalg.solve(A, B) x_vec[i] = x_vec[i - 1] + dx_vec[i - 1] * dt y_vec[i] = y_vec[i - 1] + dy_vec[i - 1] * dt z_vec[i] = z_vec[i - 1] + dz_vec[i - 1] * dt dx_vec[i] = dx_vec[i - 1] + ddx * dt dy_vec[i] = dy_vec[i - 1] + ddy * dt dz_vec[i] = dz_vec[i - 1] + ddz * dt lamb_vec[i] = lamb # This plot the horizontal position of the beads over time plt.plot(t_vec,x_vec) plt.xlabel('Time[sec]') plt.ylabel('Horizontal Position[m]') plt.title('Horizontal Position vs Time') plt.grid(True) plt.show() # This plot the vertical position of the beads over time plt.plot(t_vec,y_vec) plt.xlabel('Time[sec]') plt.ylabel('Vertical Position[m]') plt.title('Vertical Position vs Time ') plt.grid(True) plt.show()

# Animation # Set up our Figure for drawing our Pendulum fig, ax=plt.subplots() fig = plt.figure() ax = Axes3D(fig) # Create a plot on these axes, which is currently empty p,=ax.plot3D([], [],[], lw=4, c='b', marker='o') p2,=ax.plot3D([], [],[], lw=3, c='r') ax.set_xlim([-2 ,2]) # X Lim ax.set_ylim([-2 ,2]) # Y Lim ax.set_zlim([-2 ,2]) # Z Lim #ax.axis('equal') ax.set_xlabel('X') ax.set_ylabel('Y ') ax.set_zlabel('Z ') ax.set_title('Bead Simulation') video_title = "simulation03" # Define information for our Animation FPS = 20 sample_rate= int(1/(FPS*dt)) dpi = 300 # Quality of the Video writerObj = FFMpegWriter(fps=FPS) # Now we're puttin bead and wire in right place # Initialize an array containing the position of the bead over time simulation_size= len(t_vec) # number of sim time points x_b = np.zeros(simulation_size) y_b =np.zeros(simulation_size) z_b=np.zeros(simulation_size) for i in range(0,simulation_size): x_b[i]=x_vec[i] y_b[i]= y_vec[i] z_b[i]=z_vec[i] # We've computed all the pendulum positions # Now we need to update the plot in a loop and store each of the frames in a video #Plot and Create Animation with writerObj.saving(fig,video_title+".mp4",dpi): # Loop through the frames we want to animate for i in range(0,simulation_size,sample_rate): #Bead path (wire) p2.set_data_3d(x_b, y_b,z_b) # Updated the data bead p.set_data_3d(x_b[i], y_b[i],z_b[i]) #Update Drawing fig.canvas.draw() #Grab and save the frame writerObj.grab_frame() # Import the video you just made and run it on the notebook from IPython.display import Video Video('/work/simulation03.mp4', embed=True,width=640, height=480)