$$ \newcommand{\tp}{\thinspace .} $$

This chapter is taken from the book A Primer on Scientific Programming with Python by H. P. Langtangen, 5th edition, Springer, 2016.

Exercises

Exercise 1: Solve a simple ODE with function-based code

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.

Exercise 2: Solve a simple ODE with class-based code

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.

Exercise 3: Solve a simple ODE with the ODEsolver hierarchy

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.

Exercise 4: Solve an ODE specified on the command line

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.

Exercise 5: Implement a numerical method for ODEs

Implement the numerical method (36)-(37). How can you verify that the implementation is correct? Filename: Heuns_method.

Exercise 6: Solve an ODE for emptying a tank

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.

Exercise 7: Solve an ODE for the arc length

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.

Exercise 8: Simulate a falling or rising body in a fluid

A body moving vertically through a fluid (liquid or gas) is subject to three different types of forces:

(Roughly speaking, the \( F_d \) formula is suitable for medium to high velocities, while for very small velocities, or very small bodies, \( F_d \) is proportional to the velocity, not the velocity squared, see [2].)

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.

Exercise 9: Verify the limit of a solution as time grows

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.

Exercise 10: Scale the logistic equation

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.

Exercise 11: Compute logistic growth with time-varying carrying capacity

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.

Exercise 12: Solve an ODE until constant solution

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.

Exercise 13: Use a problem class to hold data about an ODE

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)
When you become familiar with such code, and appreciate it, you can call yourself a professional programmer - with a deep knowledge of how lists, dictionaries, and classes can play elegantly together to conduct scientific experiments. In the present case we perform only a few experiments that could also have been done by six separate calls to the solver functionality, but in large-scale scientific and engineering investigations with hundreds of parameter combinations, the above code is still the same, only the generation of the Problem instances becomes more involved.

Filename: pizza_cooling2.

Exercise 14: Derive and solve a scaled ODE problem

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.

Exercise 15: Clean up a file to make it a module

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)
The reason is that this 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.

Exercise 16: Simulate radioactive decay

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.

Exercise 17: Compute inverse functions by solving an ODE

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.

Exercise 18: Make a class for computing inverse functions

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'])
Here, 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.

Exercise 19: Add functionality to a class

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.

Exercise 20: Compute inverse functions by interpolation

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)
The 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.

Exercise 21: Code the 4th-order Runge-Kutta method; function

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.

Exercise 22: Code the 4th-order Runge-Kutta method; class

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.

Exercise 23: Compare ODE methods

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.

Exercise 24: Code a test function for systems of ODEs

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.

Exercise 25: Code Heun's method for ODE systems; function

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.

Exercise 26: Code Heun's method for ODE systems; class

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.

Exercise 27: Implement and test the Leapfrog method

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.

Exercise 28: Implement and test an Adams-Bashforth method

Do Exercise 27: Implement and test the Leapfrog method with the 3rd-order Adams-Bashforth method (46). Filename: AdamBashforth3.

Exercise 29: Solve two coupled ODEs for radioactive decay

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.

Exercise 30: Implement a 2nd-order Runge-Kutta method; function

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.

Exercise 31: Implement a 2nd-order Runge-Kutta method; class

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.

Exercise 32: Code the iterated midpoint method; function

a) Implement the numerical method (47)-(48) as a function

iterated_Midpoint_method(f, U0, T, n, N)
where 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.

Exercise 33: Code the iterated midpoint method; class

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.

Exercise 34: Make a subclass for the iterated midpoint method

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.

Exercise 35: Compare the accuracy of various methods for ODEs

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.

Exercise 36: Animate how various methods for ODEs converge

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.

Exercise 37: Study convergence of numerical methods for ODEs

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=[])
that returns a series of estimated \( r \) values corresponding to a series of \( \Delta t \) values given as the 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.

Exercise 38: Find a body's position along with its velocity

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.

Exercise 39: Add the effect of air resistance on a ball

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.

Exercise 40: Solve an ODE system for an electric circuit

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.

Remarks

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.

Exercise 41: Simulate the spreading of a disease by a SIR 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.

Exercise 42: Introduce problem and solver classes in the SIR model

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)
Write the complete code for class 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.

Exercise 43: Introduce vaccination in a SIR model

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.

Exercise 44: Introduce a vaccination campaign in a SIR model

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.

Exercise 45: Find an optimal vaccination period

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.

Exercise 46: Simulate human-zombie interaction

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:

  1. S: susceptible humans who can become zombies.
  2. I: infected humans, being bitten by zombies.
  3. Z: zombies.
  4. R: removed individuals, either conquered zombies or dead humans.
The corresponding functions counting how many individuals we have in each category are named \( S(t) \), \( I(t) \), \( Z(t) \), and \( R(t) \), respectively.

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:

  1. Susceptibles are infected by zombies, modeled by a term \( -\Delta t\beta SZ \), similar to the S-I interaction in the SIR model.
  2. Susceptibles die naturally or get killed and therefore enter the removed category. If the probability that one susceptible dies during a unit time interval is \( \delta_S \), the total expected number of deaths in a time interval \( \Delta t \) becomes \( \Delta t\delta_S S \).
  3. We also allow new humans to enter the area with zombies, as this effect may be necessary to successfully run a war on zombies. The number of new individuals in the S category arriving per time unit is denoted by \( \Sigma \), giving an increase in \( S(t) \) by \( \Delta t\Sigma \) during a time \( \Delta t \).
We could also add newborns to the S category, but we simply skip this effect since it will not be significant over time scales of a few days.

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:

Note that probabilities per unit time do not necessarily lie in the interval \( [0,1] \). The real probability, lying between 0 and 1, arises after multiplication by the time interval of interest.

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.

Exercise 47: Simulate a zombie movie

The movie The Night of the Living Dead has three phases:

  1. The initial phase, lasting for (say) 4 hours, where two humans meet one zombie and one of the humans get infected. A rough (and uncertain) estimation of parameters in this phase, taking into account dynamics not shown in the movie, yet necessary to establish a more realistic evolution of the S and Z categories later in the movie, is \( \Sigma =20 \), \( \beta = 0.03 \), \( \rho = 1 \), \( S(0)=60 \), and \( Z(0)=1 \). All other parameters are taken as zero when not specified.
  2. The hysterical phase, when the zombie treat is evident. This phase lasts for 24 hours, and relevant parameters can be taken as \( \beta = 0.0012 \), \( \alpha = 0.0016 \), \( \delta_I = 0.014 \), \( \Sigma =2 \), \( \rho = 1 \).
  3. The counter attack by humans, estimated to last for 5 hours, with parameters \( \alpha = 0.006 \), \( \beta = 0 \) (humans no longer get infected), \( \delta_S=0.0067 \), \( \rho = 1 \).
Use the program from Exercise 46: Simulate human-zombie interaction to simulate all three phases of the movie.

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.

Exercise 48: Simulate a war on zombies

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.

Exercise 49: Explore predator-prey population interactions

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.

Exercise 50: Formulate a 2nd-order ODE as a system

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.

Physical applications

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.

Exercise 51: Solve \( \ddot u + u =0 \)

Make a function

def rhs(u, t):
    ...
for returning a list with two elements with the two right-hand side expressions in the first-order differential equation system from Exercise 50: Formulate a 2nd-order ODE as a system. As usual, the 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.

Exercise 52: Make a tool for analyzing oscillatory solutions

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.

Exercise 53: Implement problem, solver, and visualizer classes

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
        ...
Several functions are specified as the 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
The use of 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"
The other tasks from Exercise 51: Solve \( \ddot u + u =0 \) can be tested similarly.

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
Filename: oscillator.

Exercise 54: Use classes for flexible choices of models

Some typical choices of \( f(\dot u) \), \( s(u) \), and \( F(t) \) in (90) are listed below:

Make a module 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.

Exercise 55: Apply software for oscillating systems

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)
...
This is the simplest way of making use of the objects in the 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.

Exercise 56: Model the economy of fishing

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.

References

  1. L. S. Lerner. Physics for Scientists and Engineers, Jones and Barlett, 1996.
  2. F. M. White. Fluid Mechanics, 2nd edition, McGraw-Hill, 1986.
  3. H. P. Langtangen. Object-oriented programming, \emphhttp://hplgit.github.io/primer.html/doc/pub/oo, http://hplgit.github.io/primer.html/doc/pub/oo.