# Installing all necessary libraries for simulation and animation
!apt update -y
!apt install ffmpeg -y
# Problem One : Bead on a Wire
import numpy as np
import matplotlib.pyplot as plt
import sympy
import matplotlib.animation as animation
from mpl_toolkits import mplot3d
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FFMpegWriter
from matplotlib.patches import Circle
from matplotlib.patches import Rectangle
# Bead on a Wire: Simulation
# Problem Parameters
g = 9.8 # m/s^2
m = 1 # kg
# Creating the time vector and setting the time step
dt = 0.001
t_vec = np.arange(0,10,dt)
# Initialize a vector of zeros for x and y and any constraint forces (if any)
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))
lamda_vec = np.zeros(len(t_vec))
# Set up initial conditions
x_vec[0] = -1 # initial distance, m
dx_vec[0] = 0 # initial velocity, m/s
y_vec[0] = 1 # initial distance, m
dy_vec[0] = 0 # initial velocity, m/s
#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[0, 2] = 2 * x_vec[i-1]
A[2, 0] = 2 * x_vec[i-1]
# b must be updated every iteration:
B[2] = -2*np.power(dx_vec[i-1],2)
[ddx, ddy, lamda_vec[i]] = np.linalg.solve(A, B)
# Use ddx and ddy to solve:
x_vec[i] = x_vec[i - 1] + dx_vec[i - 1] * dt
y_vec[i] = np.power(x_vec[i-1],2)
dx_vec[i] = dx_vec[i - 1] + ddx * dt
dy_vec[i] = dy_vec[i - 1] + ddy * dt
plt.plot(t_vec,x_vec)
plt.xlabel('Time(sec)')
plt.ylabel('Horizontal Position (x)')
plt.title ('Horizontal Position over Time')
plt.show()
plt.plot(t_vec,y_vec)
plt.xlabel('Time(sec)')
plt.ylabel('Vertical Position (y)')
plt.title('Vertical Position over Time')
plt.show()
plt.plot(t_vec,lamda_vec)
plt.xlabel('Time(sec)')
plt.ylabel('Constraint Forces (λ)')
plt.title('Constraint Forces over Time')
plt.show()
# Bead on a Wire: Animation
# Setup Figure: Initialize Figure / Axe Handles
fig, ax = plt.subplots() # Initialize Figure and Axes
p, = ax.plot([], [], color='cornflowerblue') # Initialize Line Plot
ax.plot(x_vec, y_vec, color='cornflowerblue')
ax.axis('equal')
ax.set_xlim([-2, 2]) # X Lim
ax.set_ylim([-2, 2]) # Y Lim
ax.set_xlabel('X') # X Label
ax.set_ylabel('Y') # Y Label
ax.set_title('Bead on a Wire:') # Plot Title
video_title = "simulation" # name of the mp4
# Initialize Patch:
c = Circle((0, 0), radius=0.01, color='red') # Draws a circle of radius 0.1 at (X, Y) location (0, 0)
ax.add_patch(c) # Add the patch to the axes
# Setup Animation Writer:
FPS = 20 # Frames per second for the video
sample_rate = int(1 / (dt * FPS)) # Rate at which we sample the simulation to visualize
dpi = 300 # Quality of the Video
writerObj = FFMpegWriter(fps=FPS) # Video Object
# We need to calculate the cartesian location of the bead:
# Initialize Array:
simulation_size = len(t_vec) # We use this alot, so lets store it.
### Change what's necessary for bead and wire animation ###
x_parab = np.zeros(simulation_size)
y_parab = np.zeros(simulation_size)
origin = (-1,1) # where the bead will start
# Set up initial conditions
x_vec[0] = -1 # initial distance, m
dx_vec[0] = 0 # initial velocity, m/s
y_vec[0] = 1 # initial distance, m
dy_vec[0] = 0 # initial velocity, m/s
# Find the X and Y trajectories of the bead:
for i in range(0, simulation_size):
x_parab[i] = x_vec[i - 1] + dx_vec[i - 1] * dt
y_parab[i] = np.power(x_vec[i-1],2)
# Plot and Create Animation:
with writerObj.saving(fig, video_title+".mp4", dpi):
# We want to create video that represents the simulation
# So we need to sample only a few of the frames:
for i in range(0, simulation_size, sample_rate):
# Update Bead:
x_data_points = [x_parab[i]]
y_data_points = [y_parab[i]]
# We want to avoid creating new plots to make an animation (Very Slow)
# Instead lets take the plot we made earlier and just update it with new data.
p.set_data(x_data_points, y_data_points) # Update plot with set_data
# Update Bead Patch:
patch_center = (x_parab[i], y_parab[i])
c.center = patch_center # Same idea here, instead of drawing a new patch update its location
# Update Drawing:
fig.canvas.draw() # Update the figure with the new changes
# Grab and Save Frame:
writerObj.grab_frame()
# Saving the video
from IPython.display import Video
Video("/work/simulation.mp4", embed=True, width=640, height=480)
# You can find the video's path by right clicking on it in Notebooks & Files
# and selecting "Copy path to clipboard"
# Problem Two : Pendulum with Mass and Spring
# Pendulum with Mass and Spring: Simulation
import sympy
# Problem Parameters
g = 9.8 # m/s^2
L = 1 # m
m1 = 1 # pendulum, kg
k = 10 # linear stiffness, N/m
m2 = 2 # box, kg
hinge = 1.5 #m
springL_0 = 1 #m
# Create Symbols for Time:
t = sympy.Symbol('t') # Creates symbolic variable t
# Create Generalized Coordinates as a function of time: q = [theta, x]
th = sympy.Function('th')(t)
x = sympy.Function('x')(t)
# Position Equation
r1 = np.array([L * sympy.sin(th), (hinge - L*sympy.cos(th))]) # Position of spring pendulum
r2 = np.array([x]) # Position box
# Velocity Equation
v1 = np.array([r1[0].diff(t), r1[1].diff(t)]) # Velocity
v2 = np.array([r2[0].diff(t)]) # Velocity
# Finding Lf
Lf_init = ((x - L*sympy.sin(th))**2) + ((hinge - L*sympy.cos(th))**2)
Lf = (Lf_init**0.5)
# Energy Equations:
T = 1/2 * m1 * np.dot(v1, v1) + 1/2 * m2 * np.dot(v2,v2) # Kinetic Energy
V = m1 * g * r1[1] + 1/2* k*((Lf)**2) # Potential Energy
L = T - V # Lagrangian
# Lagrange Terms:
dL_dth = L.diff(th)
dL_dth_dt = L.diff(th.diff(t)).diff(t) # inside () is theta dot
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
# 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)
# 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()
# 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'))
# 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)
ddtheta_f=B1/A1
xdd=B2/A4
# Generate Lambda Functions for A and B and Position Equations:
replacements = (sympy.Symbol('th'), sympy.Symbol('dth'), sympy.Symbol('x'), sympy.Symbol('dx'))
ddtheta_f = sympy.utilities.lambdify(replacements, ddtheta_f, "numpy")
xdd = sympy.utilities.lambdify(replacements, xdd, "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 # th, dth, x, dx
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])
# Euler Integration:
for i in range(1, sim_length):
# Evaluate Dynamics:
ddth = ddtheta_f(th_vec[i-1], dth_vec[i-1], x_vec[i-1], dx_vec[i-1])
ddx = xdd(th_vec[i-1], dth_vec[i-1], x_vec[i-1], dx_vec[i-1])
# Euler Step Integration:
th_vec[i] = th_vec[i-1] + dth_vec[i-1] * dt
dth_vec[i] = dth_vec[i-1] + ddth * dt
x_vec[i] = x_vec[i-1] + dx_vec[i-1] * dt
dx_vec[i] = dx_vec[i-1] + ddx * dt
plt.plot(time,th_vec)
plt.xlabel('Time(sec)')
plt.ylabel('Angle (rad)')
plt.title ('Angle over Time')
plt.show()
plt.plot(time,x_vec)
plt.xlabel('Time(sec)')
plt.ylabel('Horizontal Position (x)')
plt.title ('Horizontal Position over Time')
plt.show()
# Pendulum with Mass and Spring: Animation
# Setup Figure: Initialize Figure / Axe Handles
fig, ax = plt.subplots()
p1, = ax.plot([], [], color='black', linewidth=2)
p2, = ax.plot([], [], color='black', linewidth=2)
lb, ub = -5, 5
ax.axis('equal')
ax.set_xlim([lb, ub])
ax.set_xlabel('X') # X Label
ax.set_ylabel('Y') # Y Label
ax.set_title('Pendulum with Mass and Spring Simulation')
video_title = "simulation"
y_offset = 0
# Draw the pendulum
pendulum = Circle((0, 3), radius=0.1, color='cornflowerblue', zorder=10)
ax.add_patch(pendulum)
# Initialize Patch: (Box)
width = 1 # Width of Box
height = width / 2 # Height of Box
xy_box= (x_vec[0] - width / 2, y_offset - height / 2) # Bottom Left Corner of box
r = Rectangle(xy_box, width, height, color='cornflowerblue') # Rectangle Patch
ax.add_patch(r) # Add Patch to Plot
# Draw the Ground:
ground = ax.hlines(-height / 2, lb, ub, colors='black')
height_hatch = 0.25
width_hatch = ub - lb
xy_hatch = (lb, y_offset - height / 2 - height_hatch)
ground_hatch = Rectangle(xy_hatch, width_hatch, height_hatch, facecolor='None', linestyle='None', hatch='/')
ax.add_patch(ground_hatch)
# Setup Animation Writer:
FPS = 20 # Frames per second for the video
sample_rate = int(1 / (dt * FPS)) # Rate at which we sample the simulation to visualize
dpi = 300 # Quality of the Video
writerObj = FFMpegWriter(fps=FPS) # Video Object
# Draw Static Objects:
pin_joint_pend = Circle((0, 3), radius=0.05, color='black', zorder=10)
ax.add_patch(pin_joint_pend)
# Animate:
with writerObj.saving(fig, video_title + ".mp4", dpi):
for i in range(0, len(t_vec), sample_rate):
# Update Pendulum Arm:
x_data_points = [x_vec[i], 0]
y_data_points = [th_vec[i], 0]
p.set_data(x_data_points, y_data_points)
# Update Pendulum Patch:
patch_center = x_vec[i], th_vec[i]
c.center = patch_center
# Update Box Patch:
r.set(xy=(x_vec[i] - width / 2, y_offset - height / 2))
# Update Drawing:
fig.canvas.draw()
# Save Frame:
writerObj.grab_frame()
# Saving the video
from IPython.display import Video
Video("/work/simulation.mp4", embed=True, width=640, height=480)
# You can find the video's path by right clicking on it in Notebooks & Files
# and selecting "Copy path to clipboard"
# Problem Three : Bead on a Paraboloid
# Bead on a Paraboloid: Simulation
# Problem Parameters
g = 9.8 # m/s^2
m = 1 # kg
# Creating the time vector and setting the time step
dt = 0.001
t_vec = np.arange(0,10,dt)
# Initialize a vector of zeros for x and y and z and any constraint forces (if any)
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))
lamda_vec = np.zeros(len(t_vec))
# Set up initial conditions
x_vec[0] = 1 # initial distance, m
dx_vec[0] = 0 # initial velocity, m/s
y_vec[0] = 0 # initial distance, m
dy_vec[0] = 1 # initial velocity, m/s
z_vec[0] = 1
dz_vec[0] = 0
#Euler Simulation: Using Matrix Form (A * x = B)
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, lamda] = np.linalg.solve(A, B)
# Use ddx and ddy to solve:
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
lamda_vec[i] = lamda
plt.plot(t_vec,x_vec)
plt.xlabel('Time(sec)')
plt.ylabel('Horizontal Position (x)')
plt.title ('Horizontal Position over Time')
plt.show()
plt.plot(t_vec,y_vec)
plt.xlabel('Time(sec)')
plt.ylabel('Vertical Position (y)')
plt.title('Vertical Position over Time')
plt.show()
fig = plt.figure()
# syntax for 3-D projection
ax = plt.axes(projection ='3d')
# defining all 3 axes
x = x_vec
y = y_vec
z = z_vec
# plotting
ax.plot3D(x, y, z, 'green')
plt.show()
# Create Animation:
# ANIMATION FUNCTION
def func(num, dataSet, line):
# NOTE: there is no .set_data() for 3 dim data...
line.set_data(dataSet[0:2, :num])
line.set_3d_properties(dataSet[2, :num])
return line
# THE DATA POINTS
x = x_vec
y = y_vec
z = z_vec
dataSet = np.array([x, y, z])
numDataPoints = len(z)
# GET SOME MATPLOTLIB OBJECTS
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
# NOTE: Can't pass empty arrays into 3d version of plot()
line = plt.plot(dataSet[0], dataSet[1], dataSet[2], lw=2, c='g')[0] # For line plot
# AXES PROPERTIES]
# ax.set_xlim3d([limit0, limit1])
ax.set_xlabel('X(t)')
ax.set_ylabel('Y(t)')
ax.set_zlabel('Z(t)')
ax.set_title('Bead on a Paraboloid')
# Creating the Animation object
line_ani = animation.FuncAnimation(fig, func, frames=numDataPoints*2, fargs=(dataSet,line), interval=50, blit=False)
anim = line_ani
line_ani.save('AnimationNew.mp4')
plt.show()