#Install ffmpeg (allows us to create videos) !apt update -y !apt install ffmpeg -y import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from matplotlib.animation import FFMpegWriter from matplotlib.patches import Circle from matplotlib.patches import Rectangle import pandas as pd import numpy as np import scipy as sp import sympy #from sympy import * # Problem 1 m=1 g=9.8 #m/s^2-1 dt=0.001 t=np.arange(0,10,dt) x_vec=np.zeros(len(t)) x_vec[0]=-1 x_dot=np.zeros(len(t)) x_dot[0]=0 y_vec=np.zeros(len(t)) y_vec[0]=1 y_dot=np.zeros(len(t)) y_dot[0]=0 A=np.array([[m,0,0],[0,m,-1],[0,-1,0]]) b=np.array([0,-m*g,0]) lamb=np.zeros(len(t)) for i in range(1,len(t)): A[0,2]=2*x_vec[i-1] A[2,0]=2*x_vec[i-1] b[2]=-2*np.power(x_dot[i-1],2) [x_ddot,y_ddot,lamb[i]]=np.linalg.solve(A,b) x_dot[i]=x_dot[i-1]+x_ddot*dt x_vec[i]=x_vec[i-1]+x_dot[i-1]*dt+x_ddot*np.power(dt,2) y_dot[i]=y_dot[i-1]+y_ddot*dt y_vec[i]=np.power(x_vec[i],2) plt.xlabel('time [s]') plt.ylabel('x(t) [m]') plt.title('Plot of x(t)') plt.plot(t,x_vec) plt.show() plt.xlabel('time [s]') plt.ylabel('y(t) [m]') plt.title('Plot of y(t)') plt.plot(t,y_vec) plt.show() plt.xlabel('time [s]') plt.ylabel('lambda(t) [N]') plt.title('Plot of lambda(t)') plt.plot(t,lamb) 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='cornflowerblue') # initializes an empty plot ax.set_xlim([-3 ,3]) # X Lim ax.set_ylim([-3 ,3]) # Y Lim ax.axis('equal') ax.set_xlabel('X') ax.set_ylabel('Y ') ax.set_title('Bead Simulation') video_title = "simulation" # We want to plot stuff on those axes c= Circle((0,0), radius=0.1, color='cornflowerblue') ax.add_patch(c) # 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) # number of sim time points x_bead = np.zeros(simulation_size) y_bead =np.zeros(simulation_size) for i in range(0,simulation_size): x_bead[i]=x_vec[i] y_bead[i]= y_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 wire p.set_data(x_bead, y_bead) # Update our bob circle (path) patch_center = x_bead[i],y_bead[i] # ^ note: the commas without brackets create a LIST c.center=patch_center #^ updates the circle #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/simulation.mp4', embed=True,width=640, height=480)

# Problem 2 # System Parameters L=1 # meters m1= 1 # kg m2 = 2 #kg k = 10 # N/m g=9.8 #m/s^2 # Create Symbols for Time: t=sympy.Symbol('t') # Create Generalized Coordinates as a function of time: q=[x,theta] th1=sympy.Function('th1')(t) x=sympy.Function('x')(t) #Position Equation: r=[x,y] r1=np.array([L*sympy.sin(th1),1.5-L*sympy.cos(th1)]) #Position of pendulum r2=np.array([x,0]) spring_dist=np.power(x-L*sympy.sin(th1),2)+np.power(1.5-L*sympy.cos(th1),2) bloop=0 total_spring=np.power(np.power(spring_dist,0.5)-bloop,2) v_grav=m*g*r1[1] #Velocity Equation: d/dt(r) = [dx/dt,dy/dt] v1=np.array([r1[0].diff(t),r1[1].diff(t)]) # velocity of pendulum v2=np.array([r2[0].diff(t),0]) # velocity of block # Energy Equation: T= 1/2*m2*np.dot(v2,v2)+1/2*m2*L*np.dot(v1,v1) V=v_grav+1/2*k*total_spring L=T-V # Lagrangian # Lagrange Terms: dL_dth1=L.diff(th1) dL_dth1_dt=L.diff(th1.diff(t)).diff(t) dL_dx=L.diff(x) dl_dx_dt=L.diff(x.diff(t)).diff(t) #Euler-Lagrange Equations: dL/d1 -d/dt(dL/ddq) = 0 th1_eqn= dL_dth1-dL_dth1_dt x_eqn=dL_dx -dl_dx_dt # Replace Time Derivatives and Functions with Symbolic Variables: replacements=[(th1.diff(t).diff(t),sympy.Symbol('ddth1')),(th1.diff(t),sympy.Symbol('dth1')),(th1,sympy.Symbol('th1')), (x.diff(t).diff(t),sympy.Symbol('ddx')),(x.diff(t),sympy.Symbol('dx')),(x,sympy.Symbol('x'))] th1_eqn=th1_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] # Simplify then Force SymPy to Cancel factorization: [Sometimes needed to use .coeff()] th1_eqn=sympy.simplify(th1_eqn) x_eqn=sympy.simplify(x_eqn) th1_eqn=th1_eqn.cancel() x_eqn=x_eqn.cancel() # Solve for Coefficients for A*x=B where x=[ddth1 ddx] A1=th1_eqn.coeff(sympy.Symbol("ddth1")) A2=th1_eqn.coeff(sympy.Symbol("ddx")) A3=x_eqn.coeff(sympy.Symbol("ddth1")) A4=x_eqn.coeff(sympy.Symbol("ddx")) # Multiply remaining terms by -1 to switch to other side of equation: A * x -B = 0 -> A * x = B remainder = [(sympy.Symbol('ddth1'),0),(sympy.Symbol('ddx'),0)] B1=-1*th1_eqn.subs(remainder) B2=-1*x_eqn.subs(remainder) ddthfinal=B1/A1 ddxfinal=B2/A4 # Generate Lambda FUnctions for A and B and Position Equations: replacements = (sympy.Symbol('th1'), sympy.Symbol('dth1'), sympy.Symbol('x'), sympy.Symbol('dx')) ddthfinal=sympy.utilities.lambdify(replacements, ddthfinal, "numpy") ddxfinal=sympy.utilities.lambdify(replacements, ddxfinal, "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")

# Simulate System: x0= 1,0,0,0 dt=0.001 sim_time=10 time=np.arange(0,sim_time,dt) sim_length=len(time) # 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): # Evaluate Dynamics: #A[0, 0] = A1(th1_vec[i-1], dth1_vec[i-1], x_vec[i-1], dx_vec[i-1]) #A[0, 1] = A2(th1_vec[i-1], dth1_vec[i-1], x_vec[i-1], dx_vec[i-1]) #A[1, 0] = A3(th1_vec[i-1], dth1_vec[i-1], x_vec[i-1], dx_vec[i-1]) #A[1, 1] = A4(th1_vec[i-1], dth1_vec[i-1], x_vec[i-1], dx_vec[i-1]) #B[0] = B1(th1_vec[i-1], dth1_vec[i-1], x_vec[i-1], dx_vec[i-1]).astype('int64') #B[1] = B2(th1_vec[i-1], dth1_vec[i-1], x_vec[i-1], dx_vec[i-1]).astype('int64') ddth1=ddthfinal(th1_vec[i-1], dth1_vec[i-1], x_vec[i-1], dx_vec[i-1]) ddx=ddxfinal(th1_vec[i-1], dth1_vec[i-1], x_vec[i-1], dx_vec[i-1]) #[ddth1, ddx] = np.linalg.solve(A, B) # 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('Plot of x(t)') plt.plot(time,x_vec) plt.show() plt.xlabel('time [s]') plt.ylabel('theta(t) [m]') plt.title('Plot of theta(t)') 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='black') # initializes an empty plot ax.set_xlim([-3 ,3]) # X Lim ax.set_ylim([-3 ,3]) # Y Lim ax.axis('equal') ax.set_xlabel('X') ax.set_ylabel('Y ') ax.set_title('Cart-Pendulum Simulation') video_title = "simulation2" # 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='red') 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) # Now we're puttin bead and wire in right place # Initialize an array containing the position of the bead over time simulation_size= len(time) # number of sim time points x_pendulum_arm=np.zeros(simulation_size) y_pendulum_arm=np.zeros(simulation_size) x_block = np.zeros(simulation_size) y_block =np.zeros(simulation_size) for i in range(0,simulation_size): x_block[i]=x2_vec[i] y_block[i]= y2_vec[i] x_pendulum_arm[i]=x1_vec[i] y_pendulum_arm[i]= y1_vec[i] print("Test",x_pendulum_arm[10],y_pendulum_arm[10]) # 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): # Pendulum Arm x_data_points = [0,x_pendulum_arm[i],x_block[i]] y_data_points = [1.5, y_pendulum_arm[i],y_block[i]] # Update the data in plotting object p.set_data(x_data_points, y_data_points) #Upd