# 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()