Exercise 4: Numerical integration methods

PHY224H1F Exercise 4: Numerical integration methods
The spring and mass
You will get experimental data of a mass-spring system, and write a Python program to solve
the equation of motion. You will then use this simulation data to plot a visualization of position
vs. time and the phase plot. You will also calculate the energy and plot it. You will have to
discuss the output, and eventually improve the code.
This lab contains explicit questions for you to answer labelled with Q. However, you still need
to maintain a notebook for the lab, make plots and write Python code.
Background knowledge for exercise 4
Python lists, arrays, numpy, pyplot
Physics Harmonic oscillator, Hooke’s law. Review these from your first year text.
1 First Session
1.1 Introduction
Figure 1: A vertical spring-mass
system
Physics of springs is based on Hooke’s Law
Fspring = −ky , (1)
where y is the vertical displacement, and k is the spring constant measured in N/m. Hooke’s law can be used with Newton’s second law to
derive the equation of motion:
d
2y
dt
2
+
k
m
y = 0 . (2)
Q1 Show how equation (2) was derived. Point out the main approximations involved (if any).
The most common method used to find numerical solutions to equations of motion is by setting up a pair of coupled ordinary differential
equations.
Given a mass, m, coordinate ~q, momentum ~p, and force F~ ,
d~p
dt
= F~ (3a)
d~q
dt
=
~p
m
. (3b)
The force can be any function of ~p and ~q in general, but in this exercise, we will look at a
linear spring: F~ = F~
spring
1.2 Numerical Methods
The solution to equation (2) is well-known, but in general, differential equations can be difficult to solve. In
practice, we can calculate approximate solutions by approximating the derivatives as
dy
dt
≈
1
∆t

y(t + ∆t) − y(t)

. (4)
Page 1
PHY224H1F Exercise 4: Numerical integration methods
Therefore, equations (3) can be approximated with
~p(t + ∆t) = ~p(t) + ∆t F(~q(t)) (5a)
~q(t + ∆t) = ~q(t) + ∆t
m
~p(t). (5b)
Equations (5) are update formulae. Given initial position, momentum, and starting time t = t0, we can
determine a numerical approximation of position and momentum at times t0 + ∆t, 2∆t, . . .. The numerical
solution will approach the actual solution as ∆t → 0.
This time-stepping scheme is called the “Forward Euler” method, because it uses the derivative at the
start of an interval to extrapolate forward. We will see other time-stepping schemes in the next session.
The scheme can also be written in a recursive form,
~pi+1 = ~pi + ∆t F(~qi) (6a)
~qi+1 = ~qi +
∆t
m
~pi
, (6b)
where ~pi = ~p(ti), ~qi = ~q(ti), and ti = t0 + i∆t.
1.3 Numerical methods for the mass-spring system
The equation of motion for the vertical spring-mass system can be written in the coupled form as:
dv
dt
= −Ω
2
0
y (7a)
dy
dt
= v , (7b)
where Ω0 =
p
k/m is the angular frequency of oscillation.
The initial conditions will be y0 = A (the amplitude measured in experiment), and v0 = 0 m/s. Using (7),
the numerical approximation can be written as
yi+1 = yi + ∆t vi (8a)
vi+1 = vi − ∆t Ω
2
0
yi (8b)
for i = 0, 1, 2, . . . . This is the update scheme you will use in your first program.
Q2 Show the derivation of equations (8) using the equation of motion (2) and the definition of the Forward
Euler method (5).
1.4 The experiment
You should have a mass hanging from a spring. Underneath the mass is a motion sensor, which is plugged
into a data acquisition device (DAQ). The output of the DAQ is analyzed by a LabView application
(~\Desktop\2nd Year lab Files\MotionSensor.vi)
The goal of the lab exercise is to determing the period of oscillation of the bob, T (and therefore the
frequency Ω0). This will be used in our later analysis.
Do the following:
1. Measure the mass of the bob.
2. Open MotionSensor.vi.
3. Switch to the “Setup/Callibrate” panel, and change the samples/second to 100.
4. Switch to the “Collect” panel. Position the bob at about 20 cm above the motion sensor.
5. Start the bob oscillation by either raising or lowering it from equilibrium and letting go; a few
centimeters should be enough.
Page 2
PHY224H1F Exercise 4: Numerical integration methods
6. Click “start detector” the motion sensor should start buzzing now. Click “Collect Data” to begin the
data collection.
7. Let the bob oscillate for about 10 seconds, then stop the detector.
8. If the motion sensor is not aligned well, the data may be poor. As long as the graph shows a fairly
regular motion, it should be possible to measure the period.
9. Switch to the “Analyze” panel. Adjust the yellow cursors to the begining and end of 5 to 10 oscillations.
10. You may need to change the bounds of the plot to get a closer view. You do this by highlighting the
last value on either scale and replacing it by your value.
11. Read the time between cursors, and divide by the number of oscillations observed to get the period.
12. Record the period.
13. Open File/Position Data Set. Save the data file on your memory stick.
14. Plot Distance vs. Time with labelled axes. Ignore the ”Dist. Error” for now.. Calculate the period
from your plot.
Q3 What is the value of the spring constant, k?
1.5 Python programming
Scripts for simulations typically follow the same pattern
1. define constants/parameters
2. set the initial conditions
3. use the numerical approximation (8) to step forward in time
4. loop until the final time
5. plot the graph and interpret the result
Additional steps that may be useful include
• wrapping it in a function for repeated simulations,
• saving the data for later analysis.
The outline for our specific program is,
1. import required modules (eg. pylab, or numpy and matplotlib.pyplot)
2. define constants ∆t and Ω0.
3. calculate the time values, ti for the plot using ∆t = 0.01 s, t0 = 0 s, T = 10.0 s.
4. preallocate the array (with the zeros() function)
5. set initial conditions
6. write the time-stepping loop
7. plot the data (position vs. time) with labelled axes
Write and run the program. Save the program and its output on your memory stick.
1.6 Programming - Analysis
Q4 What do you see happening in the plot? Is this what you expect? Can you interpret the graph?
There are other graphs that are useful for analysis:
1. ‘velocity vs. time’ (time on horizontal axis, velocity on vertical axis)
2. ‘velocity vs. position’ (position on horizontal axis, velocity on vertical axis)
The latter is called a phase-plot.
Create these plots with labelled axes and save them.
Page 3
PHY224H1F Exercise 4: Numerical integration methods
1.6.1 Energy
Energy (and other conserved quantities) are very useful for analyzing dynamical systems. For the spring-mass
system, the mechanical energy is conserved,
Etot = K( ˙y) + U(y), (9)
where K is kinetic energy (K =
1
2mv2
) and U is potential energy (both elastic and gravitational). The
energy expression we will use is
Etot =
1
2
mv2 +
1
2
ky2
. (10)
The term for gavitational potential energy can be omitted by wisely choosing the reference for potential
energy. You should be able to prove this statement.
• Modify your program to calculate the energy at each step. Note that the energy is not zero at t = 0.
The mass (m) and spring constant (k) will need to be added to your code.
• Plot ‘Energy vs. time’
Q5 What does the energy plot suggest? Does this explain the strange appearance of the velocity vs. time
and velocity vs. position plots?
Q6 For a spring-mass system, the phase plot should be an ellipse. Explain why, using energy conservation.
Give an explanation for the appearance of your phase plot.
Q7 (bonus) Determine the leading error in our numerical method. Hint: use a Taylor expansion of y(t+∆t)
and find the terms we ignored in equation (5).
1.7 More on numerical integration schemes
If you properly implemented the numerical approximation in Equation (8), you should have found that the
simulation was bad: the amplitude increased, and energy was not conserved. The reason for the failure was
the integration method. We will discuss other numerical integrators, and implement a symplectic scheme,
which is suitable for conservative systems.
In the first program, we used the most primitive numerical integration method called the Forward Euler
or Explicit Euler method,
y[i + 1] = y[i] + ∆t v[i]
v[i + 1] = v[i] − ∆t Ω
2
0
y[i] .
Forward Euler extrapolates the position through an interval ∆t using only the (approximate) derivative
at the begining of the interval. It can be proved [1] that Forward Euler is unstable, which means that the
amplitude and total energy increases in time. In phase space, the increase in energy looks like a spiral
outward. The accuracy of the Forward Euler method can be improved by decreasing the time
step, but the energy will always increase.
A similar method to Forward Euler is Backward Euler or Implicit Euler. This method uses the derivative
at the end of the interval (making an implicit scheme) to calculate the next point,
y[i + 1] = y[i] + ∆t v[i + 1] (11a)
v[i + 1] = v[i] − ∆t Ω
2
y[i + 1] . (11b)
Backward Euler is absolutely stable, which means that longer timesteps can be used (especially for “stiff”
or rapidly changing ODEs). However, it artificially dissipates energy.
Page 4
PHY224H1F Exercise 4: Numerical integration methods
A simple remedy is to combine the Forward and Backward Euler methods to create the Euler-Cromer or
Symplectic Euler Method,
y[i + 1] = y[i] + ∆t v[i] (12a)
v[i + 1] = v[i] − ∆t
k
m
y[i + 1] . (12b)
The result is an explicit method, which conserves approximate energy. Do the following:
• Insert the new code lines for the symplectic method (12) into your previous program, and plot the
“Energy vs. time” and phase plots.
You may want to simulate using both methods, and plot both results on the same axes.
• Compare and discuss the plots from the Forward Euler and Symplectic Euler Methods. Use the same
time step for both methods in the comparison.
Q8 How do the energy and phase plots change using the symplectic method?
Keep in mind for the future that it is essential to use a method that conserves energy (e.g.
symplectic methods) when simulating conservative problems.
Your submission for this exercise should include
• All plots marked above (position, energy, and phase space plots for both Forward and Symplectic
Euler Methods)
• answers to questions 1-8,
• the final version of your program,
• values from your experiment (period, amplitude, frequency, and spring constant).
Page 5
PHY224H1F Exercise 4: Numerical integration methods
2 Second session
In the first part of this exercise, we analyzed the equation of motion of a spring-mass system and wrote
Python code to numerically integrate the equation. In this session, we will add a physical damping term,
and qualitatively compare the numerical simulation to the experiment.
2.1 Adding a damping term
In the real world, mechanical energy cannot be perfectly conserved. So far, the equations of motion you have
used model the system without damping (physical dissipation). In general, the damping force exerted on a
body depends on the velocity ~v of the body. Drag is one source of damping caused by velocity of a body
relative to the surrounding medium (air or water) commonly modelled as quadratic drag,
F~
d = −
1
2
CρA|~v|~v , (13)
where
C is the drag coefficient (dimensionless),
ρ is the density of the medium (kg/m3
),
A is the cross-sectional area perpendicular to the flow (m2
), and
~v is the velocity of the body relative to the medium (m/s).
The direction of the damping force is always opposite to the direction of velocity. The drag coefficient is not
a universal constant: it depends on viscosity of the medium, shape of the body, and roughness of its surface.
It also depends on the velocity through the Reynolds number.
The Reynolds number Re, is a very useful dimensionless quantity used to characterize flows in fluid
mechanics. In our case, it will characterize the dependence of the drag coefficient on velocity. The Reynolds
number is the ratio of the inertial force of the medium to the viscous force:
Re = ρvl
η
(14)
where
ρ is the density of the medium (kg/m3
),
v is the velocity of the body relative to the medium (m/s),
l is the characteristic length in the direction of flow (m), and
η is the dynamic viscosity (N·s/m2
).
If the flow is laminar, the Reynolds number takes small values (Re < 2300) and the drag coefficient is
inversely proportional to the velocity. This makes the drag force directly proportional to velocity:
F~
d = −γ~v , (15)
where γ is the damping coefficient.
When the flow is turbulent, the Reynolds number is large (Re > 4000) the drag coefficient is approximately
constant and the drag force is dependent on the square of the velocity.
Q9 What is the Reynolds number for the spring-mass system? Hint: this may depend on the amplitude of
the oscillation.
We expect that the motion and speed of a typical mass-spring system in air correspond to small Reynolds
numbers. Therefore, the equation of motion of our system would be written as:
d
2y
dt
2
+ γ
dy
dt
+ Ω2
0
y = 0 . (16)
Page 6
PHY224H1F Exercise 4: Numerical integration methods
2.2 Repeat the experiment
The effects of damping were likely not seen in the last session, because the damping for the bob is quite
small. To make damping more observable, we can attach a disk to the bob to increase drag.
• Attach the damping disk to the bob.Weigh the new system.
• Repeat the experiment as outlined in Section 1.4 for a longer period of time so that you can see
the decay – about 2 minutes with the fin attached. Note: because the experiment setup has
changed, you should measure the frequency and mass again this session.
• Using the cursors, take 5-6 readings of amplitude. Assuming an exponential decay of the oscillation
envelope, estimate the decay constant γ (γ/2 is the inverse of time for which amplitude falls to 1/e
of the initial value). You do not have to use Python to do this calculation. Please verify that the
exponent in position decay is γ/2 while the exponent in energy decay is γ.
• Save the position vs. time plot.
• Export the data in case you need it later.
2.3 The new Python program
The coupled equations we need to formulate for our new Python program are
dv
dt
= −Ω
2
0
y − γv (17a)
dy
dt
= v (17b)
Modify your program from part 1.3 to include drag with (17) and plot ‘Position vs. time’ and ‘energy vs.
time’. Use the decay coefficient determined experimentally.
Q10 Discuss the plot. What effect does damping has on the energy?
2.4 Compare with experimental data
The practice of Physics often includes these steps:
• come up with a (mathematical) model for a phenomenon,
• make predictions with that model,
• compare those predictions against experiment.
The comparisons can be qualitative (e.g. do both the predicted and measured amplitudes decay with time?),
or quantitative (e.g. what is the difference in period between the prediction and experiment?). For now,
we will just do a qualitative comparison. Fitting this complex function with curve_fit is very hard and
requires a carefully designed fitting routine to find the correct answer.
Q11 Compare the experimental data with the output of your program. What do you think makes them
different?
Your submission for this exercise should include:
• plots for position and energy of the damped system,
• a plot of the experimental data,
• answers to questions 9-11,
• the final version of your program,
• values from your experiment (period, frequency, decay rate, and spring constant).
Page 7
PHY224H1F Exercise 4: Numerical integration methods
References
[1] William H Press, Saul A Teukolsky, William T Vetterling, and Brian P Flannery. Numerical recipes in
C, volume 2. Cambridge university press Cambridge, 1996.
Prepared by John Ladan and Ruxandra Serbanescu, 2016. Edited by Christopher Lee, 2017.
Page 8