PHY224H1F Numerical integration methods
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit
Exercise 4: Numerical integration methods
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 ﬁrst year text.
1 First Session
Physics of springs is based on Hooke’s Law
Fspring = .ky ,
where y is the vertical displacement, and k is the spring constant mea- sured in N/m. Hooke’s law can be used with Newton’s second law to derive the equation of motion:
Q1 Show how equation (2) was derived. Point out the main approxima- tions involved (if any).
Given a mass, m, coordinate -q, momentum p-, and force ,
dt m .
Figure 1: A vertical spring-mass system
The force can be any function of p- and -q in general, but in this exercise, we will look at a linear spring: = spring
1.2 Numerical Methods
The solution to equation (2) is well-known, but in general, diﬀerential equations can be diﬃcult to solve. In practice, we can calculate approximate solutions by approximating the derivatives as
dt s ∆t ┌y(t + ∆t) . y(t)┐ .
Therefore, equations (3) can be approximated with
p-(t + ∆t) = p-(t) + ∆t F (-q(t))
-q(t + ∆t) = -q(t) + p t) .
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,
p-i+1 = p-i + ∆t F (q-i )
where p-i = p-(ti ), q-i = -q(ti ), and ti = t0 + i∆t.
The equation of motion for the vertical spring-mass system can be written in the coupled form as:
where Ω0 = ′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
vi+1 = vi . ∆t Ω0(2)yi
for i = 0, 1, 2, . . . . This is the update scheme you will use in your ﬁrst program.
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.
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 oscilla- tions.
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 ﬁle 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. deﬁne constants/parameters
2. set the initial conditions
3. use the numerical approximation (8) to step forward in time
4. loop until the ﬁnal 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 speciﬁc program is,
1. import required modules (eg. pylab, or numpy and matplotlib.pyplot)
2. deﬁne 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.
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 = mv2 ) and U is potential energy (both elastic and gravitational). The energy expression we will use is
Etot = mv2 + 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 ﬁnd 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 ﬁrst 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 Ω0(2)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  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]
v[i + 1] = v[i] . ∆t Ω2 y[i + 1] .
Backward Euler is absolutely stable, which means that longer timesteps can be used (especially for “stiﬀ” or rapidly changing ODEs). However, it artiﬁcially dissipates energy.
Symplectic Euler Method,
y[i + 1] = y[i] + ∆t v[i]
v[i + 1] = v[i] . ∆t y[i + 1] .
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 ﬁnal version of your program,
● values from your experiment (period, amplitude, frequency, and spring constant).
2 Second session
In the ﬁrst 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,
d = . CρAl-vl-v , (13)
C is the drag coeﬃcient (dimensionless),
ρ is the density of the medium (kg/m3 ),
A is the cross-sectional area perpendicular to the ﬂow (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 coeﬃcient 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 ﬂows in ﬂuid mechanics. In our case, it will characterize the dependence of the drag coeﬃcient on velocity. The Reynolds number is the ratio of the inertial force of the medium to the viscous force:
ρ 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 ﬂow (m), and
η is the dynamic viscosity (N ·s/m2 ).
If the ﬂow is laminar, the Reynolds number takes small values (Re < 2300) and the drag coeﬃcient is inversely proportional to the velocity. This makes the drag force directly proportional to velocity:
d = .γ-v , (15)
where γ is the damping coeﬃcient.
When the ﬂow is turbulent, the Reynolds number is large (Re > 4000) the drag coeﬃcient 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:
+ γ + Ω0(2)y = 0 .
2.2 Repeat the experiment
The eﬀects 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 ﬁn 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.
The coupled equations we need to formulate for our new Python program are
Q10 Discuss the plot. What eﬀect does damping has on the energy?