## Transcribed Text

1. Introduction
This tutorial is designed to familiarise you with solving ordinary differential equations using the Python function odeint. Note, in the lectures the scipy class ode is used to get a more elegant interface. Here we cover the function odeint from scipy.integrate.
2. Scalar Ordinary Differential Equation (ODEs) Suppose we want to find the solution u(t) of the ordinary differential equation
(2.1) du = ku(t)(1 − u(t)), t ≥ 0. dt
This equation is a simple model for the spread of a disease. Consider a herd of P animals, and let u(t) be the proportion of animals infected by the disease after t days. Thus 1 − u(t) is the proportion of uninfected animals. New infections occur when infected and uninfected animals meet. The number of such meetings is proportional to u(t)(1 − u(t)). Thus a simple model for the spread of the disease is that the new infections per day, du/dt, is proportional to u(t)(1 − u(t)). That is, we obtain Equation (2.1) where the constant k depends on the density of the animals and the infectiousness of the disease. As a specific example, suppose here P = 10000 animals and k = 0.2. Initially at time t = 0 there is just one infectious animal, and so we add the initial condition
(2.2) u(0) = 1/P = 1/10000.
We want to determine the proportion of infected individuals over the next 100 days, that is to plot u(t) from t = 0 to t = 100. To solve this problem we need to define a function to evaluate the right hand side of (2.1), and then use the Python function odeint. Create a Python function which specifies the right hand side of the differential equation. Make sure you get the calling sequence correct. For instance
def f(y,t,k):
return k∗y∗(1−y)
(Although t is not used in calculating the result, the Python function odeint requires the function f must have at least two input arguments (in the correct order). Use help odeint for more information.)
Having set up the problem, the differential equation is solved as follows:
from scipy . integrate import odeint from scipy import arange
from pylab import plot , show
# define the domain
xs = arange (0 , 101)
# solve the ode
1
2
ys = odeint(f, 1./10000, xs, args = (.2,))
# note: for args = (.2,), the com is required
# plot the solution
plot(xs, ys[:,0]) show()
The function odeint also has parameters which can be adjusted to control the accuracy of the computation. It is always a good idea to check the original solution by resolving the problem with higher accuracy requirements. For example, we can use odeint so that the estimated absolute and relative errors are kept less than 1 × 10−9. This is done with the command
# solve the ode to higher degree of accuracy
ys2 = odeint(f,1./10000, xs, args=(.2,), rtol =1e−9,atol=1e−9) # plot the two solutions
plot(xs, ys[:,0]) plot(xs,ys2[: ,0]) show()
We immediately see that the solution computed in the first call to odeint lies on top of the more accurate solution curve. That is the two solutions coincide to graphical accuracy.
For this simple model it is possible to find an analytic formula for u(t). In fact (2.3) u(t) = 1 , where c = P − 1.
1 + c exp(kt)
We can add this formula to the above graph as a final check.
from scipy import arange , exp c = 10000−1
k=.2
# define the exact solution
def exactu(t ):
return 1.0/(1.0+c∗exp(−k∗t ));
# plot the approximate solution
plot(xs, ys[:,0]);
# plot the exact solution
plot(xs, exactu(xs), ’r+’) show( )
The exact formula gives a curve that is the same as the numerical solution. In fact the difference between the exact solution and ys is at most 5.54 × 10−5 and the difference between the exact solution and ys2 is at most 3.52 × 10−6. So we can be confident in using odeint to solve the problems below.
Exercise 1
3
(a) Check that (2.3) is a solution to (2.1). That is, when we differentiate
this formula for u(t) we get the right side of (2.1). And this formula satisfies (2.2). Do this by hand.
(b) From the graph of u(t), estimate how many days are needed for half the population to be infected?
Hint: You can inspect the plot to get a solution up to the closest day. The com- mand pylab.grid will make this easier to answer. To narrow the range of the plots you may want to use the pylab.xlim and pylab.ylim commands.
(c) Suppose now that after the infection is noted at day t = 0, an in-oculation program is started from day 5, and 200 of the uninfected animals are inoculated per day. After t days, the proportion of animals inoculated per day is 200(t − 5)+/10000. Here (t − 5)+ means max(t − 5, 0). These inoculated animals cannot be infected, an new infections occur when an uninfected, un-inoculated animal meets an infected animal. Thus the rate of increase of u is now modelled
by the differential equation
du = ku(t)(1 − u(t) − .02(t − 5)+)+;
dt
Now, how many animals are uninfected after 100 days with the inoculation pro- gram?
3. System and higher order equations Suppose we want to solve the nonlinear 2nd-order initial value problem
y′′ +2yy′ +y = sin(t); y(0) = 0, y′(0) = 2.
We can do this with off-the-shelf methods in Python. We first have to write it as a system
of first-order ODEs. Setting w(t) = y′(t), we have
y′ =w, w′ =sin(t)−2yw−y
together with initial conditions y(0) = 0, w(0) = 2. To solve this numerically we have to define a function accepting two arguments: (1) a vector V whose first and second components V [0] and V [1] are the values of y and w, and (2) the value of the independent variable t. We call this function F; it should return a 2-vector containing the derivatives y′, w′. We type:
def F(V,t):
returnarray([V[1], sin(t)−2∗V[0]∗V[1]−V[0] ])
We then define the initial condition and the t-values where we wish to know the solution:
from scipy import array , linspace , pi , sin IC = array([0.,2.])
t = linspace(0,3∗pi/2,200)
Finally, we type
from scipy . integrate import odeint Y = odeint(F, IC, t)
4
This returns a 200 × 2 vector Y containing the values of y in the first column and those of w in the second column.
Exercise 2.
(a) The Lotka-Volterra equations, also known as the predator-prey equa-
tions, are a pair of first-order, non-linear, differential equations frequently used to describe the dynamics of biological systems in which two species interact, one as a predator and the other as prey. The populations change through time according to the pair of equations
dx = αx − βxy dt
dy = δxy − γy dt
where x is the number of prey (for example, rabbits), y is the number of some predator (for example, foxes); dy and dx represent the growth rates of the two
dt dt
populations over time, t represents time, and α, β, γ, δ are positive real parameters describing the interaction of the two species. Assume α = 0.7, β = δ = 0.007 and γ = 1. Modify the above example program to solve the Lokta-Volterra equations with the initial conditions x1(0) = 70 and x2(0) = 50 for 30 units of time with step size 0.1.
(b) Consider the initial value problem of the third order ordinary differ- ential equation
x′′′ +4x′′ +5x′ +2x = 2t2 +10t+8,
x(0) = 1, x′(0) = −1, x′′(0) = 3.
Check, by hand, that x(t) = e−t + t2 is the solution of the problem. Write the differential equation as a system of first-order differential equations. Modify the above example program to solve this initial value problem on the interval [0,10] with step size 0.1. Plot the computational results and the exact solution for comparison.

These solutions may offer step-by-step problem-solving explanations or good writing examples that include modern styles of formatting and construction
of bibliographies out of text citations and references. Students may use these solutions for personal skill-building and practice.
Unethical use is strictly forbidden.