In this lecture, we will explore applying the techniques for solving initial value problems that we learned in lecture 11. We will also learn how to use numerical integration of ODEs to solve boundary value problems.
Learning objectives: After completing this lecture, you should be able to:
# Notebook code
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from scipy.signal import find_peaks
import sys
sys.path.append('resource/asnlib/public/')
sys.path.append('../answer_checking/')
from validate_answers import *
# Notebook code
plt.rcParams['figure.dpi'] = 100
The first physical problem we will be interested in simulating is the driven, damped Harmonic oscillator.
The physics of a harmonic oscillator is described by the motion of a mass free to move in one dimension that is attached to a rigid wall by a spring exerts a restoring force:
$$ F_r = -kx $$In general, a mass will also experience a friction drag force proportional to and opposing its velocity:
$$ F_f = - c \frac{dx}{dt} $$where $c$ is the friction coefficient. Including including a time dependent external force $F_{e}(t)$ applied to the mass, the total force becomes:
$$ F_{tot} = - k x - c \frac{dx}{dt} + F_0 \cos \omega t $$Using $F = ma$, we find the following second order differential equation for x:
$$ m \frac{d^2x}{dx^2} + c\frac{dx}{dt} + k x = F_0 \cos(\omega t) $$For this equation, there is a steady-state solution for $x(t)$ given by:
$$ x(t) = A(\omega) \cos \left[ \omega t + \theta(\omega) \right] $$where $A(\omega)$ and $\theta(\omega)$ are an amplitude and phase that depend on the driving frequency, and are given by:
\begin{eqnarray} A(\omega) & = & \frac{F_0 / m}{\sqrt{\gamma^2 \omega^2 + (\omega^2 - \omega_0^2)^2}} \\ \theta(\omega) &=& \arctan \left( \frac{\omega \gamma}{\omega^2 - \omega_0^2} \right) - \pi \end{eqnarray}where $\omega_0 = \sqrt{k/m}$ is the natural frequency of the harmonic oscillator and $\gamma = c/m$ is the damping rate.
In this problem, you will perform numerical integration of the driven damped Harmonic oscillator (DDHO) to find its steady-state response,, and compare what you find to the analytical formulas above.
We will consider the case of $m = 1$ kg, $k = 1$ N/m, $c = 0.1$ Ns/m, and $F_0 = 1$ N.
Exercise 1(a) Consider the case in which the mass is at rest and is at position $x=0$ at $t=0$, and a driving force at frequency $w = w_0$. Use the solve_ivp() routine to calculate $x(t)$ for $t$ from 0 to 200 seconds with 1000 points. Make a plot of your solution and discuss if it behaves in the way that you would expect.
m = 1
k = 1
c = 0.1
F0 = 1
w = 1
# I will chose y[0] = x and y[1] = v
def dydt(t,y):
###
### Your code here
###
T = 200
N = 1000
t = np.linspace(0,T,N)
x0 = 0
v0 = 0
# sol = solve_ivp(....)
# t = ...
# x = ...
###
### Your code here
###
# Now a plot
###
### Your code here
###
answer_12_1a_1 = np.copy(x)
question = "answer_12_1a"
num = 1
to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
res, msg = check_answer(eval(var), var)
passed = passed and res
print(msg); feedback += msg + "n"
assert passed == True, feedback
If you've done things correctly, you should see an oscillating signal that grows in amplitude and then exponentially approaches a fixed amplitude for long times (this is referred to as the "ring-up" of the oscillator). At steady state, it looks like the signal is oscillating as we expect. Qualitatively, it looks good, but does it quantitatively agree with the analytical formulas from above for the steady state behaviour?
If I want to check this, I need to extract the steady state response from the time trace above. How do I know after which time my simulated trace has reached steady-state? Technically, it only approaches steady state exponentially, so it will actually never completely reach steady-state in any case.
For estimating the steady state amplitude, we could do this by calculating the height of the last maximum of the oscillations.
Exercise 1(b): Find the steady state amplitude of the calculated time trace by finding the value of the amplitude of the last peak in the trace. For this, you can use the find_peaks() routine of scipy:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html
to find all of the peaks of the oscillations. Make a plot of the peak values vs. time, and then extract an estimate of the steady state amplitude using the value of the last peak.
Does it agree with the prediction of the analytical formulas above? To check this, create a function amp(w, m, c, k, F0) (with w as $\omega$) to calculate the amplitude from the formula above. Use this to calculate the theoretically predicted amplitude based on the parameters you have used in the simulation.
def amp(w, m, c, k, F0):
###
### Your code here
###
# peak_indices, _ = ....
# t_p = ...
# x_p = ...
###
### Your code here
###
# A plot
###
### Your code here
###
# Extract the calculated amplitude
# amp_calc = ....
###
### Your code here
###
# Calculate the theory prediction
# amp_theory = ....
###
### Your code here
###
print("Calculated steady state amplitude: ", amp_calc)
print("Predicted steady state amplitude : ", amp_theory)
answer_12_1b_1 = amp_calc
answer_12_1b_2 = amp_theory
question = "answer_12_1b"
num = 2
to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
res, msg = check_answer(eval(var), var)
passed = passed and res
print(msg); feedback += msg + "n"
assert passed == True, feedback
How can I check though that the frequency of my steady-state oscillations is correct? One way we can do this is by using the Fourier transform.
Exercise 1(c): Calculate the power spectrum of the second half of the calculated DDHO response trace calculated in 1(a). Make a plot of the power spectrum, and find the peak frequency.
In your power spectrum, keep only the part of the spectrum at positive frequencies.
# Pick out the second half of the trace to analyze for the steady state
# xss = ...
# Tss = ...
# ...
# power = ...
# f = ...
###
### Your code here
###
# Now keep only the positive frequencies
end = int(N/4)
f_pow = f[0:end]
power = power[0:end]
# Make a plot
###
### Your code here
###
# Find the frequency of the peak
# ...
# fmax = ...
# fmax_expected = ...
# ...
###
### Your code here
###
print("Calculated frequency: %.5f Hz (%.5f Rad/s)" %(f[max_index], f[max_index]*2*np.pi))
print("Theoretical frequency: %.5f Hz (%.5f Rad/s)" % (1/2/np.pi, 1))
answer_12_1c_1 = power
answer_12_1c_2 = f_pow
answer_12_1c_3 = fmax
question = "answer_12_1c"
num = 3
to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
res, msg = check_answer(eval(var), var)
passed = passed and res
print(msg); feedback += msg + "n"
assert passed == True, feedback
For the simple Harmonic oscillator, the steady state solution does not depend on the initial conditions (this is not true for non-linear restoring forces, which we will see soon!).
Exercise 1(d): Perform the calculation from Exercise 1(a) with an initial condition $x(0) = 30$. Make a plot of $x(t)$. Write a function find_ss_amp(x) to find the steady state amplitude as we did above in 1(b) by finding the peak value of the last peak. Use this function to the steady state amplitude and show that this is the same (or at least close) to the value for $x(0) = 0$.
def find_ss_amp(x):
###
### Your code here
###
x0 = 30
v0 = 0
# sol = solve_ivp(....)
# t = ...
# x = ...
###
### Your code here
###
# Make a plot
###
### Your code here
###
print("Calculated steady state amplitude: ", find_ss_amp(x))
answer_12_1d_1 = np.copy(x)
question = "answer_12_1d"
num = 1
to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
res, msg = check_answer(eval(var), var)
passed = passed and res
print(msg); feedback += msg + "n"
assert passed == True, feedback
Until now, we have been finding the solutions of differential equations given a particular initial condition: for example, for a second order differential equation, we have been specifying the initial velocity and initial position.
There are some times where you want to constrain your solution of the differential equation not by the initial conditions, but by a combination of the initial condition and the final condition. An example of this is for a second order ODE: instead of finding the solution corresponding to a given initial position $x_i$ and initial velocity $v_i$, one might want to find a solution that has a given initial position $x_i$ and a given final position $x_f$. This is an example of a boundary value problem: we are constrained by the values of the function on the boundaries, not by the initial velocity.
How does one approach such a problem?
One method for solving this is a technique known as the shooting method. In the shooting method, one makes a guess at an initial velocity $v_i$ and then calculates the final position $x_f$. If you do not get it right in the first guess, you keep redoing the calculation with different $v_i$, using a technique such as binary search to find the initial velocity $v_i$ that gives the desired $x_f$.
(The inspiration for the name of this technique comes from how one would hit a target with a cannon: a simple way is to take a trial shot, see where it lands, and then adjust the angle up and down until you get it right; or at least close enough to hit your enemy!)
As an example, we consider the following problem: I am going to throw a ball straight upwards in the air, and I want the ball to land back down at the same spot in exactly 10 seconds. What velocity do I need to give the ball?
Neglecting air resistance, this problem is easy to solve (you can solve it with high school physics!). Once we add air resistance to the problem (which we will at the end of this notebook), the problem becomes very difficult to solve analytically, which is where numerical methods become very useful.
We have learned already how to perform RK4 numerical integration of differential equations, and earlier in the course, we have also learned how to find zeros of a function efficiently using binary. With the skills you have learned, you could easily write code to perform the shooting method directly yourself (and you may be asked to in the exam!).
However, here, we will explore how to implement the shooting method using some more advanced features of the solve_ivp() routine of the scipy library. To do this, we will make use of the optional parameter events of the solve_ivp() routine. This is a special way to tell solve_ivp() to detect conditions on the integration it is performing and stop the integration when that condition is satisfied.
The wayevents= parameter works is that you need define a function that gets passed your variables y and t and which should return a number. When your "event" function undergoes a zero crossing, then you can have solve_ivp() stop the integration automatically. You can also configure if this zero crossing is a "rising" or a "falling" edge trigger (see the assignment from week 1).
How do I use this events parameter in practice? The events parameter needs to take a object (or a list of objects in case you want to track many events).
Wait a minute: what the heck is an "object"?
Objects are an advanced feature of Python. In fact, most of the things you have been working with are already objects, even though you didn't notice it! Objects are special variables that can "contain" stuff: they can contain, for example, both functions that can perform operations and data that stores values.
An example you have already seen is numpy arrays. For example, consider the following numpy array:
This array has a "field" (data stored inside the object) called shape that tells you what shape the array is:
x.shape
The object x also has functions built into it that can excute actions. For example, there is a function copy built into x that makes a copy of the array:
Cool! But how do I make my own "object"? It sounds scary! In most languages, it requires some more detailed of knowledge of the language to do so... However, python has some great shortcuts that makes building your own objects very easily, on the fly, using a technique in python with a funny name called "monkey patching".
You can construct the "object" that the solve_ivp() parameter events= requires by first creating a function that does what you need and then turn it into an "object" by adding some fields with the correct names. To be concrete, for our example of throwing a ball vertically in the air, we will want to stop the numerical integration once the height of the ball falls back down to zero. For this, I will need an event function that returns the height of the ball. If I choose my y variable array such that y[0] is the vertical position x and y[1] is the vertical velocity v, then my function would look like this:
def myevent(t,y):
return y[0]For solve_ivp(), I then need to make it an object by adding two extra fields. One field is called terminal, which tells solve_ivp() if should stop when this event is triggered, or if it should keep on going. As we are done once the ball hits the ground, we can make our even terminal:
myevent.terminal = TrueMagico! Now myevent is no longer just a function that can do stuff, but an object that contains a field that stores some information! In addition, we need to add another field to our myevent object that tells solve_ivp() if it should look for a "rising edge zero crossing" (direction = 1) or a "falling edge zero crossing" (direction = -1). We want the latter:
myevent.direction = -1When finished, the sol object that solve_ivp() returns will contain a field t_events that is a list of all the times that an event occurred at: for us, this list will contain one entry, corresponding to the time at which the ball came back to the place we threw it from.
With this, it becomes very easy to implement the shooting method:
One problem with this technique is that we need to get some reasonable values for the initial guesses. For this, it is handy to run the numerical integration a few times with different values of initial velocity to see what range we should choose for our binary search.
Exercise 2(a): Consider a ball with mass $m = 0.1$ kg thrown upwards into the air with an initial velocity of 1 m/s. Use solve_ivp() to calculate how long it takes for the ball to reach the ground again. Assume that the $g = 9.81$ m/s$^2$. For the numerical integration, you need to pick a time endpoint that will be longer than the time it will take for the ball to come back down. A time limit of 100 seconds should probably be fine.
g = 9.81
# y[0] is height (which I will call "x"), y[1] is velocity ("v")
def myevent(t,y):
###
### Your code here
###
###
### Your code here
###
def dydt(t,y):
###
### Your code here
###
vi = 1
# sol = solve_ivp(....)
###
### Your code here
###
t_f = sol.t_events[0]
print("An initial velocity of 1 m/s gave a total time of %f seconds" % t_f)
answer_12_2a_1 = t_f
question = "answer_12_2a"
num = 1
to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
res, msg = check_answer(eval(var), var)
passed = passed and res
print(msg); feedback += msg + "n"
assert passed == True, feedback
Exercise 2(b): Implement binary search (as described in the textbook) to find the initial velocity needed to have $t_f = 10$ seconds with a target accuracy of 1 ms. For your initial search range, choose $v_i$ = 1 m/s and 100 m/s. Compare to what you would expect theoretically from the acceleration due to gravity with no air resistance.
def find_tf(vi):
###
### Your code here
###
v1 = 1
t1 = find_tf(v1)
v2 = 200
t2 = find_tf(v2)
target = 1e-3
# while np.abs(t2-t1) > target:
# ...
###
### Your code here
###
# vi = ...
# vi_theory = ...
###
### Your code here
###
print(vi)
print(vi_theory)
answer_12_2b_1 = vi
answer_12_2b_2 = vi_theory
question = "answer_12_2b"
num = 2
to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
res, msg = check_answer(eval(var), var)
passed = passed and res
print(msg); feedback += msg + "n"
assert passed == True, feedback
Exercise 2(c): We will now add air resistance to our calculation. To make life simple, we will assume that the drag coefficient of our ball results in a friction force with the magnitude:
$$ |F_f| = Cv^2 $$For a 3 cm diameter sphere, Gary's estimate of the constant $C$ in this equation (based on the formulas on wikipedia) is that it has a value on the order of $C = 10^{-3}$ Ns$^2$/m$^2$. Since the force always opposes the direction of the velocity, we have to somehow account for its sign before incorporating it into an equation of motion. One way to do this is as follows:
$$ F_f = - \frac{Cv^3}{|v|} $$Repeat question 2(a) but now including air resistance in your calculation.
Before you start: do you think it will take less time or more time for the ball to fall for the same velocity?
On the one hand, when falling, the ball will move more slowly because the air resistance is holding it back. On the other hand, the ball will not go as high with air resistance.
(Give your dydt function with air resistance a different name so we can compare the two in the next question...)
your_guess = "comes back in less time /OR/ takes more time to come back (delete one!)"
g = 9.81
m = 0.1
C = 1e-3
# y[0] is height (which I will call "x"), y[1] is velocity ("v")
def myevent(t,y):
###
### Your code here
###
###
### Your code here
###
def dydt2(t,y):
###
### Your code here
###
vi = 1
# sol = solve_ivp(....)
###
### Your code here
###
t_f = sol.t_events[0]
print("An initial velocity of 1 m/s gave a total time of %f seconds" % t_f)
answer_12_2c_1 = t_f
question = "answer_12_2c"
num = 1
to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
res, msg = check_answer(eval(var), var)
passed = passed and res
print(msg); feedback += msg + "n"
assert passed == True, feedback
It seems like adding air resistance slightly decreased the amount of time it takes for the ball to come back down, but not by much.
Exercise 2(d): Make a calculation of the two times as a function of the initial velocity, with $v_i$ varying from 1 m/s to 200 m/s (720 km/h!) with 100 points. Make a plot of $t_f$ vs $v_i$ for the two cases.
Ni = 200
vi = np.linspace(1,200,Ni)
tf1 = np.empty(Ni)
tf2 = np.empty(Ni)
for i in range(Ni):
###
### Your code here
###
for i in range(Ni):
###
### Your code here
###
# Make some plots
###
### Your code here
###
answer_12_2d_1 = vi.copy()
answer_12_2d_2 = tf1.copy()
answer_12_2d_3 = tf2.copy()
question = "answer_12_2d"
num = 1
to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
res, msg = check_answer(eval(var), var)
passed = passed and res
print(msg); feedback += msg + "n"
assert passed == True, feedback
Exercise 2(e): Calculate the initial velocity required to achieve $t_f$ of 10 seconds including air resistance.
def find_tf2(vi):
###
### Your code here
###
v1 = 1
t1 = find_tf2(v1)
v2 = 200
t2 = find_tf2(v2)
target = 1e-3
# while np.abs(...) > target:
# ...
###
### Your code here
###
# vi = ...
###
### Your code here
###
print(vi)
answer_12_2e_1 = vi
question = "answer_12_2e"
num = 1
to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
res, msg = check_answer(eval(var), var)
passed = passed and res
print(msg); feedback += msg + "n"
assert passed == True, feedback
Exercise 2(f): Including air resistance, does the ball spend more time on the upwards trajectory, on the downward tragectory, or does it spend the same amount of time going up as it does going down?
Use the solve_ivp() routine to calculate the time going up and the time going down for the initial velocity you found in 2(e).
Hint: an easy way to do this is to create an additional terminal = False event function that checks when the velocity crosses zero!
# y[0] is height (which I will call "x"), y[1] is velocity ("v")
def myevent(t,y):
###
### Your code here
###
###
### Your code here
###
def myevent2(t,y):
###
### Your code here
###
###
### Your code here
###
# sol = solve_ivp(..., events=[myevent2, myevent])
###
### Your code here
###
t1 = sol.t_events[0][0]
t2 = sol.t_events[1][0]
t_down = t2-t1
print("Time up: %f" % t1)
print("Time down: %f" % (t2-t1))
answer_12_2f_1 = t1
answer_12_2f_2 = t_down
question = "answer_12_2f"
num = 2
to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
res, msg = check_answer(eval(var), var)
passed = passed and res
print(msg); feedback += msg + "n"
assert passed == True, feedback