# MOD HW 1
# Question 1
# 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 matplotlib.pyplot as plt
from matplotlib.animation import FFMpegWriter
from matplotlib.patches import Circle
# Parameters
m = 1 # kg
g = 9.8 # m/s^2
# 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))
y_vec = np.zeros(len(t_vec))
lam_vec = np.zeros(len(t_vec))
# Set our initial condition
x_vec[0] = -1 # initial x position (m)
dx_vec[0] = 0 # initial x velocity (m/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] + (-2*x_vec[i-1]*(2*dx_vec[i-1]**2+g))/(4*x_vec[i-1]**2+1)*dt
lam_vec[i] = (2*m*dx_vec[i-1]**2+g*m)/(4*x_vec[i-1]**2+1)
# calculating y positions
for i in range(0, len(t_vec)):
y_vec[i] = x_vec[i]**2
# 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('y')
plt.title('y(t)')
plt.plot(t_vec,y_vec)
plt.show()
# Plotting conservative forces
plt.xlabel('t')
plt.ylabel('λ')
plt.title('Conservative Forces')
plt.plot(t_vec,lam_vec)
plt.show()
# Setup Figure: Initialize Figure / Axe Handles
fig, ax = plt.subplots()
p, = ax.plot([], [], color='cornflowerblue')
ax.axis('equal')
ax.set_xlim([-2.5, 2.5])
ax.set_ylim([-0.5, 1.5])
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_title('Question 1 Simulation:')
video_title = "Q1simulation"
# Now we want to plot stuff on those axes
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)
x_path = np.linspace(-1.5, 1.5, 1000)
y_path = np.zeros(len(x_path))
for i in range(0, len(x_path)):
y_path[i] = x_path[i]**2
fig1 = plt.plot(x_path,y_path)
plt.show()
# Plot and Create Animation:
with writerObj.saving(fig, video_title+".mp4", dpi):
for i in range(0, sim_length, sample_rate):
# Update Pendulum Patch:
patch_center = x_vec[i], y_vec[i]
c.center = patch_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)