Goal

Effectively model simple harmonic motion using the Euler method compared to actual data.

Overview

In this lab we collected data from a spring oscillating when a mass is applied to it, and used the data to model it using plots. We also used the Euler method to model this as well to see its effectiveness.

Part 1: Evaluating using data

Data Collection and uploading

To collect our data we used a ring stand with a spring attached to it, we then added a 100g weight to the spring and oscillated it. We collected the data using a censor and it was recorded on our computer.

After collecting the data we then ran it in python using the following code.

from google.colab import files

uploaded = files.upload()

for fn in uploaded.keys():
  print('User uploaded file "{name}" with length {length} bytes'.format(
      name=fn, length=len(uploaded[fn])))

Then to read in our data we used this code.

from astropy.table import Table
data = Table.read("csv-export.csv")
data.colnames

And assigned the time and position arrays using our data.

time = data['Data Set 1:Time(s)']
position = data['Data Set 1:Position(m)']

# subtract the mean of the position from the position
position_sub=position-np.mean(position)

Visualizing and Analyzing the data

To start we made a plot of position vs time which as shown below. We estimated our amplitude and period to be .82m and 1.2s, respectively, by looking at our plot and roughly estimating where they appeared to be.

We had to calculate $\omega$ using the equation $\omega=\frac{2\pi}{T}$. We used this to get our analytic solution.

To determine our analytical solution we created a function to calculate it.

def calc_sho(A,omega,phi,t):
  """This function takes in variables A, omega, phi, and t and returns position for the spring"""
  position= A*np.cos(omega*t-phi)

  return position

We used trial and error to find the right $\phi$ value, which was found to be $\pi$. and plotted the analytic solution with our data.

We found our k to be 2.4, our $x_0$ to be -.1, and our $v_i$ to be -.07

Part 2: Using the Euler method

Functions

First we wrote a function to calculate acceleration.

def acceleration(k,m,dx):
  """
  description:This equation calculates acceleration using Hooke's law

  parameters:Spring constant (k), mass (m), displacement from equilibrium (dx)

  return: acceleration


  """
  # your code here
  a=(-k/m)*dx
  return a

We then wrote a function to initialize time, position and velocity arrays.

def initialize(x0,v0,tmax,dt):
  """
  description: This function calculates the number of the time steps, and creates arrays for time, position, and velocity

  parameters: Initial x position (x0), intital velocity (v0), total time(tmax), and time steps (dt)

  return: Time array, position array, velocity array


  """
  nsteps=int(tmax/dt)
  pos=np.zeros(nsteps)
  vel=np.zeros(nsteps)
  time=np.zeros(nsteps)
  pos[0]=x0
  vel[0]=v0

  return time, pos, vel

We then wrote a function to calculate our Euler solution.

def calculate(time,pos,vel,k,m,dt):
  """
  description: Calculates velocity time and position arrays

  parameters: Time (s), position (m), velocity (m/s), spring constant (k), mass (kg), time steps (dt)

  return:Time, position, velocity arrays
  """

  # your code here
  for i in range(1,len(pos)):
    pos[i]=pos[i-1]+vel[i-1]*dt
    vel[i]=vel[i-1]+acceleration(k,m,pos[i-1])*dt
    time[i]=time[i-1]+dt

  return time, pos, vel

Main Function

In the main function we assign variables with their respective values, and call our initialize and calculate functions in order to make our plot.

# define the values for the initial position and velocity
# use the values you determined from Part 1
x0 =-.1
vx0 =-.07

# set the amplitude, mass and spring constant to the values from Part 1
A =.1
m =.1
k =2.4

# Enter the period you determined from Part 1
T =1.28
# set the time step to a small fraction of the period
dt = T/100

# set tmax to 5x the period
tmax =5*T


# call initialize

t,pos,vel=initialize(x0,vx0,tmax,dt)

# call calculate
t_calc,pos_calc,vel_calc=calculate(t,pos,vel,k,m,dt)

# plot position versus time
# set the figure size to (8,6)

plt.figure()
plt.plot(t_calc,pos_calc,label='Motion of mass on a spring')
plt.title('Position vs Time')
plt.xlabel('Time (s)')
plt.ylabel('Position (m)')

We then overplot this with our analytical solution function.

def calc_sho(A,omega,phi,t):
  """This function takes in variables A, omega, phi, and t and returns position for the spring"""
  position= A*np.cos(omega*t-phi)

  return position
T = 1.28
omega = np.pi*2/T
phi = np.pi
A = .1
x_predict=calc_sho(A,omega,phi,t_calc)
plt.plot(t_calc,x_predict,label='Analytic Solution')
plt.legend()

The final plot looks as follows.

Conservation of energy

For a one-dimensional oscillator that is moving only in the x direction, the total energy is: $E=\frac{1}{2}mv_x^2 +\frac{1}{2}kx^2$ We then used the equation to plot energy vs time which gave us the following. This shows that the Euler solution does not conserve energy, and that the Euler method is not a good approach for modeling oscillatory motion.

Euler-Cromer solution

For this we made a new function for the Euler-Cromer calculation.

def euler_cromer_calculate(time,pos,vel,k,m,dt):
  """
  description: Calculates velocity time and position arrays

  parameters: Time (s), position (m), velocity (m/s), spring constant (k), mass (kg), time steps (dt)

  return:Time, position, velocity arrays
  """
  for i in range(1,len(pos)):
    pos[i]=pos[i-1]+vel[i]*dt
    vel[i]=vel[i-1]+acceleration(k,m,pos[i-1])*dt
    time[i]=time[i-1]+dt
  return time, pos, vel

The difference is in the position array, which loops over the same velocity as position instead of 1 before it.

We then called our Eulor-Cromer function.

# copy your main function cell from above
# call euler_cromer_calculate instead of calculate

# define the values for the initial position and velocity
# use the values you determined from Part 1
x0 =-.1
vx0 =-.07

# set the amplitude, mass and spring constant to the values from Part 1
A =.1
m =.1
k =2.4

# Enter the period you determined from Part 1
T =1.28
# set the time step to a small fraction of the period
dt = T/100

# set tmax to 5x the period
tmax =5*T


# call initialize

t,pos,vel=initialize(x0,vx0,tmax,dt)

# call calculate
t_calc1,pos_calc1,vel_calc1=euler_cromer_calculate(t,pos,vel,k,m,dt)

We then plotted our Euler-Cromer and our Euler plots side by side in order to compare them.

By observing these we can clearly see that the Euler-Cromer conserves energy while the Euler does not.