# Simulation and animation of a bead on a wire constrained to y=x^2 # Model based robotics HW1 Problem 1 # Spring 2022 # Code adapted from class pendulum example # 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 g = 9.8 m = 1 # We need an array of time points from 0 to 10 in increments of 0.01 seconds dt = 0.001 t_vec = np.arange(0,10,dt) # 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)) lambda_vec = np.zeros(len(t_vec)) # Set our initial condition x_vec[0] = -1 # initial position dx_vec[0] = 0 # initial velocity y_vec[0] = x_vec[0]**2 ddx = -(m*g - 2*m*dx_vec[0]**2)/(2*m*x_vec[0] + 0.5*m/x_vec[0]) lambda_vec[0] = m*ddx/(2*x_vec[0]) # Loop through time # Euler's Method (approximately integrates the differential equation) for i in range(1, len(t_vec)): ddx = -(m*g - 2*m*dx_vec[i-1]**2)/(2*m*x_vec[i-1] + 0.5*m/x_vec[i-1]) x_vec[i] = x_vec[i-1] + dx_vec[i-1]*dt dx_vec[i] = dx_vec[i-1] + ddx*dt y_vec[i] = x_vec[i]**2 ddx = -(m*g - 2*m*dx_vec[i]**2)/(2*m*x_vec[i] + 0.5*m/x_vec[i]) lambda_vec[i] = m*ddx/(2*x_vec[i]) plt.plot(t_vec,x_vec) plt.title("X over time") plt.show() plt.plot(t_vec,y_vec) plt.title("Y over time") plt.show() plt.plot(t_vec,lambda_vec) plt.title("Constraint forces over time") plt.show()

# Set up our figure for drawing our pendulum fig, ax = plt.subplots() # create a plot on those 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') # x label ax.set_ylabel('Y') # y label ax.set_title('Bead on a Wire Simulation') video_title = "Simulation" # Now we want to plot stuff 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 putting the rod and bob in the right place # initialize an array containing the positions of the pendulum simulation_size = len(t_vec) # 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] = x_bead[i]**2 # 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): # update the pendulum arm x_data_points = [x_bead[i]] y_data_points = [y_bead[i]] # p.set_data(x_data_points,y_data_points) # update our bob circle patch_center = x_bead[i], y_bead[i] # ^ note: the commas without brackets create a LIST c.center = patch_center # update drawing fig.canvas.draw() # grab and save 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)