Goal

We want to visualize the motion of a spring-mass harmonic oscillator. To do so we import data of a mass-spring system and analyze it to determine necessary constants. Using those constant we apply the Euler Method to determine the position and velocity. This method gives us an incorrect approximation, so we use the Euler-Cramer method instead.

Overview

In order to visualize the motion of a mass-spring system, we attempt to utilize the Euler Method, then use the Euler-Cramer Method to obtain position and velocity. This has been done in the order indicated below:

• Visualize Data and Plot Solution for Simple Harmonic Oscillator
• Model Motion with Euler Method
• Model Motion with the Euler-Cromer Method
• Dampening of System

Visualize Data

We import data of a mass-spring system and use it to determine realistic constants we will use later. First we import the data:

!pip install wget
import wget

!wget https://raw.githubusercontent.com/rfinn/computational_physics/refs/heads/main/data/shm_mass_spring_good.csv

# use astropy to read in your data file
from astropy.io import ascii
data = ascii.read("shm_mass_spring_good.csv",delimiter=',')

# print data so you can see the column names
print(data)

# save the time and position data
t = data["Time(s)"]
x = data["Position(m)"]
v = data["Velocity(m/s)"]
# subtract the mean of the position from the position
x = x - np.mean(x)

Then we can plot this data:

# code to plot your position vs time data here
plt.plot(t,x,"bo")
plt.xlabel("t (s)")
plt.ylabel("x (m)")

Resulting in:

We can then determine the necessary constants we will use later: A, T, and omega. Amplitude (A) is the height from the initial position (0) to the peak. The period (t) is the length of time between the peaks. And omega is calculated using: $\omega = \frac{2\pi}{T}$. We estimated these values to be: 0.04, 0.8, and 7.85 respectively. We can then use these values to calculate the position analytically using the formula:

$$x = Acos(\omega t-\phi)$$

We make a function to calculate this, then plot the output value.

# write a function to calculate the position of the SHO
# input = A, omega, phi, t
# return = position
def pos(A,omega,phi,t):
  '''
  calculates position (x)
  '''
  x = A*np.cos(omega*t-phi)
  return x

# enter your code to plot your data and analytic solution
# be sure to use your function
phi = 0
pos = pos(A,omega,phi,t)
plt.figure()
plt.plot(t,pos,"bo",label="Analytic")
plt.plot(t,x,"ko",label="Data")
plt.xlabel("t (s)")
plt.ylabel("x (m)")
plt.legend()

Resulting plot:

As we can see, the analytic solution is more precise than the data.

Model Motion with Euler Method

Using the knowledge of the Euler Method we applied in previous labs, we have to create an initialize and a calculate function for this system, then we can plot the output values.

Initialize function:

# write your initialize function
# be sure to include your multi-line comment string

def initialize(x0,v0,tmax,dt):
  """
  description:
  sets empty arrays of (nsteps) indexes for t,x,v
  parameters:
  initial position (x0),initial velocity (v0),tmax,dt
  return:
  time (t), position (x), velocity (v) 

  """
  nsteps = int(tmax/dt)
  t = np.zeros(nsteps)
  x = np.zeros(nsteps)
  v = np.zeros(nsteps)
  x[0] = x0
  v[0] = v0
  return t,x,v

Calculate function:

# write your calculate function here
def calculate(t,x,v,k,m,dt):
  """
  description:
  fill empty arrays (t,x,v) with correct values 
  parameters:
  time (t), position (x), velocity (v), spring constant (k), mass (m), dt
  return:
  time (t), position (x), velocity (v)
  """

  # your code here
  for i in range(1,len(t)):
    t[i] = t[i-1]+dt
    v[i] = v[i-1]-acceleration(k,m,x[i-1])*dt
    x[i] = x[i-1]+v[i-1]*dt-0.5*acceleration(k,m,x[i-1])*dt**2
  return t,x,v

Setting initial values and plotting:

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

# set the amplitude, mass and spring constant to the values from Part 1
A = 0.04
m = 0.15
k = 9.25

# Enter the period you determined from Part 1
T = 0.8
# 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
t1,x1,v1 = initialize(x0,vx0,tmax,dt)


# call calculate
t1,x1,v1 = calculate(t1,x1,v1,k,m,dt)

# plot position versus time
# set the figure size to (8,6)
plt.figure(figsize=(8,6))
plt.plot(t1,x1,"ko",label="Euler")



# overplot your analytic solution from Part 1
plt.plot(t,x,"bo",label="Analytic")
plt.xlabel("time (s)")
plt.ylabel("position (m)")

plt.legend()

Resulting plot:

In this plot we see that the system is amplifying over time. This does not reflect reality and a more clearer example can be shown with total energy, which is given by:

$E = \frac{1}{2}mv^2+\frac{1}{2}kx^2$

We can calculate this value in our code and plot the total energy over time:

# Your plot of energy vs time here
E = 0.5*m*v1**2+0.5*k*x1**2
plt.plot(t1,E,"k.")
plt.xlabel("time (s)")
plt.ylabel("Energy (J)")

# Euler method is not a great approach. Total energy should not increase over time (not conserving).

We see that the total energy increases over time, which disobeys the law of conservation of energy.

Model Motion with the Euler-Cromer Method

As observed in the previous section, the Euler method failed to retrieve us accurate results for this system. To fix our results, we use the Euler-Cramer Method. This method involves a slight change to the Euler Method. In this case, an equation like this:

$x = x_0 + vt + \frac{1}{2}at^2$

becomes:

$x = x_0 +v_0t+\frac{1}{2}at^2$

We revise our calculate function with this new method:

# your euler_cromer_calculate

def euler_cromer_calculate(t,x,v,k,m,dt):
  """
  description:
  fill empty arrays (t,x,v) with revised values 
  parameters:
  time (t), position (x), velocity (v), spring constant (k), mass (m), dt
  return:
  time (t), position (x), velocity (v)
  """
  for i in range(1,len(t)):
    t[i] = t[i-1]+dt
    v[i] = v[i-1]-acceleration(k,m,x[i-1])*dt
    x[i] = x[i-1]+v[i]*dt-0.5*acceleration(k,m,x[i-1])*dt**2
  return t,x,v

Then we apply this function to the problem:

# 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 =0.04
vx0 = -0.32

# set the amplitude, mass and spring constant to the values from Part 1
A = 0.04
m = 0.15
k = 9.25

# Enter the period you determined from Part 1
T = 0.8
# 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
t2,x2,v2 = initialize(x0,vx0,tmax,dt)


# call calculate
t2,x2,v2 = euler_cromer_calculate(t2,x2,v2,k,m,dt)

# plot position versus time
# set the figure size to (8,6)
plt.figure(figsize=(8,6))
plt.plot(t2,x2,"ko")
plt.xlabel("time (s)")
plt.ylabel("position (m)")

Resulting plot:

And for total energy (replacing values from previous code):

We see in these plots now that the oscillation dampens over time and the total energy goes to zero. This reflects reality because we expect the mass-spring system to eventually return to rest.

Dampening the oscillator

We redefine the acceleration function to include a dampening force:

$F_{damp} = -bv$

def accel_damp(k,m,b,x,v):
  """
  description:
  calculates acceleration of the mass with dampening force
  parameters:
  spring constant (k), mass (m), position (x)
  return:
  acceleration (a)

  """
  # your code here
  a = k*x/m+b*v/m
  return a

Then we make a small edit to the calculate function to use the new acceleration function:

def euler_cromer_calculate_damp(t,x,v,k,m,dt):
  """
  description:
  fill empty arrays (t,x,v) with revised values for dampened
  parameters:
  time (t), position (x), velocity (v), spring constant (k), mass (m), dt
  return:
  time (t), position (x), velocity (v)
  """
  for i in range(1,len(t)):
    t[i] = t[i-1]+dt
    v[i] = v[i-1]-accel_damp(k,m,b,x[i-1],v[i-1])*dt
    x[i] = x[i-1]+v[i]*dt-0.5*accel_damp(k,m,b,x[i-1],v[i-1])*dt**2
  return t,x,v

And finally we set initial values, run the functions, and plot the results:

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

# set the amplitude, mass and spring constant to the values from Part 1
A = 0.04
m = 0.15
k = 9.25

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

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

# call initialize
t3,x3,v3 = initialize(x0,vx0,tmax,dt)


# call calculate
t3,x3,v3 = euler_cromer_calculate_damp(t3,x3,v3,k,m,dt)

# plot position versus time
# set the figure size to (8,6)
plt.figure(figsize=(8,6))
plt.plot(t3,x3,"ko",label="b=0.2")
plt.xlabel("time (s)")
plt.ylabel("position (m)")

b = 0.4

# call initialize
t3,x3,v3 = initialize(x0,vx0,tmax,dt)


# call calculate
t3,x3,v3 = euler_cromer_calculate_damp(t3,x3,v3,k,m,dt)

# plot position versus time
# set the figure size to (8,6)
plt.plot(t3,x3,"ro",label="b=0.4")
plt.xlabel("time (s)")
plt.ylabel("position (m)")

b = 1

# call initialize
t3,x3,v3 = initialize(x0,vx0,tmax,dt)


# call calculate
t3,x3,v3 = euler_cromer_calculate_damp(t3,x3,v3,k,m,dt)

# plot position versus time
# set the figure size to (8,6)
plt.plot(t3,x3,"bo",label="b=1")
plt.xlabel("time (s)")
plt.ylabel("position (m)")

plt.legend()

Resulting plot:

By varying the dampening constant, $b$, we observe that our revision is correct (higher constant equals more dampening of the system).

Chaotic Motion

Mass-Spring

We can add a dampening force and driving force to the system simultaneously to observe the resultant motion. Total acceleration for the mass-spring system thus becomes:

$$a_x=-\frac{k}{m}x-qv_x+F_Dsin(\Omega_Dt)$$

Then we edit our acceleration and calculate functions to reflect the new acceleration:

# acceleration function here
# write the multi-line comment at the beginning
def acceleration_fd(t,omega, x, q=0, vx=0, fd=0, omegad=0):
  """
  description:
  acceleration of the spring based on displacement with optional dampening, initial velocity, and driving force
  parameters:
  angular frequency (omega), displacement (x), dampening parameter (q), velocity (vx), driving force (fd), and driving angular frequency (omegad)
  return:
  acceleration (a)

  """
  # your code here
  a = -omega**2*x-q*vx+fd*np.sin(omegad*t)
  return a

def euler_cromer_calculate_fd(t,x,vx,dt,omega,q=0,fd=0,omegad=0):
  """
  description:
  fill empty arrays (t,x,v) with revised values
  parameters:
  time (t), position (x), velocity (v), dt, angular frequency (omega), dampening parameter (q), driving force (fd), and driving angular frequency (omegad)
  return:
  time (t), position (x), velocity (v)
  """
  for i in range(1,len(t)):
    t[i] = t[i-1]+dt
    v[i] = v[i-1]+acceleration_fd(t[i-1],omega,x[i-1],q=myq,vx=v[i-1],fd=FD,omegad=OmegaD)*dt
    x[i] = x[i-1]+v[i]*dt
  return t,x,v

Now we can use these functions and generate the position vs time graph and phase diagram:

# enter the parameters listed above
# set the amplitude, mass and spring constant
A = .05
m = 1
k = 1
# set the initial conditions to easy values
x0 = A
vx0 = 0
phi = 0

# calculate omega from k and m
omega = np.sqrt(k/m)
# calculate the period from omega
T = 2*np.pi/omega

# set the time step to a small fraction of the period
dt = T/100

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

# add a figure with a figsize = (10,4)
plt.figure(figsize=(10,4))
# adjust the spacing between plots
plt.subplots_adjust(wspace=.3)

# use q = 0.5 while we vary the values of FD
myq = 0.5

# set the frequency of the damping force, OmegaD, to 1/3 of the
# natural frequency omega
OmegaD = 1/3*omega

# create a list of values to loop over
FDlist = [0, 0.5, 1.2]
for FD in FDlist:
  # call initialize and save returned arrays as t, x, and v
  t,x,v = initialize(x0,vx0,tmax,dt)

  # call calculate
  t,x,v = euler_cromer_calculate_fd(t,x,v,dt,omega,q=myq,fd=FD,omegad=OmegaD)

  plt.subplot(1,2, 1)
  # plot position versus time, label each plot with value of q and FD
  plt.plot(t,x,label=f"q={myq}, fd={FD}")

  plt.subplot(1,2,2)
  # plot velocity vs position to create a phase-space diagram
  # label each plot with value of q and FD
  plt.plot(x,v,label=f"q={myq}, fd={FD}")


# overplot your analytic solution from Part 1
plt.subplot(1,2,1)
t_analytic = np.linspace(0,tmax,1000)
xeqn = A*np.cos(omega*t_analytic)
plt.plot(t_analytic,xeqn,ls='--',lw=1,label='Analytic Solution')

# add legend, labels and title for the left plot
plt.legend()
plt.xlabel("time (s)")
plt.ylabel("position (m)")
plt.title("dampening and driving force of mass-spring")

plt.subplot(1,2,2)
# add legend, labels and title for the right plot
plt.legend()
plt.xlabel("postion (m)")
plt.ylabel("velocity (m/s)")
plt.title("phase diagram")

Resultant plot:

The position vs time graph does not indicate chaotic motion, but it is shown in the phase diagram.

Pendulum

For a pendulum system we calculate a new acceleration:

$$a_(\theta)=-\omega^2sin(\theta)-qv_(\theta)+F_Dsin(\omega_Dt)$$

where $\omega^2=\frac{g}{l}$. Once again we change our acceleration and calculate functions:

# your accelerationP function here
def accelerationP(t,v, x, omega,q=0, fd=0, omegad=0):
  """
  description:
  acceleration of the pendulum based on displacement with optional dampening, initial velocity, and driving force
  parameters:
  angular frequency (omega), angle (theta), dampening parameter (q), driving force (fd), and driving angular frequency (omegad)
  return:
  acceleration (a)

  """
  # your code here
  a = -omega**2*np.sin(x)-q*v+fd*np.sin(omegad*t)
  return a

# write your calculate function here
def calculateP(t,x,v,omega, dt, q=0, fd=0, omegad=0):
  """
  description:

  parameters:

  return:
  """

  for i in range(1,len(t)):
    t[i] = t[i-1] + dt
    v[i] = v[i-1] + accelerationP(t[i-1],v[i-1], x[i-1],omega, q=myq, fd=FD, omegad=OmegaD) * dt
    x[i] = x[i-1] + v[i] * dt
    #if i == 1:
    #  print(i, x[i])
    if (x[i] > np.pi):
      x[i] = x[i] - 2*np.pi
    if (x[i] < -np.pi):
      x[i] = x[i] + 2*np.pi
    #if i == 1:
    #  print(i, x[i])
  return t, x, v

Then we create our plot:

# your main code here

# parameters from above
# initial conditions for the pendulum
theta0 = np.pi/6
omega0 = 0

# parameters of the pendulum
g = 9.8
l = 9.8
omega = np.sqrt(g/l)


# calculate period T from omega
T = 2*np.pi/omega

# set the time step to a small fraction of the period
dt = T/100

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

# we will use a damping constant of 0.5 and explore the impact of
# different driving forces
myq = 0.5

# we set the frequency of the driving force to 2/3 of the natural frequency
OmegaD = 2/3*omega

###   PLOTTING
# make a figure with figsize=(10,4)
plt.figure(figsize=(10,4))
# adjust spacing between plots
plt.subplots_adjust(wspace=.3)


# calculate the result for the following damping forces
FDlist = [0, 0.5, 1.2]

# loop over values of FD
for FD in FDlist:
  # call initialize
  t,x,v = initialize(theta0,omega0,tmax,dt)
  # call calculateP
  t,x,v = calculateP(t,x,v,omega,dt,q=myq,fd=FD,omegad=OmegaD)

  # plot position versus time in plt.subplot(1,2,1)
  plt.subplot(1,2,1)
  plt.plot(t,x,label=f"q = {myq} FD={FD}")
  # plot velocity vs time in plt.subplot(1,2,2)
  plt.subplot(1,2,2)
  plt.plot(x,v,label=f"q = {myq} FD={FD}")
# overplot your analytic solution
plt.subplot(1,2,1)
t_analytic = np.linspace(0,tmax,1000)
xeqn = A*np.cos(omega*t_analytic)

# label axes for subplot(1,2,1) using 16 point font
plt.xlabel("time (s)",fontsize=16)
plt.ylabel("position (m)",fontsize=16)
# add a title and legend for suplot(1,2,1)
plt.title("dampening and driving force of pendulum")
plt.legend()
plt.subplot(1,2,2)
# label axes for subplot(1,2,2) using 16 point font
plt.xlabel("position (m)",fontsize=16)
plt.ylabel("velocity (m/s)",fontsize=16)
# add a title (Phase Space Diagram) and legend for suplot(1,2,2)
plt.title("phase diagram")
plt.legend()

Resultant plot:

Chaotic motion is observed in both the position vs time plot and the phase diagram.