This chapter is taken from the book A Primer on Scientific Programming with Python by H. P. Langtangen, 5th edition, Springer, 2016.
This exercise aims to solve the ODE problem \( u - 10 u^{\prime}=0 \), \( u(0)=0.2 \), for \( t\in [0,20] \).
a) Identify the mathematical function \( f(u,t) \) in the generic ODE form \( u^{\prime}=f(u,t) \).
b) Implement the \( f(u, t) \) function in a Python function.
c)
Use the ForwardEuler function from
the section Function implementation
to compute a numerical solution of the ODE problem.
Use a time step \( \Delta t =5 \).
d) Plot the numerical solution and the exact solution \( u(t)=0.2e^{0.1t} \).
e) Save the numerical solution to file. Decide upon a suitable file format.
f) Perform simulations for smaller \( \Delta t \) values and demonstrate visually that the numerical solution approaches the exact solution.
Filename: simple_ODE_func.
Solve the same ODE problem as in
Exercise 1: Solve a simple ODE with function-based code, but use
the ForwardEuler class from
the section Class implementation.
Implement the right-hand side function \( f \) as a class too.
Filename: simple_ODE_class.
Solve the same ODE problem as in
Exercise 1: Solve a simple ODE with function-based code, but use
the ForwardEuler class in the ODESolver hierarchy from
the section The ODESolver class hierarchy.
Filename: simple_ODE_class_ODESolver.
We want to make a program odesolver_cml.py which
accepts an ODE problem to be specified on the command line. The
command-line arguments are f u0 dt T, where f is the right-hand
side \( f(u,t) \) specified as a string formula,
u0 is the initial condition, dt is the
time step, and T is the final time of the simulation. A fifth
optional argument can be given to specify the name of the
numerical solution method (set any method of your choice as default value).
A curve plot of the solution versus time should be produced and stored
in a file plot.png.
Hint 1.
Use StringFunction from scitools.std for convenient conversion of
a formula on the command line to a Python function.
Hint 2.
Use the ODESolver hierarchy to solve the ODE and let the fifth
command-line argument be the class name in the ODESolver hierarchy.
Filename: odesolver_cml.
Implement the numerical method
(36)-(37).
How can you verify that the implementation is correct?
Filename: Heuns_method.
A cylindrical tank of radius \( R \) is filled with water to a height \( h_0 \). By opening a valve of radius \( r \) at the bottom of the tank, water flows out, and the height of water, \( h(t) \), decreases with time. We can derive an ODE that governs the height function \( h(t) \).
Mass conservation of water requires that the reduction in height balances the outflow. In a time interval \( \Delta t \), the height is reduced by \( \Delta h \), which corresponds to a water volume of \( \pi R^2\Delta h \). The water leaving the tank in the same interval of time equals \( \pi r^2v\Delta t \), where \( v \) is the outflow velocity. It can be shown (from what is known as Bernoulli's equation) [1] [2] that $$ \begin{equation*} v(t) = \sqrt{2gh(t) + h'(t)^2},\end{equation*} $$ where \( g \) is the acceleration of gravity. Note that \( \Delta h > 0 \) implies an increase in \( h \), which means that \( -\pi R^2\Delta h \) is the corresponding decrease in volume that must balance the outflow loss of volume \( \pi r^2 v \Delta t \). Elimination of \( v \) and taking the limit \( \Delta t\rightarrow 0 \) lead to the ODE $$ \begin{equation*} {dh\over dt} = -\left( {r\over R}\right)^2 \left( 1 - \left({r\over R}\right)^4\right)^{-1/2}\sqrt{2gh}\tp \end{equation*} $$ For practical applications \( r\ll R \) so that \( 1-(r/R)^4\approx 1 \) is a reasonable approximation, since other approximations are done as well: friction is neglected in the derivation, and we are going to solve the ODE by approximate numerical methods. The final ODE then becomes $$ \begin{equation} {dh\over dt} = -\left( {r\over R}\right)^2 \sqrt{2gh}\tp \tag{61} \end{equation} $$ The initial condition follows from the initial height of water, \( h_0 \), in the tank: \( h(0)=h_0 \).
Solve (61) by a numerical method of your choice in
a program. Set \( r = 1\hbox{ cm} \), \( R=20 \hbox{ cm} \), \( g=9.81\hbox{
m/s}^2 \), and \( h_0=1\hbox{ m} \). Use a time step of 10 seconds. Plot
the solution, and experiment to see what a proper time interval for
the simulation is. Make sure to test for \( h < 0 \) so that you do not
apply the square root function to negative numbers. Can you find an
analytical solution of the problem to compare the numerical solution
with?
Filename: tank_ODE.
Given a curve \( y=f(x) \), the length of the curve from \( x=x_0 \) to some
point \( x \) is given by the function \( s(x) \), which solves the problem
$$
\begin{equation}
{ds\over dx} = \sqrt{1 + [f'(x)]^2},\quad s(x_0)=0\tp
\tag{62}
\end{equation}
$$
Since \( s \) does not enter the right-hand side, (62) can
immediately be integrated from \( x_0 \) to
\( x \).
However, we shall solve
(62) as an ODE. Use the Forward Euler method and
compute the length of a straight line (for easy verification)
and a parabola:
\( f(x)=\frac{1}{2} x + 1 \), \( x\in [0,2] \);
\( f(x)=x^2 \),
\( x\in [0,2] \).
Filename: arclength_ODE.
A body moving vertically through a fluid (liquid or gas) is subject to three different types of forces:
Newton's second law applied to the body says that the sum of the listed forces must equal the mass of the body times its acceleration \( a \): $$ \begin{equation*} F_g + F_d + F_b = ma,\end{equation*} $$ which gives $$ \begin{equation*} -mg -\frac{1}{2}C_D\varrho A|v|v + \varrho gV = ma\tp\end{equation*} $$ The unknowns here are \( v \) and \( a \), i.e., we have two unknowns but only one equation. From kinematics in physics we know that the acceleration is the time derivative of the velocity: \( a = dv/dt \). This is our second equation. We can easily eliminate \( a \) and get a single differential equation for \( v \): $$ \begin{equation*} -mg -\frac{1}{2}C_D\varrho A|v|v + \varrho gV = m{dv\over dt}\tp \end{equation*} $$ A small rewrite of this equation is handy: we express \( m \) as \( \varrho_bV \), where \( \varrho_b \) is the density of the body, and we isolate \( dv/dt \) on the left-hand side, $$ \begin{equation} {dv\over dt} = -g\left(1 - {\varrho\over\varrho_b}\right) -\frac{1}{2}C_D{\varrho A\over\varrho_b V}|v|v\tp \tag{63} \end{equation} $$ This differential equation must be accompanied by an initial condition: \( v(0)=V_0 \).
a) Make a program for solving (63) numerically, using any numerical method of your choice.
Hint.
It is not strictly necessary, but it is an elegant Python solution to
implement the right-hand side of (63) in
the __call__ method of a class where the parameters
\( g \), \( \varrho \), \( \varrho_b \), \( C_D \), \( A \), and \( V \) are data attributes.
b) To verify the program, assume a heavy body in air such that the \( F_b \) force can be neglected, and assume a small velocity such that the air resistance \( F_d \) can also be neglected. Mathematically, setting \( \varrho =0 \) removes both these terms from the equation. The solution is then \( v(t)=y'(t)=v_0 - gt \). Observe through experiments that the linear solution is exactly reproduced by the numerical solution regardless of the value of \( \Delta t \). (Note that if you use the Forward Euler method, the method can become unstable for large \( \Delta t \), see the section Example: Exponential decay. and time steps above the critical limit for stability cannot be used to reproduce the linear solution.) Write a test function for automatically checking that the numerical solution is \( u_k=v_0-gk\Delta t \) in this test case.
c) Make a function for plotting the forces \( F_g \), \( F_b \), and \( F_d \) as functions of \( t \). Seeing the relative importance of the forces as time develops gives an increased understanding of how the different forces contribute to changing the velocity.
d) Simulate a skydiver in free fall before the parachute opens. We set the density of the human body as \( \varrho_b = 1003 \hbox{ kg}/\hbox{m}^3 \) and the mass as \( m=80 \) kg, implying \( V=m/\varrho_b = 0.08\hbox{ m}^3 \). We can base the cross-sectional area \( A \) the assumption of a circular cross section of diameter 50 cm, giving \( A= \pi R^2 = 0.9\hbox{ m}^2 \). The density of air decreases with height, and we here use the value 0.79 \( \hbox{kg/m}^3 \) which is relevant for about 5000 m height. The \( C_D \) coefficient can be set as 0.6. Start with \( v_0=0 \).
e) A ball with the size of a soccer ball is placed in deep water, and we seek to compute its motion upwards. Contrary to the former example, where the buoyancy force \( F_b \) is very small, \( F_b \) is now the driving force, and the gravity force \( F_g \) is small. Set \( A=\pi a^2 \) with \( a=11 \) cm. The mass of the ball is 0.43 kg, the density of water is 1000 \( \hbox{kg/m}^3 \), and \( C_D \) can be taken as 0.4. Start with \( v_0=0 \).
Filename: body_in_fluid.
The solution of (63) often tends to a constant velocity, called the terminal velocity. This happens when the sum of the forces, i.e., the right-hand side in (63), vanishes.
a) Compute the formula for the terminal velocity by hand.
b)
Solve the ODE using class ODESolver
and call the solve method with a terminate function that
terminates the computations when a constant velocity is reached, that
is, when \( |v(t_n)-v(t_{n-1})|\leq\epsilon \), where \( \epsilon \) is a
small number.
c) Run a series of \( \Delta t \) values and make a graph of the terminal velocity as a function of \( \Delta t \) for the two cases in Exercise 8: Simulate a falling or rising body in a fluid d) and e). Indicate the exact terminal velocity in the plot by a horizontal line.
Filename: body_in_fluid_termvel.
Consider the logistic model (5): $$ \begin{equation*} u^{\prime}(t)=\alpha u(t)\left( 1-\frac{u(t)}{R}\right),\quad u(0)=U_0\tp\end{equation*} $$ This problem involves three input parameters: \( U_0 \), \( R \), and \( \alpha \). Learning how \( u \) varies with \( U_0 \), \( R \), and \( \alpha \) requires much experimentation where we vary all three parameters and observe the solution. A much more effective approach is to scale the problem. By this technique the solution depends only on one parameter: \( U_0/R \). This exercise tells how the scaling is done.
The idea of scaling is to introduce dimensionless versions of the independent and dependent variables: $$ \begin{equation*} v = \frac{u}{u_c},\quad \tau = \frac{t}{t_c},\end{equation*} $$ where \( u_c \) and \( t_c \) are characteristic sizes of \( u \) and \( t \), respectively, such that the dimensionless variables \( v \) and \( \tau \) are of approximately unit size. Since we know that \( u \rightarrow R \) as \( t\rightarrow\infty \), \( R \) can be taken as the characteristic size of \( u \).
Insert \( u=Rv \) and \( t=t_c\tau \) in the governing ODE and choose \( t_c = 1/\alpha \). Show that the ODE for the new function \( v(\tau) \) becomes $$ \begin{equation} {dv\over d\tau} = v(1-v),\quad v(0)=v_0\tp \tag{64} \end{equation} $$ We see that the three parameters \( U_0 \), \( R \), and \( \alpha \) have disappeared from the ODE problem, and only one parameter \( v_0=U_0/R \) is involved.
Show that if \( v(\tau) \) is computed, one can recover \( u(t) \) by $$ \begin{equation} u(t) = Rv(\alpha t)\tp \tag{65} \end{equation} $$ Geometrically, the transformation from \( v \) to \( u \) is just a stretching of the two axis in the coordinate system.
Make a program logistic_scaled.py where you compute
\( v(\tau) \), given \( v_0=0.05 \), and then you use (65)
to plot \( u(t) \) for \( R=100,500,1000 \) and \( \alpha=1 \) in one figure,
and \( u(t) \) for \( \alpha=1, 5, 10 \) and \( R=1000 \) in another figure.
Note how effectively you can generate \( u(t) \) without needing to solve
an ODE problem, and also note how varying \( R \) and \( \alpha \) impacts
the graph of \( u(t) \).
Filename: logistic_scaled.
Use classes Problem2 and AutoSolver
from the section Example: The logistic equation with problem and solver classes to study logistic growth
when the carrying capacity of
the environment, \( R \), changes periodically with time:
\( R=500 \) for \( it_s\leq t < (i+1)t_s \) and \( R=200 \) for
\( (i+1)t_s\leq t < (i+2)t_s \), with \( i=0, 2, 4, 6, \ldots \).
Use the same data as in the section Example: The logistic equation with problem and solver classes, and
find some relevant sizes of the period of variation, \( t_s \), to
experiment with.
Filename: seasonal_logistic_growth.
Newton's law of cooling, $$ \begin{equation} {dT\over dt} = -h(T-T_s) \tag{66} \end{equation} $$ can be used to see how the temperature \( T \) of an object changes because of heat exchange with the surroundings, which have a temperature \( T_s \). The parameter \( h \), with unit \( \hbox{s}^{-1} \) is an experimental constant (heat transfer coefficient) telling how efficient the heat exchange with the surroundings is. For example, (66) may model the cooling of a hot pizza taken out of the oven. The problem with applying (66) is that \( h \) must be measured. Suppose we have measured \( T \) at \( t=0 \) and \( t_1 \). We can use a rough Forward Euler approximation of (66) with one time step of length \( t_1 \), $$ \begin{equation*} \frac{T(t_1)-T(0)}{t_1} = -h(T(0)-T_s),\end{equation*} $$ to make the estimate $$ \begin{equation} h = \frac{T(t_1)-T(0)}{t_1(T_s-T(0))}\tp \tag{67} \end{equation} $$
a) The temperature of a pizza is 200 C at \( t=0 \), when it is taken out of the oven, and 180 C after 50 seconds, in a room with temperature 20 C. Find an estimate of \( h \) from the formula above.
b) Solve (66) numerically by a method of your choice to find the evolution of the temperature of the pizza. Plot the solution.
Hint.
You may solve the ODE the way you like, but the solve method in
the classes in the ODESolver hierarchy
accepts an optional terminate function that can be used to
terminate the solution process when \( T \) is sufficiently close to \( T_s \).
Reading the first part of the section Example: The logistic equation with problem and solver classes may be useful.
Filename: pizza_cooling1.
We address the same physical problem as in Exercise 12: Solve an ODE until constant solution, but we will now provide a class-based implementation for the user's code.
a)
Make a class Problem containing the parameters \( h \), \( T_s \), \( T(0) \),
and \( \Delta t \) as data attributes. Let these parameters be set in the
constructor.
The right-hand side of the ODE can be
implemented in a __call__ method. If you use a class from the
ODESolver hierarchy to solve the ODE, include the terminate
function as a method in class Problem.
Create a stand-alone function estimate_h(t0, t1, T0, T1)
which applies (67) from Exercise 12: Solve an ODE until constant solution
to estimate the \( h \) parameter based on the initial temperature and
one data point \( (t_1,T(t_1)) \). You can use this function to estimate
a value for \( h \) prior to calling the constructor in the Problem class.
Hint.
You may want to read the section Example: The logistic equation with problem and solver classes to see why and how
a class Problem is constructed.
b)
Implement a test function test_Problem() for testing that class
Problem works. (It is up to you to define how to test the class.)
c)
What are the advantages and disadvantages with class Problem
compared to using plain functions (in your view)?
d) We now want to run experiments with different values of some parameters: \( T_s=15, 22, 30 \) C and \( T(0)=250, 200 \) C. For each \( T(0) \), plot \( T \) for the three \( T_s \) values. The estimated value of \( h \) in Exercise 12: Solve an ODE until constant solution can be reused here.
Hint.
The typical elegant Python way to solve such a problem goes as follows.
Write a function solve(problem) that takes a Problem object
with name problem as argument and performs what it takes to
solve one case (i.e., solve must solve the ODE and plot the solution).
A dictionary can for each \( T_0 \) value hold a list of the cases to
be plotted together. Then we loop over the problems in the
dictionary of lists and call solve for each problem:
# Given h and dt
problems = {T_0: [Problem(h, T_s, T_0, dt)
for T_s in 15, 22, 30] for T_0 in 250, 200}
for T_0 in problems:
hold('off')
for problem in problems[T_0]:
solve(problem)
hold('on')
savefig('T0_%g'.pdf % T_0)
Problem instances becomes more involved.
Filename: pizza_cooling2.
Use the scaling approach outlined in Exercise 10: Scale the logistic equation to
"scale away" the parameters in the ODE in Exercise 12: Solve an ODE until constant solution. That is, introduce a new unknown
\( u=(T-T_s)/(T(0)-T_s) \) and a new time scale \( \tau = th \). Find the ODE
and the initial condition that governs the \( u(\tau ) \) function. Make
a program that computes \( u(\tau ) \) until \( |u| < 0.001 \). Store the
discrete \( u \) and \( \tau \) values in a file u_tau.dat if that file is
not already present (you can use os.path.isfile(f) to test if a file
with name f exists). Create a function T(u, tau, h, T0, Ts) that
loads the \( u \) and \( \tau \) data from the u_tau.dat file and returns
two arrays with \( T \) and \( t \) values, corresponding to the computed
arrays for \( u \) and \( \tau \). Plot \( T \) versus \( t \). Give the parameters
\( h \), \( T_s \), and \( T(0) \) on the command line. Note that this program is
supposed to solve the ODE once and then recover any \( T(t) \) solution by
a simple scaling of the single \( u(\tau ) \) solution.
Filename: pizza_cooling3.
The ForwardEuler_func.py file is not well organized to be used as a module, say for doing
>>> from ForwardEuler_func import ForwardEuler
>>> u, t = ForwardEuler(lambda u, t: -u**2, U0=1, T=5, n=30)
import statement will execute
a main program in the ForwardEuler_func.py file,
involving plotting of the solution in an example.
Also, the verification tests are run, which in more complicated problems
could take considerable time and thus make the import statement
hang until the tests are done.
Go through the file and modify it such that it becomes a module where
no code is executed unless the module file is run as a program.
Filename: ForwardEuler_func2.
The equation \( u^{\prime} = -au \) is a relevant model for radioactive decay, where \( u(t) \) is the fraction of particles that remains in the radioactive substance at time \( t \). The parameter \( a \) is the inverse of the so-called mean lifetime of the substance. The initial condition is \( u(0)=1 \).
a)
Introduce a class Decay to hold information about the physical
problem: the parameter \( a \) and a __call__ method for computing
the right-hand side \( -au \) of the ODE.
b)
Initialize an instance of class
Decay with \( a =\ln (2)/5600 \) 1/y. The unit 1/y means one divided
by year, so time is here measured in years, and the particular value
of \( a \) corresponds to the Carbon-14 radioactive isotope whose decay is
used extensively in dating organic material that is tens of thousands
of years old.
c) Solve \( u^\prime = -au \) with a time step of 500 y, and simulate the radioactive decay for \( T=20,000 \) y. Plot the solution. Write out the final \( u(T) \) value and compare it with the exact value \( e^{-aT} \).
Filename: radioactive_decay.
The inverse function \( g \) of some function \( f(x) \) takes the value of \( f(x) \) back to \( x \) again: \( g(f(x))=x \). The common technique to compute inverse functions is to set \( y=f(x) \) and solve with respect to \( x \). The formula on the right-hand side is then the desired inverse function \( g(y) \).
We can formulate a general procedure for computing inverse functions from an ODE problem. If we differentiate \( y=f(x) \) with respect to \( y \), we get \( 1 = f'(x)\frac{dx}{dy} \) by the chain rule. The inverse function we seek is \( x(y) \), but this function then fulfills the ODE $$ \begin{equation} x'(y) = \frac{1}{f'(x)}\tp \tag{68} \end{equation} $$ That \( y \) is the independent coordinate and \( x \) a function of \( y \) can be a somewhat confusing notation, so we might introduce \( u \) for \( x \) and \( t \) for \( y \): $$ \begin{equation*} u^{\prime}(t) = \frac{1}{f'(u)}\tp\end{equation*} $$ The initial condition is \( x(0)=x_r \) where \( x_r \) solves the equation \( f(x_r)=0 \) (\( x(0) \) implies \( y=0 \) and then from \( y=f(x) \) it follows that \( f(x(0))=0 \)).
Make a program that can use the described method to
compute the inverse function of \( f(x) \), given \( x_r \).
Use any numerical method of your choice for solving the ODE problem.
Verify the implementation for \( f(x)=2x \).
Apply the method for \( f(x)=\sqrt{x} \) and
plot \( f(x) \) together with its inverse function.
Filename: inverse_ODE.
The method for computing inverse functions described in
Exercise 17: Compute inverse functions by solving an ODE is very general. The purpose now is
to use this general approach
to make a more reusable utility, here called Inverse, for
computing the inverse of some Python function f(x) on
some interval I=[a,b].
The utility can be used as follows to calculate the inverse
of \( \sin x \) on \( I=[0, \pi/2] \):
def f(x):
return sin(x)
# Compute the inverse of f
inverse = Inverse(f, x0=0, I=[0, pi/2], resolution=100)
x, y = Inverse.compute()
plot(y, x, 'r-',
x, f(x), 'b-',
y, asin(y), 'go')
legend(['computed inverse', 'f(x)', 'exact inverse'])
x0 is the value of x at 0, or in general
at the left point of the interval: I[0].
The parameter resolution tells how many equally sized intervals
\( \Delta y \) we use in the numerical integration of the ODE.
A default choice of 1000 can be used if it is not given by the user.
Write class Inverse and put it in a module. Include a test function
test_Inverse() in the module for verifying that class Inverse
reproduces the exact solution in the test problem \( f(x)=2x \).
Filename: Inverse1.
Extend the module in Exercise 18: Make a class for computing inverse functions such that
the value of \( x(0) \) (x0 in class Inverse's constructor)
does not need to be provided by the user.
Hint.
Class Inverse can compute a value of \( x(0) \) as the root
of \( f(x)=0 \).
You may use the
Bisection method from the section ref{sec:input:summarizingex},
Newton's method from the section ref{sec:diffeq:Newtonsmethod:sec},
or the Secant method from ref{sec:diffeq:ex14b} to
solve \( f(x)=0 \). Class Inverse should figure out a suitable
initial interval for the Bisection method or start values for the
Newton or Secant methods. Computing \( f(x) \) for \( x \) at many
points and examining these may help in solving \( f(x)=0 \) without
any input from the user.
Filename: Inverse2.
Instead of solving an ODE for computing the inverse function \( g(y) \)
of some function \( f(x) \), as explained in Exercise 17: Compute inverse functions by solving an ODE,
one may use a simpler approach based on ideas from
the section From discrete to continuous solution.
Say we compute discrete values of \( x \) and \( f(x) \), stored in
the arrays x and y. Doing a plot(x, y)
shows \( y=f(x) \) as a function of \( x \), and doing
plot(y, x) shows \( x \) as a function of \( y \), i.e., we can
trivially plot the inverse function \( g(y) \) (!).
However, if we want the inverse function of \( f(x) \) as some
Python function g(y) that we can call for any y,
we can use the tool wrap2callable from
the section From discrete to continuous solution to turn the discrete
inverse function, described by the arrays y (independent
coordinate) and x (dependent coordinate), into a
continuous function g(y):
from scitools.std import wrap2callable
g = wrap2callable((y, x))
y = 0.5
print g(y)
g(y) function applies linear interpolation in each
interval between the points in the y array.
Implement this method in a program.
Verify the implementation for \( f(x)=2x \), \( x\in [0,4] \), and apply
the method to \( f(x)=\sin x \) for \( x\in [0, \pi/2] \).
Filename: inverse_wrap2callable.
Use the file ForwardEuler_func.py from the section Function implementation as starting point for implementing the famous
and widely used 4th-order Runge-Kutta method
(41)-(45). Use the test
function involving a linear \( u(t) \) for verifying the
implementation. Exercise 23: Compare ODE methods suggests an application
of the code.
Filename: RK4_func.
Carry out the steps in Exercise 21: Code the 4th-order Runge-Kutta method; function, but base
the implementation on the file ForwardEuler.py from the section Class implementation.
Filename: RK4_class.
Investigate the accuracy of the 4th-order Runge-Kutta method and the
Forward Euler scheme for solving the (challenging)
ODE problem
$$
\begin{equation}
{dy\over dx} = {1\over 2(y-1)},\quad y(0)=1+\sqrt{\epsilon},\quad x\in [0,4],
\tag{69}
\end{equation}
$$
where \( \epsilon \) is a small number, say \( \epsilon = 0.001 \).
Start with four steps in \( [0,4] \) and reduce the step size repeatedly by
a factor of two until you find the solutions sufficiently accurate.
Make a plot of the numerical solutions along with
the exact solution \( y(x) = 1 + \sqrt{x+\epsilon} \) for each step size.
Filename: yx_ODE_FE_vs_RK4.
The ForwardEuler_func.py
file from the section Function implementation does not contain any test
function for verifying the implementation. We can use the fact that
linear functions of time will be exactly reproduced by most numerical
methods for ODEs. A simple system of two ODEs with linear solutions
\( v(t)=2+3t \) and \( w(t)=3+4t \) is
$$
\begin{align}
v' &= 3 + (3 + 4t - w)^3,
\tag{70} \\
w' &= 4 + (2 + 3t - v)^4
\tag{71}
\end{align}
$$
Write a test function test_ForwardEuler() for comparing the
numerical solution of this system with the exact solution.
Filename: ForwardEuler_sys_func2.
Use the file ForwardEuler_sys_func.py from the section Function implementation as starting point for implementing
Heun's method (36)-(37)
for systems of ODEs.
Verify the solution using the test function suggested in
Exercise 24: Code a test function for systems of ODEs.
Filename: Heun_sys_func.
Carry out the steps in Exercise 25: Code Heun's method for ODE systems; function, but
make a class implementation based on the file ForwardEuler_sys.py from the section Class implementation.
Filename: Heun_sys_class.
a)
Implement the Leapfrog method specified in formula
(35) from
the section Numerical methods in a subclass of
ODESolver. Place the code in a separate module file Leapfrog.py.
b) Make a test function for verifying the implementation.
Hint. Use the fact that the method will exactly produce a linear \( u \).
c) Make a movie that shows how the Leapfrog method, the Forward Euler method, and the 4th-order Runge-Kutta method converge to the exact solution as \( \Delta t \) is reduced. Use the model problem \( u^{\prime}=u \), \( u(0)=1 \), \( t\in [0,8] \), with \( n=2^k \) intervals, \( k=1,2\ldots,14 \). Place the movie generation in a function.
d) Repeat c) for the model problem \( u^{\prime}=-u \), \( u(0)=1 \), \( t\in [0,5] \), with \( n=2^k \) intervals, \( k=1,2\ldots,14 \). In the movie, start with the finest resolution and reduce \( n \) until \( n=2 \). The lessons learned is that Leapfrog can give completely wrong, oscillating solutions if the time step is not small enough.
Filename: Leapfrog.
Do Exercise 27: Implement and test the Leapfrog method with the
3rd-order Adams-Bashforth method (46).
Filename: AdamBashforth3.
Consider two radioactive substances A and B. The nuclei in substance A decay to form nuclei of type B with a mean lifetime \( \tau_A \), while substance B decay to form type A nuclei with a mean lifetime \( \tau_B \). Letting \( u_A \) and \( u_B \) be the fractions of the initial amount of material in substance A and B, respectively, the following system of ODEs governs the evolution of \( u_A(t) \) and \( u_B(t) \): $$ \begin{align} u_A' &= u_B/\tau_B - u_A/\tau_A, \tag{72}\\ u_B' &= u_A/\tau_A - u_B/\tau_B, \tag{73} \end{align} $$ with \( u_A(0)=u_B(0)=1 \).
a)
Introduce a problem class, which
holds the parameters \( \tau_A \) and \( \tau_B \) and offers
a __call__ method to compute the right-hand side
vector of the ODE system, i.e.,
\( (u_B/\tau_B - u_A/\tau_A,u_A/\tau_A - u_B/\tau_B) \).
b)
Solve for \( u_A \) and \( u_B \) using a subclass in the ODESolver hierarchy
and the parameter choices \( \tau_A=8 \) minutes, \( \tau_B=40 \) minutes, and
\( \Delta t = 10 \) seconds.
c) Plot \( u_A \) and \( u_B \) against time measured in minutes.
d) From the ODE system it follows that the ratio \( u_A/u_B\rightarrow \tau_A/\tau_B \) as \( t\rightarrow\infty \) (assuming \( u_A'=u_B'=0 \) in the limit \( t\rightarrow\infty \)). Extend the problem class with a test method for checking that two given solutions \( u_A \) and \( u_B \) fulfill this requirement. Verify that this is indeed the case with the computed solutions in b).
Filename: radioactive_decay2.
Implement the 2nd-order Runge-Kutta method specified in formula
(38). Use a plain function RungeKutta2 of the
type shown in the section The Forward Euler scheme for the Forward
Euler method. Construct a test problem where you know the analytical
solution, and plot the difference between the numerical and analytical
solution. Demonstrate that the numerical solution approaches the
exact solution as \( \Delta t \) is reduced.
Filename: RungeKutta2_func.
a)
Make a new subclass RungeKutta2 in the ODESolver hierarchy from
the section The ODESolver class hierarchy for solving ordinary differential
equations with the 2nd-order Runge-Kutta method specified in formula
(38).
b) Construct a test problem where you can find an exact solution. Run different values of \( \Delta t \) and demonstrate in a plot that the numerical solution approaches the exact solution as \( \Delta t \) is decreased. Put the code that creates the plot in a function.
c)
Make a test function test_RungeKutta2_against_hand_calc()
where you do the computations of \( u_1 \) and
\( u_2 \) i Python based on the mathematical formulas. Thereafter, run the
RungeKutta2 class for two time steps and check that the two solutions
are equal (within a small tolerance).
Use an ODE where the right-hand side depends on \( t \) as well as \( u \) such
that you can test that RungeKutta2 treats the \( t \) argument
in \( f(u,t) \) correctly.
d)
Make a module
out of the RungeKutta2 class and the associated functions.
Call the functions from a test block in the module file.
Filename: RungeKutta2.
a) Implement the numerical method (47)-(48) as a function
iterated_Midpoint_method(f, U0, T, n, N)
f is a Python implementation of \( f(u,t) \), U0 is the initial
condition \( u(0)=U_0 \), T is the final time of the simulation, n is
the number of time steps, and N is the parameter \( N \) in the method
(47). The iterated_Midpoint_method should
return two arrays: \( u_0,\ldots,u_n \) and \( t_0,\ldots,t_n \).
Hint. You may want to build the function on the software described in the section Function implementation.
b)
To verify the implementation, calculate by hand \( u_1 \) and \( u_2 \) when \( N=2 \) for
the ODE \( u^{\prime}=-2u \), \( u(0)=1 \), with \( \Delta t = 1/4 \). Compare
your hand calculations with the results of the program. Make a
test function test_iterated_Midpoint_method() for automatically
comparing the hand calculations with the output of the function in a).
c) Consider the ODE problem \( u^{\prime}=-2(t-4)u \), \( u(0)=e^{-16} \), \( t\in [0,8] \), with exact solution \( u=e^{-(t-4)^2} \). Write a function for comparing the numerical and exact solution in a plot. Enable setting of \( \Delta t \) and \( N \) from the command line and use the function to study the behavior of the numerical solution as you vary \( \Delta t \) and \( N \). Start with \( \Delta t=0.5 \) and \( N=1 \). Continue with reducing \( \Delta t \) and increasing \( N \).
Filename: MidpointIter_func.
The purpose of this exercise is to implement the numerical method
(47)-(48) in a class
MidpointIter, like the ForwardEuler class from the section Class implementation. Also make a test function
test_MidpointIter() where you apply the verification technique from
Exercise 32: Code the iterated midpoint method; functionb.
Filename: MidpointIter_class.
Implement the numerical method
(47)-(48) in a subclass in
the ODESolver hierarchy. The code should reside in a separate file
where the ODESolver class is imported. One can either fix \( N \) or
introduce an \( \epsilon \) and iterate until the change in
\( |v_q-v_{q-1}| \) is less than \( \epsilon \). Allow the constructor to take
both \( N \) and \( \epsilon \) as arguments. Compute a new \( v_q \) as long as
\( q\leq N \) and \( |v_q-v_{q-1}|>\epsilon \). Let \( N=20 \) and \( \epsilon =
10^{-6} \) by default. Store \( N \) as an attribute such that the user's
code can access what \( N \) was in the last computation. Also write a
test function for verifying the implementation.
Filename: MidpointIter.
We want to see how various numerical methods treat the following
ODE problem:
$$ u^\prime = -2(t-4)u,\quad u(0)=e^{-16},\quad t\in (0,10]\tp$$
The exact solution is a Gaussian function: \( u(t)=e^{-(t-4)^2} \).
Compare the Forward Euler method with other methods of
your choice in the same plot. Relevant methods are
the 2nd-order Runge-Kutta method, Heun's method, the
4th-order Runge-Kutta method, or an implicit method.
Put the value of \( \Delta t \) in the title of the plot.
Perform experiments with \( \Delta t = 0.3, 0.25, 0.1, 0.05, 0.01, 0.001 \)
and report how the various methods behave.
Filename: methods4gaussian.
Make a movie for illustrating how three selected numerical methods
converge to the exact solution for the problem described in Exercise 35: Compare the accuracy of various methods for ODEs as \( \Delta t \) is reduced. Start with
\( \Delta t=1 \), fix the \( y \) axis in \( [-0.1, 1.1] \), and reduce \( \Delta t \)
by a quite small factor, say 1.5, between each frame in the movie. The
movie must last until all methods have their
curves visually on top of the exact solution.
Filename: animate_methods4gaussian.
The approximation error when solving an ODE numerically is usually of the form \( C\Delta t^r \), where \( C \) and \( r \) are constants that can be estimated from numerical experiments. The constant \( r \), called the convergence rate, is of particular interest. Halving \( \Delta t \) halves the error if \( r=1 \), but if \( r=3 \), halving \( \Delta t \) reduces the error by a factor of 8.
The exercise "Compute convergence rates of numerical integration methods" in the document Object-oriented programming [3] describes a method for estimating \( r \) from two consecutive experiments. Make a function
ODE_convergence(f, U0, u_e, method, dt=[])
dt list. The argument f is a Python implementation
of \( f(u,t) \) in the ODE \( u^{\prime}=f(u,t) \).
The initial condition is \( u(0)=U_0 \), where \( U_0 \) is given
as the U0 argument, u_e is the exact solution \( u_e(t) \)
of the ODE, and method is the name of a class in the ODESolver
hierarchy. The error between the exact solution \( u_e \) and the computed
solution \( u_0,u_1,\ldots,u_n \) can be defined as
$$
\begin{equation*} e = \left(\Delta t\sum_{i=0}^n (u_e(t_i) - u_i)^2\right)^{1/2}\tp\end{equation*}
$$
Call the ODE_convergence function for some numerical methods
and print the estimated \( r \)
values for each method. Make your own choice of the ODE
problem and the collection of numerical methods.
Filename: ODE_convergence.
In Exercise 8: Simulate a falling or rising body in a fluid we compute the velocity
\( v(t) \). The position of the
body, \( y(t) \), is related to the velocity by \( y'(t)=v(t) \).
Extend the program from Exercise 8: Simulate a falling or rising body in a fluid to solve
the system
$$
\begin{align*}
{dy\over dt} &= v,\\
{dv\over dt} &= -g\left(1 - {\varrho\over\varrho_b}\right) -
-\frac{1}{2}C_D{\varrho A\over\varrho_b V}|v|v\tp
\end{align*}
$$
Filename: body_in_fluid2.
The differential equations governing the horizontal and vertical motion of a ball subject to gravity and air resistance read $$ \begin{align} {d^2 x\over dt^2} = -\frac{3}{8}C_D\bar\varrho a^{-1}\sqrt{\left(\frac{dx}{dt}\right)^2 + \left({dy\over dt}\right)^2} \frac{dx}{dt}, \tag{74}\\ {d^2 y\over dt^2} = -g - \frac{3}{8}C_D\bar\varrho a^{-1}\sqrt{\left(\frac{dx}{dt}\right)^2 + \left({dy\over dt}\right)^2} {dy\over dt}, \tag{75} \end{align} $$ where \( (x,y) \) is the position of the ball (\( x \) is a horizontal measure and \( y \) is a vertical measure), \( g \) is the acceleration of gravity, \( C_D=0.2 \) is a drag coefficient, \( \bar\varrho \) is the ratio of the density of air and the ball, and \( a \) is the radius of the ball.
Let the initial condition be \( x=y=0 \) (start position in the origin) and $$ \begin{equation*}dx/dt = v_0\cos\theta, \quad dy/dt = v_0\sin\theta,\end{equation*} $$ where \( v_0 \) is the magnitude of the initial velocity and \( \theta \) is the angle the velocity makes with the horizontal.
a) Express the two second-order equations above as a system of four first-order equations with four initial conditions.
b)
Implement the right-hand side in a problem class
where the physical parameters \( C_D \), \( \bar\varrho \), \( a \), \( v_0 \), and
\( \theta \) are stored along with the initial conditions.
You may also want to add a terminate method in this class for
checking when the ball hits the ground and then terminate the
solution process.
c) Simulate a hard football kick where \( v_0=120 \) km/h and \( \theta \) is 30 degrees. Take the density of the ball as 0.017 \( \hbox{hg}/\hbox{m}^3 \) and the radius as 11 cm. Solve the ODE system for \( C_D=0 \) (no air resistance) and \( C_D=0.2 \), and plot \( y \) as a function of \( x \) in both cases to illustrate the effect of air resistance. Make sure you express all units in kg, m, s, and radians.
Filename: kick2D.
An electric circuit with a resistor, a capacitor, an inductor,
and a voltage source can be described by
the ODE
$$
\begin{equation}
L{dI\over dt} + RI + {Q\over C} = E(t),
\tag{76}
\end{equation}
$$
where \( LdI/dt \) is the voltage drop across the inductor,
\( I \) is the current (measured in amperes, A),
\( L \) is the inductance (measured in henrys, H),
\( R \) is the resistance (measured in ohms, \( \Omega \)),
\( Q \) is the charge on the
capacitor (measured in coulombs, C),
\( C \) is the capacitance (measured in farads, F), \( E(t) \) is the
time-variable voltage source (measured in volts, V),
and \( t \) is time (measured in seconds, s).
There is a relation between
\( I \) and \( Q \):
$$
\begin{equation}
{dQ\over dt} = I\tp
\tag{77}
\end{equation}
$$
Equations (76)-(77)
is a system two ODEs.
Solve these for \( L=1 \) H,
\( E(t)=2\sin\omega t \) V, \( \omega^2=3.5\hbox{ s}^{-2} \),
\( C=0.25 \) C, \( R=0.2\ \Omega \), \( I(0)=1 \) A, and \( Q(0)=1 C \).
Use the Forward Euler scheme with \( \Delta t = 2\pi/(60\omega) \).
The solution will, after some time, oscillate with the same period as
\( E(t) \), a period of \( 2\pi/\omega \). Simulate 10 periods.
Filename: electric_circuit.
It turns out that the Forward Euler scheme overestimates the amplitudes of the oscillations. The more accurate 4th-order Runge-Kutta method is much better for this type of differential equation model.
We shall in this exercise model epidemiological diseases such as measles or swine flu. Suppose we have three categories of people: susceptibles (S) who can get the disease, infected (I) who have developed the disease and who can infect susceptibles, and recovered (R) who have recovered from the disease and become immune. Let \( S(t) \), \( I(t) \), and \( R(t) \) be the number of people in category S, I, and R, respectively. We have that \( S+I+R=N \), where \( N \) is the size of the population, assumed constant here for simplicity.
When people mix in the population there are \( SI \) possible pairs of susceptibles and infected, and a certain fraction \( \beta SI \) per time interval meets with the result that the infected "successfully" infect the susceptible. During a time interval \( \Delta t \), \( \beta SI\Delta t \) get infected and move from the S to the I category: $$ \begin{equation*} S(t+\Delta t) = S(t) - \beta SI\Delta t\tp\end{equation*} $$ We divide by \( \Delta t \) and let \( \Delta \rightarrow 0 \) to get the differential equation $$ \begin{equation} S'(t) = -\beta SI\tp \tag{78} \end{equation} $$ A fraction \( \nu I \) of the infected will per time unit recover from the disease. In a time \( \Delta t \), \( \nu I\Delta t \) recover and move from the I to the R category. The quantity \( 1/\nu \) typically reflects the duration of the disease. In the same time interval, \( \beta SI\Delta t \) come from the S to the I category. The accounting for the I category therefore becomes $$ \begin{equation*} I(t+\Delta t) = I(t) + \beta SI\Delta t - \nu I \Delta t,\end{equation*} $$ which in the limit \( \Delta t\rightarrow 0 \) becomes the differential equation $$ \begin{equation} I'(t) = \beta SI - \nu I\tp \tag{79} \end{equation} $$ Finally, the R category gets contributions from the I category: $$ \begin{equation*} R(t+\Delta t) = R(t) + \nu I\Delta t\tp\end{equation*} $$ The corresponding ODE for \( R \) reads $$ \begin{equation} R'(t) = \nu I\tp \tag{80} \end{equation} $$ In case the recovered do not become immune, we do not need the recovered category, since the recovered go directly out of the I category to the S category again. This gives a contribution \( \nu I \) to the equation for \( S \) and we end up with a system for \( S \) and \( I \).
The system (78)-(80) is known as a SIR model in epidemiology (which is the name of the scientific field studying the spreading of epidemic diseases).
Make a function for solving the differential equations in the SIR model by any numerical method of your choice. Make a separate function for visualizing \( S(t) \), \( I(t) \), and \( R(t) \) in the same plot.
Adding the equations shows that
\( S'+I'+R'=0 \), which means that \( S+I+R \) must be constant. Perform
a test at each time level for checking that \( S+I+R \) equals
\( S_0+I_0+R_0 \) within some small tolerance.
If a subclass of ODESolver
is used to solve the ODE system, the test can be implemented as a
user-specified terminate function that is called by the solve
method a every time level (simply return True for termination if
\( S+I+R \) is not sufficiently constant).
A specific population has 1500 susceptibles and one infected. We are
interested in how the disease develops.
Set \( S(0)=1500 \), \( I(0)=1 \),
and \( R(0)=0 \). Choose \( \nu = 0.1 \), \( \Delta t = 0.5 \), and \( t\in [0, 60] \).
Time \( t \) here counts days.
Visualize first how the disease develops when
\( \beta = 0.0005 \). Certain precautions, like
staying inside, will reduce \( \beta \). Try \( \beta = 0.0001 \) and comment
from the plot how a reduction in \( \beta \) influences \( S(t) \). (Put the comment
as a multi-line string in the bottom of the program file.)
Filename: SIR.
The parameters \( \nu \) and \( \beta \) in the SIR model in Exercise 41: Simulate the spreading of a disease by a SIR model can be constants or functions of time. Now we shall make an implementation of the \( f(u,t) \) function specifying the ODE system such that \( \nu \) and \( \beta \) can be given as either a constant or a Python function. Introduce a class for \( f(u,t) \), with the following code sketch:
class ProblemSIR(object):
def __init__(self, nu, beta, S0, I0, R0, T):
"""
nu, beta: parameters in the ODE system
S0, I0, R0: initial values
T: simulation for t in [0,T]
"""
if isinstance(nu, (float,int)): # number?
self.nu = lambda t: nu # wrap as function
elif callable(nu):
self.nu = nu
# same for beta and self.beta
...
# store the other parameters
def __call__(self, u, t):
"""Right-hand side function of the ODE system."""
S, I, R = u
return [-self.beta(t)*S*I, # S equation
..., # I equation
self.nu(t)*I] # R equation
# Example:
problem = ProblemSIR(beta=lambda t: 0.0005 if t <= 12 else 0.0001,
nu=0.1, S0=1500, I0=1, R0=0, T=60)
solver = ODESolver.ForwardEuler(problem)
ProblemSIR based on the sketch
of ideas above. The \( \nu \) parameter is usually not varying with time
as \( 1/\nu \) is a characteristic size of the period a person is sick,
but introduction of new medicine during the disease
might change the picture such that
time dependence becomes relevant.
We can also make a class SolverSIR for solving the problem
(see the section Example: The logistic equation with problem and solver classes for similar examples):
class SolverSIR(object):
def __init__(self, problem, dt):
self.problem, self.dt = problem, dt
def solve(self, method=ODESolver.RungeKutta4):
self.solver = method(self.problem)
ic = [self.problem.S0, self.problem.I0, self.problem.R0]
self.solver.set_initial_condition(ic)
n = int(round(self.problem.T/float(self.dt)))
t = np.linspace(0, self.problem.T, n+1)
u, self.t = self.solver.solve(t)
self.S, self.I, self.R = u[:,0], u[:,1], u[:,2]
def plot(self):
# plot S(t), I(t), and R(t)
After the breakout of a disease, authorities often start campaigns for
decreasing the spreading of the disease. Suppose a massive campaign
telling people to wash their hands more frequently is launched, with
the effect that \( \beta \) is significantly reduced after some days.
For the specific case simulated in Exercise 41: Simulate the spreading of a disease by a SIR model, let
$$
\begin{equation*} \beta(t) = \left\lbrace\begin{array}{ll}
0.0005,& 0\leq t\leq 12,\\
0.0001,& t > 12\end{array}\right.\end{equation*}
$$
Simulate this scenario with the Problem and Solver classes.
Report the maximum number of infected people and compare it to the
case where \( \beta(t) = 0.0005 \).
Filename: SIR_class.
We shall now extend the SIR model in Exercise 41: Simulate the spreading of a disease by a SIR model with a vaccination program. If a fraction \( p \) of the susceptibles per time unit is being vaccinated, and we say that the vaccination is 100 percent effective, \( pS\Delta t \) individuals will be removed from the S category in a time interval \( \Delta t \). We place the vaccinated people in a new category V. The equations for \( S \) and \( V \) becomes $$ \begin{align} S' &= -\beta SI -pS, \tag{81}\\ V' &= pS\tp \tag{82} \end{align} $$ The equations for \( I \) and \( R \) are not affected. The initial condition for \( V \) can be taken as \( V(0)=0 \). The resulting model is named SIRV.
Try the same parameters as in Exercise 41: Simulate the spreading of a disease by a SIR model in combination with \( p=0.1 \) and compute the evolution of \( S(t) \), \( I(t) \), \( R(t) \), and \( V(t) \). Comment on the effect of vaccination on the maximum number of infected.
Hint. You can of course edit the code from Exercise 42: Introduce problem and solver classes in the SIR model, but it is much better to avoid duplicating code and use object-oriented programming to implement the extensions in the present exercise as subclasses of the classes from Exercise 42: Introduce problem and solver classes in the SIR model.
Filename: SIRV.
Let the vaccination campaign in Exercise 43: Introduce vaccination in a SIR model
start 6 days after the outbreak of the disease and
let it last for 10 days,
$$
\begin{equation*} p(t) = \left\lbrace\begin{array}{ll}
0.1,& 6\leq t\leq 15,\\
0,& \hbox{otherwise} \end{array}\right.\end{equation*}
$$
Plot the corresponding solutions \( S(t) \), \( I(t) \), \( R(t) \), and \( V(t) \).
(It is clearly advantageous to have the SIRV model implemented
as an extension to the classes in Exercise 42: Introduce problem and solver classes in the SIR model.)
Filename: SIRV_varying_p.
Let the vaccination campaign in Exercise 44: Introduce a vaccination campaign in a SIR model
last for \( V_T \) days:
$$
\begin{equation*} p(t) = \left\lbrace\begin{array}{ll}
0.1,& 6\leq t\leq 6 + V_T,\\
0,& \hbox{otherwise} \end{array}\right.\end{equation*}
$$
Compute the maximum number of infected people, \( \max_t I(t) \), as
a function of \( V_T\in [0,31] \), by running the model for
\( V_T=0,1,2\ldots,31 \).
Plot this function. Determine from the plot
the optimal \( V_T \), i.e., the smallest
vaccination period
\( V_T \) such that increasing \( V_T \) has negligible effect on the maximum
number of infected people.
Filename: SIRV_optimal_duration.
Suppose the human population is attacked by zombies. This is quite a common happening in movies, and the "zombification" of humans acts much like the spreading of a disease. Let us make a differential equation model, inspired by the SIR model from Exercise 41: Simulate the spreading of a disease by a SIR model, to simulate how humans and zombies interact.
We introduce four categories of individuals:
The type of zombies considered here is inspired by the standard for modern zombies set by the classic movie The Night of the Living Dead, by George A. Romero from 1968. Only a small extension of the SIR model is necessary to model the effect of human-zombie interaction mathematically. A fraction of the human susceptibles is getting bitten by zombies and moves to the infected category. A fraction of the infected is then turned into zombies. On the other hand, humans can conquer zombies.
Now we shall precisely set up all the dynamic features of the human-zombie populations we aim to model. Changes in the S category are due to three effects:
The balance of the S category is then $$ \begin{equation*} S' = \Sigma - \beta SZ - \delta_S S,\end{equation*} $$ in the limit \( \Delta t\rightarrow 0 \).
The infected category gets a contribution \( \Delta t\beta SZ \) from the S category, but loses individuals to the Z and R category. That is, some infected are turned into zombies, while others die. Movies reveal that infected may commit suicide or that others (susceptibles) may kill them. Let \( \delta_I \) be the probability of being killed in a unit time interval. During time \( \Delta t \), a total of \( \delta_I\Delta t I \) will die and hence be transferred to the removed category. The probability that a single infected is turned into a zombie during a unit time interval is denoted by \( \rho \), so that a total of \( \Delta t\rho I \) individuals are lost from the I to the Z category in time \( \Delta t \). The accounting in the I category becomes $$ \begin{equation*} I' = \beta SZ -\rho I -\delta_I I\tp\end{equation*} $$
The zombie category gains \( -\Delta t\rho I \) individuals from the I category. We disregard the effect that any removed individual can turn into a zombie again, as we consider that effect as pure magic beyond reasonable behavior, at least according to what is observed in the Romero movie tradition. A fundamental feature in zombie movies is that humans can conquer zombies. Here we consider zombie killing in a "man-to-man" human-zombie fight. This interaction resembles the nature of zombification (or the susceptible-infective interaction in the SIR model) and can be modeled by a loss \( -\alpha SZ \) for some parameter \( \alpha \) with an interpretation similar to that of \( \beta \). The equation for \( Z \) then becomes $$ \begin{equation*} Z' = \rho I - \alpha SZ\tp\end{equation*} $$
The accounting in the R category consists of a gain \( \delta S \) of natural deaths from the S category, a gain \( \delta I \) from the I category, and a gain \( \alpha SZ \) from defeated zombies: $$ \begin{equation*} R' = \delta_S S + \delta_I I + \alpha SZ\tp\end{equation*} $$
The complete SIZR model for human-zombie interaction can be summarized as $$ \begin{align} S' &= \Sigma - \beta SZ - \delta_S S, \tag{83}\\ I' &= \beta SZ -\rho I -\delta_I I, \tag{84}\\ Z' &= \rho I - \alpha SZ , \tag{85}\\ R' &= \delta_S S + \delta_I I + \alpha SZ\tp \tag{86} \end{align} $$ The interpretations of the parameters are as follows:
Implement the SIZR model with a Problem and Solver class as
explained in Exercise 42: Introduce problem and solver classes in the SIR model, allowing parameters to vary
in time. The time variation is essential to make a realistic model
that can mimic what happens in movies.
Test the implementation with the following data: \( \beta = 0.0012 \),
\( \alpha = 0.0016 \), \( \delta_I = 0.014 \), \( \Sigma =2 \), \( \rho = 1 \),
\( S(0)=10 \), \( Z(0)=100 \), \( I(0)=0 \), \( R(0)=0 \), and simulation time \( T=24 \)
hours. Other parameters can be set to zero. These values are
estimated from the hysterical phase of the movie The Night of the
Living Dead. The time unit is hours. Plot the \( S \), \( I \), \( Z \), and \( R \)
quantities.
Filename: SIZR.
The movie The Night of the Living Dead has three phases:
Hint. It becomes necessary to work with piecewise constant functions in time. These can be hardcoded for each special case, our one can employ a ready-made tool for such functions:
from scitools.std import PiecewiseConstant
# Define f(t) as 1.5 in [0,3], 0.1 in [3,4] and 1 in [4,7]
f = PiecewiseConstant(domain=[0, 7],
data=[(0, 1.5), (3, 0.1), (4, 1)])
Filename: Night_of_the_Living_Dead.
A war on zombies can be implemented through large-scale effective attacks. A possible model is to increase \( \alpha \) in the SIZR model from Exercise 46: Simulate human-zombie interaction by some additional amount \( \omega (t) \), where \( \omega (t) \) varies in time to model strong attacks at \( m+1 \) distinct points of time \( T_0 < T_1 < \cdots < T_m \). Around these \( t \) values we want \( \omega \) to have a large value, while in between the attacks \( \omega \) is small. One possible mathematical function with this behavior is a sum of Gaussian functions: $$ \begin{equation} \omega (t) = a\sum_{i=0}^m \exp{\left(-\frac{1}{2}\left({t - T_i\over\sigma} \right)^2\right)}, \tag{87} \end{equation} $$ where \( a \) measures the strength of the attacks (the maximum value of \( \omega(t) \)) and \( \sigma \) measures the length of the attacks, which should be much less than the time between the points of attack: typically, \( 4\sigma \) measures the length of an attack, and we must have \( 4\sigma \ll T_i-T_{i-1} \) for \( i=1,\ldots,m \). We should choose \( a \) significantly larger than \( \alpha \) to make the attacks in the war on zombies much stronger than the usual "man-to-man" killing of zombies.
Modify the model and the implementation from Exercise 46: Simulate human-zombie interaction to include a war on zombies. We start out with 50
humans and 3 zombies and \( \beta=0.03 \). This leads to rapid
zombification. Assume that there are some small resistances against
zombies from the humans, \( \alpha = 0.2\beta \), throughout the
simulations. In addition, the humans implement three strong attacks,
\( a=50\alpha \), at 5, 10, and 18 hours after the zombification
starts. The attacks last for about 2 hours (\( \sigma = 0.5 \)). Set
\( \delta_S=\Delta_I=\Sigma=0 \), \( \beta =0.03 \), and \( \rho=1 \), simulate
for \( T=20 \) hours, and see if the war on zombies modeled by the
suggested \( \omega(t) \) is sufficient to save mankind.
Filename: war_on_zombies.
Suppose we have two species in an environment: a predator and a prey.
How will the two populations interact and change with time?
A system of ordinary differential equations can give insight into
this question. Let \( x(t) \) and \( y(t) \) be the size of the prey and
the predator populations, respectively. In the absence of a predator,
the population of the prey will follow the ODE
$$
\begin{equation*} \frac{dx}{dt} = rx,\end{equation*}
$$
with \( r>0 \), assuming there are enough resources for exponential growth.
Similarly, in the absence of prey, the predator population will just
experience a death rate \( m>0 \):
$$
\begin{equation*} {dy\over dt} = -my\tp\end{equation*}
$$
In the presence of the predator, the prey population will experience a
reduction in the growth proportional to \( xy \). The number of
interactions (meetings) between \( x \) and \( y \) numbers of animals is
\( xy \), and in a certain fraction of these interactions the predator
eats the prey. The predator population will correspondingly
experience a growth in the population because of the \( xy \) interactions
with the prey population. The adjusted growth of both populations can
now be expressed as
$$
\begin{align}
\frac{dx}{dt} &= rx - axy,
\tag{88}\\
{dy\over dt} &= -my + bxy,
\tag{89}
\end{align}
$$
for positive constants \( r \), \( m \), \( a \), and \( b \). Solve this system and
plot \( x(t) \) and \( y(t) \) for \( r=m=1 \), \( a=0.3 \), \( b=0.2 \), \( x(0)=1 \), and
\( y(0)=1 \), \( t\in [0,20] \). Try to explain the dynamics of the population
growth you observe. Experiment with other values of \( a \) and \( b \).
Filename: predator_prey.
In this and subsequent exercises we shall deal with the following second-order ordinary differential equation with two initial conditions: $$ \begin{equation} m\ddot u + f(\dot u) + s(u) = F(t),\quad t>0,\quad u(0)=U_0,\ \dot u(0)=V_0 \tp \tag{90} \end{equation} $$ The notation \( \dot u \) and \( \ddot u \) means \( u^{\prime}(t) \) and \( u^{\prime\prime}(t) \), respectively. Write (90) as a system of two first-order differential equations. Also set up the initial condition for this system.
Equation (90) has a wide range of applications throughout science and engineering. A primary application is damped spring systems in, e.g., cars and bicycles: \( u \) is the vertical displacement of the spring system attached to a wheel; \( \dot u \) is then the corresponding velocity; \( F(t) \) resembles a bumpy road; \( s(u) \) represents the force from the spring; and \( f(\dot u) \) models the damping force (friction) in the spring system. For this particular application \( f \) and \( s \) will normally be linear functions of their arguments: \( f(\dot u)=\beta\dot u \) and \( s(u)=ku \), where \( k \) is a spring constant and \( \beta \) some parameter describing viscous damping.
Equation (90) can also be used to describe the motions of a moored ship or oil platform in waves: the moorings act as a nonlinear spring \( s(u) \); \( F(t) \) represents environmental excitation from waves, wind, and current; \( f(\dot u) \) models damping of the motion; and \( u \) is the one-dimensional displacement of the ship or platform.
Oscillations of a pendulum can be described by (90): \( u \) is the angle the pendulum makes with the vertical; \( s(u)=(mg/L)\sin (u) \), where \( L \) is the length of the pendulum, \( m \) is the mass, and \( g \) is the acceleration of gravity; \( f(\dot u) = \beta |\dot u|\dot u \) models air resistance (with \( \beta \) being some suitable constant, and \( F(t) \) might be some motion of the top point of the pendulum.
Another application is electric circuits with \( u(t) \) as the charge, \( m=L \) as the inductance, \( f(\dot u)=R\dot u \) as the voltage drop across a resistor \( R \), \( s(u)=u/C \) as the voltage drop across a capacitor \( C \), and \( F(t) \) as an electromotive force (supplied by a battery or generator).
Furthermore, Equation (90) can act as a (very) simplified model of many other oscillating systems: aircraft wings, lasers, loudspeakers, microphones, tuning forks, guitar strings, ultrasound imaging, voice, tides, the El Ni\ {n}o phenomenon, climate changes – to mention some.
Make a function
def rhs(u, t):
...
u argument is an array or list with the two solution components
u[0] and u[1] at some time t.
Inside rhs, assume that you have access to three global Python functions
friction(dudt), spring(u), and external(t) for
evaluating \( f(\dot u) \), \( s(u) \), and \( F(t) \), respectively.
Test the rhs function in combination with the functions
\( f(\dot u)=0 \), \( F(t)=0 \), \( s(u)=u \), and the choice \( m=1 \).
The differential equation then reads \( \ddot u + u = 0 \). With
initial conditions \( u(0)=1 \) and \( \dot u(0)=0 \), one can show that
the solution is given by \( u(t)=\cos(t) \). Apply three
numerical methods: the 4th-order
Runge-Kutta method and the Forward Euler method from the ODESolver
module developed in
the section The ODESolver class hierarchy, as well as the
2nd-order Runge-Kutta method developed in Exercise 31: Implement a 2nd-order Runge-Kutta method; class.
Use a time step \( \Delta t=\pi/20 \).
Plot \( u(t) \) and \( \dot u(t) \) versus \( t \) together with the exact
solutions. Also make a plot of \( \dot u \) versus \( u \)
(plot(u[:,0], u[:,1]) if u is the array returned from the solver's
solve method). In the latter case, the exact plot should be a
circle because the points on the curve are \( (\cos t, \sin t) \), which
all lie on a circle as \( t \) is varied. Observe that the ForwardEuler
method results in a spiral and investigate how the spiral develops as
\( \Delta t \) is reduced.
The kinetic energy \( K \) of the motion is given by \( \frac{1}{2}m\dot u^2 \),
and the potential energy \( P \) (stored in the spring) is given by
the work done by the spring force: \( P = \int_0^u s(v)dv = \frac{1}{2}u^2 \).
Make a plot with \( K \) and \( P \) as functions
of time for both the 4th-order Runge-Kutta method and the Forward Euler
method, for the same physical problem described above.
In this test case, the sum of the kinetic and potential
energy should be constant. Compute this constant analytically and plot
it together with the sum \( K+P \) as calculated
by the 4th-order Runge-Kutta method and the Forward Euler method.
Filename: oscillator_v1.
The solution \( u(t) \) of the equation (90) often exhibits an oscillatory behavior (for the test problem in Exercise 51: Solve \( \ddot u + u =0 \) we have that \( u(t)=\cos t \)). It is then of interest to find the wavelength of the oscillations. The purpose of this exercise is to find and visualize the distance between peaks in a numerical representation of a continuous function.
Given an array \( (y_0,\ldots,y_{n-1}) \) representing a function \( y(t) \)
sampled at various points \( t_0,\ldots,t_{n-1} \),
a local maximum of \( y(t) \) occurs at \( t=t_k \) if \( y_{k-1} < y_k>y_{k+1} \).
Similarly, a local minimum of \( y(t) \) occurs at \( t=t_k \) if
\( y_{k-1}> y_k < y_{k+1} \). By iterating over the
\( y_1,\ldots,y_{n-2} \) values and making the two tests,
one can collect local maxima and minima
as \( (t_k,y_k) \) pairs. Make a function minmax(t, y) which returns
two lists, minima and maxima, where each list holds pairs
(2-tuples) of \( t \) and \( y \) values of local
minima or maxima. Ensure that the \( t \) value increases from one
pair to the next. The arguments t and y in minmax
hold the coordinates \( t_0,\ldots,t_{n-1} \) and
\( y_0,\ldots,y_{n-1} \), respectively.
Make another function wavelength(peaks) which
takes a list peaks of 2-tuples with \( t \) and \( y \) values for
local minima or maxima as argument and returns an array of distances
between consecutive \( t \) values, i.e., the distances between the peaks.
These distances reflect the local wavelength of the computed \( y \) function.
More precisely, the first element in the returned array is
peaks[1][0]-peaks[0][0], the next element is
peaks[2][0]-peaks[1][0], and so forth.
Test the minmax and wavelength functions on \( y \) values
generated by \( y=e^{t/4}\cos (2t) \) and
\( y=e^{-t/4}\cos (t^2/5) \) for \( t\in [0,4\pi] \).
Plot the \( y(t) \) curve in each case, and mark the local minima and maxima
computed by minmax with circles and
boxes, respectively. Make a separate plot with the array
returned from the wavelength function (just plot the array
against its indices - the point is to see if the wavelength
varies or not). Plot only the wavelengths corresponding to maxima.
Make a module with the minmax and wavelength function, and
let the test block perform the tests specified above.
Filename: wavelength.
The user-chosen functions \( f \), \( s \), and \( F \) in
Exercise 51: Solve \( \ddot u + u =0 \) must be coded with particular names.
It is then difficult to have several functions for \( s(u) \) and
experiment with these. A much more flexible code arises if we
adopt the ideas of a problem and a solver class as
explained in the section Example: The logistic equation with problem and solver classes.
Specifically, we shall here make use of class Problem3
in the section Example: The logistic equation with problem and solver classes to store information
about \( f(\dot u) \), \( s(u) \), \( F(t) \), \( u(0) \), \( \dot u(0) \),
\( m \), \( T \), and the exact solution (if available).
The solver class can store parameters related to the numerical
quality of the solution, i.e., \( \Delta t \) and the name of
the solver class in the ODESolver hierarchy.
In addition we will make a visualizer class for producing plots
of various kinds.
We want all parameters to be set on the command line, but also have
sensible default values. As in the section Example: The logistic equation with problem and solver classes,
the argparse module is used to read data from the command line.
Class Problem can be sketched as follows:
class Problem(object):
def define_command_line_arguments(self, parser):
"""Add arguments to parser (argparse.ArgumentParser)."""
parser.add_argument(
'--friction', type=func_dudt, default='0',
help='friction function f(dudt)',
metavar='<function expression>')
parser.add_argument(
'--spring', type=func_u, default='u',
help='spring function s(u)',
metavar='<function expression>')
parser.add_argument(
'--external', type=func_t, default='0',
help='external force function F(t)',
metavar='<function expression>')
parser.add_argument(
'--u_exact', type=func_t_vec, default='0',
help='exact solution u(t) (0 or None: now known)',
metavar='<function expression>')
parser.add_argument(
'--m', type=evalcmlarg, default=1.0, help='mass',
type=float, metavar='mass')
...
return parser
def set(self, args):
"""Initialize parameters from the command line."""
self.friction = args.friction
self.spring = args.spring
self.m = args.m
...
def __call__(self, u, t):
"""Define the right-hand side in the ODE system."""
m, f, s, F = \
self.m, self.friction, self.spring, self.external
...
type argument
to parser.add_argument for turning strings into proper
objects, in particular StringFunction objects with different
independent variables:
def evalcmlarg(text):
return eval(text)
def func_dudt(text):
return StringFunction(text, independent_variable='dudt')
def func_u(text):
return StringFunction(text, independent_variable='u')
def func_t(text):
return StringFunction(text, independent_variable='t')
def func_t_vec(text):
if text == 'None' or text == '0':
return None
else:
f = StringFunction(text, independent_variable='t')
f.vectorize(globals())
return f
evalcmlarg is essential: this function runs
the strings from the command line through eval, which means
that we can use mathematical formulas like --T '4*pi'.
Class Solver is relatively much shorter
than class Problem:
class Solver(object):
def __init__(self, problem):
self.problem = problem
def define_command_line_arguments(self, parser):
"""Add arguments to parser (argparse.ArgumentParser)."""
# add --dt and --method
...
return parser
def set(self, args):
self.dt = args.dt
self.n = int(round(self.problem.T/self.dt))
self.solver = eval(args.method)
def solve(self):
self.solver = self.method(self.problem)
ic = [self.problem.initial_u, self.problem.initial_dudt]
self.solver.set_initial_condition(ic)
time_points = linspace(0, self.problem.T, self.n+1)
self.u, self.t = self.solver.solve(time_points)
The Visualizer class holds references to a
Problem and Solver instance and creates plots.
The user can specify plots in an interactive dialog in the terminal
window. Inside a loop, the user is repeatedly asked to specify a plot
until the user responds with quit. The specification of a plot
can be one of the words u, dudt, dudt-u,
K, and wavelength which means
a plot of \( u(t) \) versus \( t \), \( \dot u(t) \) versus \( t \),
\( \dot u \) versus \( u \), \( K \) (\( =\frac{1}{2}m\dot u^2 \), kinetic energy)
versus \( t \), and \( u \)'s wavelength versus
its indices, respectively. The wavelength can be computed
from the local maxima of \( u \)
as explained in Exercise 52: Make a tool for analyzing oscillatory solutions.
A sketch of class Visualizer is given
next:
class Visualizer(object):
def __init__(self, problem, solver):
self.problem = problem
self.solver = solver
def visualize(self):
t = self.solver.t # short form
u, dudt = self.solver.u[:,0], self.solver.u[:,1]
# Tag all plots with numerical and physical input values
title = 'solver=%s, dt=%g, m=%g' % \
(self.solver.method, self.solver.dt, self.problem.m)
# Can easily get the formula for friction, spring and force
# if these are string formulas.
if isinstance(self.problem.friction, StringFunction):
title += ' f=%s' % str(self.problem.friction)
if isinstance(self.problem.spring, StringFunction):
title += ' s=%s' % str(self.problem.spring)
if isinstance(self.problem.external, StringFunction):
title += ' F=%s' % str(self.problem.external)
# Let the user interactively specify what
# to be plotted
plot_type = ''
while plot_type != 'quit':
plot_type = raw_input('Specify a plot: ')
figure()
if plot_type == 'u':
# Plot u vs t
if self.problem.u_exact is not None:
hold('on')
# Plot self.problem.u_exact vs t
show()
savefig('tmp_u.pdf')
elif plot_type == 'dudt':
...
elif plot_type == 'dudt-u':
...
elif plot_type == 'K':
...
elif plot_type == 'wavelength':
...
Make a complete implementation of the three proposed classes. Also
make a main function that (i) creates a problem, solver, and
visualizer, (ii) calls the functions to define command-line arguments
in the problem and solver classes, (iii) reads the command line, (iv)
passes on the command-line parser object to the problem and solver
classes, (v) calls the solver, and (vi) calls the visualizer's
visualize method to create plots. Collect the classes and
functions in a module oscillator, which has a call to main
in the test block.
The first task from Exercise 51: Solve \( \ddot u + u =0 \) can now be run as
oscillator.py --method ForwardEuler --u_exact "cos(t)" \
--dt "pi/20" --T "5*pi"
Explore some of the possibilities of specifying several functions on the command line:
oscillator.py --method RungeKutta4 --friction "0.1*dudt" \
--external "sin(0.5*t)" --dt "pi/80" \
--T "40*pi" --m 10
oscillator.py --method RungeKutta4 --friction "0.8*dudt" \
--external "sin(0.5*t)" --dt "pi/80" \
--T "120*pi" --m 50
oscillator.
Some typical choices of \( f(\dot u) \), \( s(u) \), and \( F(t) \) in (90) are listed below:
functions where each of the choices above are
implemented as a class with a __call__ special method. Also add a
class Zero for a function whose value is always zero. It is natural
that the parameters in a function are set as arguments to the
constructor. The different classes for spring functions can all have
a common base class holding the \( k \) parameter as data attribute.
Filename: functions.
The purpose of this exercise is to demonstrate the use of the classes from Exercise 54: Use classes for flexible choices of models to solve problems described by (90).
With a lot of models for \( f(\dot u) \), \( s(u) \), and \( F(t) \) available
as classes in functions.py, the initialization of self.friction,
self.spring, etc., from the command line does not work, because
we assume simple string formulas on the command line. Now we
want to write things like --spring 'LinearSpring(1.0)'.
There is a quite simple remedy: replace all the special conversion
functions to StringFunction objects by evalcmlarg in
the type specifications in the parser.add_argument calls.
If a from functions import * is also performed in the
oscillator.py file, a simple eval will turn
strings like 'LinearSpring(1.0)' into living objects.
However, we shall here follow a simpler approach, namely dropping initializing parameters on the command line and instead set them directly in the code. Here is an example:
problem = Problem()
problem.m = 1.0
k = 1.2
problem.spring = CubicSpring(k)
problem.friction = Zero()
problem.T = 8*pi/sqrt(k)
...
functions module.
Note that the set method in classes Solver and
Visualizer is unaffected
by the new objects from the functions module, so flexible initialization
via command-line arguments works as before for --dt,
--method, and plot. One may also dare to call
the set method in the problem object to set parameters
like m, initial_u, etc., or one can choose the safer
approach of not calling set but initialize all data attributes
explicitly in the user's code.
Make a new file say oscillator_test.py
where you import class Problem, Solver, and Visualizer,
plus all classes from the functions module.
Provide a main1 function for solving the following problem:
\( m=1 \), \( u(0)=1 \), \( \dot u(0)=0 \),
no friction (use class Zero),
no external forcing (class Zero),
a linear spring \( s(u)=u \), \( \Delta t=\pi/20 \), \( T=8\pi \), and
exact \( u(t)=\cos (t) \). Use the Forward Euler method.
Then make another function main2 for the
case with
\( m=5 \), \( u(0)=1 \), \( \dot u(0)=0 \),
linear friction \( f(\dot u)=0.1\dot u \), \( s(u)=u \), \( F(t)=\sin (\frac{1}{2}t) \),
\( \Delta t =\pi/80 \), \( T=60\pi \), and no knowledge of an exact solution.
Use the 4-th order Runge-Kutta method.
Let a test block use the first command-line argument to indicate
a call to main1 or main2.
Filename: oscillator_test.
A population of fish is governed by the differential equation $$ \begin{equation} \begin{array}{c} \frac{dx}{dt}=\frac{1}{10}x\left( 1-\frac{x}{100}\right) -h, \quad x(0) =500, \tag{91} \end{array} \end{equation} $$ where \( x(t) \) is the size of the population at time \( t \) and \( h \) is the harvest.
a) Assume \( h=0 \). Find an exact solution for \( x(t) \). For which value of \( t \) is \( \frac{dx}{dt} \) largest? For which value of \( t \) is \( \frac{1}{x}\frac{dx}{dt} \) largest?
b) Solve the differential equation (91) by the Forward Euler method. Plot the numerical and exact solution in the same plot.
c) Suppose the harvest \( h \) depends on the fishers' efforts, \( E \), in the following way: \( h=qxE \), with \( q \) as a constant. Set \( q=0.1 \) and assume \( E \) is constant. Show the effect of \( E \) on \( x(t) \) by plotting several curves, corresponding to different \( E \) values, in the same figure.
d) The fishers' total revenue is given by \( \pi=ph-\frac{c}{2}E^{2} \), where \( p \) is a constant. In the literature about the economy of fisheries, one is often interested in how a fishery will develop in the case the harvest is not regulated. Then new fishers will appear as long as there is money to earn (\( \pi > 0 \)). It can (for simplicity) be reasonable to model the dependence of \( E \) on \( \pi \) as $$ \begin{equation} {dE\over dt} = \gamma\pi, \tag{92} \end{equation} $$ where \( \gamma \) is a constant. Solve the system of differential equations for \( x(t) \) and \( E(t) \) by the 4th-order Runge-Kutta method, and plot the curve with points (\( x(t),E(t) \)) in the two cases \( \gamma =1/2 \) and \( \gamma\rightarrow\infty \). Choose \( c=0.3 \), \( p=10 \), \( E(0)=0.5 \), and \( T=1 \).
Filename: fishery.