!apt update -y
!apt install ffmpeg -y
import numpy as np
import sympy
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from matplotlib.animation import FFMpegWriter
# Parameters
m = 1 # kg
g = 9.8 # m/s^2
# time vector
dt = 0.001 # sec
t_vec = np.arange(0,10,dt)
# State variable vectors initialized to zero
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))
lam = np.zeros(len(t_vec))
# Initial Conditions
x_vec[0] = -1 # initial x-position (m)
dx_vec[0] = 0 # x-velocity (m/s)
y_vec[0] = (x_vec[0])**2 # initial y-position (m)
dy_vec[0] = 0 # intial y-velocity (m/s)
# Initialize A and B:
A = np.array([[m,0,-2*x_vec[0]],[0,m,-1],[2*x_vec[0],-1,0]])
B = np.array([[0], [-m*g], [-2*(dx_vec[0])**2]])
# Euler Simulation with A*x=b:
for i in range(1,len(t_vec)):
x_vec[i] = x_vec[i-1]+dx_vec[i-1]*dt
y_vec[i] = y_vec[i-1]+dy_vec[i-1]*dt
# update A and B where necessary
A[0,2] = -2*x_vec[i-1]
A[2,0] = -2*x_vec[i-1]
B[2] = -2*(dx_vec[i-1])**2
# solve for ddx, ddy, and lam
[ddx,ddy,lam[i-1]] = np.linalg.solve(A,B)
dx_vec[i] = dx_vec[i-1] + (-2*lam[i-1]*x_vec[i-1])*dt
dy_vec[i] = dy_vec[i-1] + ((-m*g+lam[i-1])/m)*dt
plt.plot(t_vec,x_vec)
plt.show()
plt.plot(t_vec,y_vec)
plt.show()
plt.plot(t_vec,lam)
plt.show()
# Animation
# Initialize Figure / Axe Hanndles
fig, ax = plt.subplots() # Initialize Figure and Axes
p, = ax.plot([],[],color='cornflowerblue') # Initialize Empy Plot
ax.axis('equal')
ax.set_xlim([-2,2])
ax.set_ylim([0,2])
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_title('Question 1 Simulation:')
video_title = "MOD-HW1-Q1-Simulation"
# Initialize Patch:
c = Circle((-1,1), 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)
simulation_size = len(t_vec)
#plotting parabola line
x = np.linspace(-10, 10, 500)
y = x**2
parabola_line, = plt.plot(x,y)
with writerObj.saving(fig,video_title+".mp4",dpi):
for i in range(0,simulation_size,sample_rate):
patch_center = x_vec[i], y_vec[i]
c.center = patch_center
fig.canvas.draw()
writerObj.grab_frame()
from IPython.display import Video
Video("/work/MOD-HW1-Q1-Simulation.mp4",embed=True,width=640,height=480)