Lab 5 - Curlss/Computational_physics GitHub Wiki
Our main objective in this lab was to observe and interpret experimental data related to a mass vibrating on a spring. To do this, we had to identify the analytic solution that described the motion, which we did by using the Euler and Euler-Cromer methods. We were able to model our oscillation motion using the Euler and Euler-Cromer numerical methods. By comparing the numerical solutions to the analytic solution, we were able to determine which of the two methods was better at modeling harmonic motion and how the Euler-Cromer method conserves energy.
This lab is divided into three sections. In part 1, we visualize data and figure out and depict the analytical solution for a basic harmonic oscillator. Our goal is to use the Euler technique to model our motion in part two and the Euler-Cromer approach to represent our motion in part three.
Part 1
To get an understanding of the facts and values related to basic harmonic motion, we began part 1 by importing actual experimental data for a mass-spring system. We were able to see the oscillating pattern by using the data in the position vs. time display.
We obtained values for the amplitude, angular frequency, and phase shift—the defining parameters of the SHM equation, which can be seen by x(t) = A cos(ωt - ϕ)—by examining the graph. ω is determined using the formula ω = 2π/T, ϕ is the oscillation's starting phase offset, and A is the peak displacement from equilibrium.
# write a function to calculate the position of the SHO
# input = A, omega, phi, t
# return = position
def pos_calculate(A, omega, phi, t):
'''
description: calculates the position of the SHO
parameters: A- amplitude, omega- angular frequency, phi- phase constant, t- time
return: position
'''
position = A*np.cos(omega*t - phi)
phi = np.pi/2
return position
In order to verify the correctness of the solution, we also plotted the analytical model of the oscillator that was created using these parameters against the experimental data.
analytic = pos_calculate(A,omega,-0.9,time)
plt.figure()
plt.plot(time,pos, 'bo', label= 'position')
plt.plot(time,analytic, 'r-', label= 'analytic')
plt.xlabel('time(s)')
plt.ylabel('position (m)')
plt.title('Position vs Time')
plt.legend()
plt.show()
We then calculated the initial position and velocity of the mass as well as the angular frequency and the spring constant divid by mass.
| Parameter | Value |
|---|---|
| $x_0$ | 0.29 |
| $v_0$ | -0.32 |
| $\omega$ | 8.38 |
| $k/m$ | 70.18 |
Part 2
In the second section, we modeled the mass on a spring using the Euler technique to see whether it could demonstrate energy conservation. Using arrays for each, we were able to estimate each subsequent interval by calculating acceleration, velocity, and position step-by-step over time using Newton's second law and Hooke's law (F = -kx) in order to compute the Euler method.
def calculate(time, position, velocity, k, m, dt):
"""
description: calculate the time, position, and velocity arrays using the Euler method
parameters: time, position, velocity, k, m, dt
return: time, position, velocity arrays
"""
# your code here
for i in range(1,len(position)):
a= acceleration(k,m,position[i-1])
position[i] = position[i-1] + velocity[i-1]*dt
velocity[i] = velocity[i-1] + a*dt
time[i] = time[i-1] + dt
return time, position, velocity
We set values and called our function to create a graph, to assess the accuracy of the Euler method, we used our analytical solution form. Our graph made it clear that the Euler method was inaccurate and not a suitable technique for simulating oscillatory motion.
define the values for the initial position and velocity
# use the values you determined from Part 1
x0 =0.029
vx0 =-0.32
# set the amplitude, mass and spring constant to the values from Part 1
A = 0.045
m = .15 #kg
k = 10.53
# Enter the period you determined from Part 1
T = 0.75
# 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
time, position, velocity = initialize(x0,vx0,tmax,dt)
# call calculate
time, position, velocity = calculate(time, position, velocity, k, m, dt)
# plot position versus time
# set the figure size to (8,6)
plt.figure(figsize=(8,6))
plt.plot(time,position, 'b.', label= 'Euler method')
plt.xlabel('time(s)')
plt.ylabel('position (m)')
plt.title('Position vs Time')
plt.legend()
plt.show()
# overplot your analytic solution from Part 1
plt.plot(time,pos_calculate(A,omega,-0.9,time), 'r.', label= 'analytic')
plt.plot(time,position, 'b.', label= 'Euler method')
plt.xlabel('time(s)')
plt.ylabel('position (m)')
plt.title('Position vs Time')
plt.legend()
plt.show()
In order to determine if the Euler approach conserves energy, we computed the system's mechanical energy (the sum of its kinetic and potential energy) over time. This experiment allowed us to see that the Euler technique did not conserve energy. The progressive rise or fall of energy using the Euler technique demonstrates that energy conservation does not exist. Use the Euler numerical approach.
E= 0.5*m*velocity**2 + 0.5*k*position**2
plt.figure()
plt.plot(time,E, 'b.', label= 'Euler method')
plt.xlabel('time(s)')
plt.ylabel('Energy (J)')
plt.title('Energy vs Time')
plt.legend()
plt.show()
Part 3
To determine if the Euler-Cromer approach would be more effective at simulating basic harmonic motion and save energy, we substituted it for the Euler technique in part 3. The Euler-Cromer technique updates velocity first, then updates position using the new velocity, in contrast to the Euler method, which updates position using the previous velocity. This modification significantly alters the situation, and we discovered that the Euler-Cromer approach is, in fact, dependable for simulating basic harmonic motion and energy conservation when we plotted it against the analytical answer. A constant energy number indicates a steady approach, which is what we saw with our Euler-Cromer technique. The energy vs. time plot monitors the system's total mechanical energy (kinetic + potential) during the course of the simulation.
def euler_cromer_calculate(time, position, velocity, k, m, dt):
"""
Description: calculate the time, position, and velocity arrays using the Euler-Cromer method
parameters: time, position, velocity, k, m, dt
return: time, position, velocity arrays
"""
for i in range(1,len(position)):
a= acceleration(k,m,position[i-1])
velocity[i] = velocity[i-1] + a*dt
time[i] = time[i-1] + dt
position[i] = position[i-1] + velocity[i]*dt
return time, position, velocity
We then plotted the velocity and position of the mass oscillating on a spring.
There was a final part in this lab that was counted as a bonus, and you were to include a damping force in the Euler-Cromer method from above. The code is similar to the above, with slightly differences in the equations.
def acceleration_d(k,m,b, position, velocity):
a = (-k * position - b * velocity)/m
return a
---------------------------------------------------------------------------------
def euler_cromer_damping_calculate(time, position, velocity, k, m, dt, b):
for i in range(1,len(position)):
a= acceleration_d(k,m,b,position[i-1], velocity[i-1])
velocity[i] = velocity[i-1] + a*dt
time[i] = time[i-1] + dt
position[i] = position[i-1] + velocity[i]*dt
return time, position, velocity
---------------------------------------------------------------------------------
def initialize(x0,vx0,tmax, dt):
nstep= int(tmax/dt) + 1
time = np.zeros(nstep)
position = np.zeros(nstep)
velocity = np.zeros(nstep)
position[0] = x0
velocity[0] = vx0
return time, position, velocity
---------------------------------------------------------------------------------
#define the values for the initial position and velocity
# use the values you determined from Part 1
x0 =0.029
vx0 =-0.32
# set the amplitude, mass and spring constant to the values from Part 1
A = 0.045
m = .15 #kg
k = 10.53
b = 5
# Enter the period you determined from Part 1
T = 0.75
# 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
time, position, velocity = initialize(x0,vx0,tmax,dt)
# call calculate
time, position, velocity = euler_cromer_damping_calculate(time, position, velocity, k, m, dt, b)
# plot position versus time
# set the figure size to (8,6)
plt.figure(figsize=(8,6))
plt.plot(time,position, 'b.', label= 'Euler method')
plt.xlabel('time(s)')
plt.ylabel('position (m)')
plt.title('Position vs Time')
plt.legend()
plt.show()
# overplot your analytic solution from Part 1
plt.plot(time,pos_calculate(A,omega,-0.9,time), 'r.', label= 'analytic')
plt.plot(time,position, 'b.', label= 'Euler method')
plt.xlabel('time(s)')
plt.ylabel('position (m)')
plt.title('Position vs Time')
plt.legend()
plt.show()