# Ordinary Differential Equations

Learning Objectives: Understand the numerical solution of ODEs and use scipy.integrate.odeint to solve and explore ODEs numerically.

## Imports



In :

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import seaborn as sns



## Overview of ODEs

Many of the equations of Physics, Chemistry, Statistics, Data Science, etc. are Ordinary Differential Equation or ODEs. An ODE is a differential equation with the form:

$$\frac{d\vec{y}}{dt} = \vec{f}\left(\vec{y}(t), t\right)$$

The goal is usually to solve for the $N$ dimensional state vector $\vec{y}(t)$ at each time $t$ given some initial condition:

$$\vec{y}(0) = \vec{y}_0$$

In this case we are using $t$ as the independent variable, which is common when studying differential equations that depend on time. But any independent variable may be used, such as $x$. Solving an ODE numerically usually involves picking a set of $M$ discrete times at which we wish to know the solution:



In :

tmax = 10.0 # The max time
M = 100     # Use 100 times between [0,tmax]
t = np.linspace(0,tmax,M)
t




Out:

array([  0.        ,   0.1010101 ,   0.2020202 ,   0.3030303 ,
0.4040404 ,   0.50505051,   0.60606061,   0.70707071,
0.80808081,   0.90909091,   1.01010101,   1.11111111,
1.21212121,   1.31313131,   1.41414141,   1.51515152,
1.61616162,   1.71717172,   1.81818182,   1.91919192,
2.02020202,   2.12121212,   2.22222222,   2.32323232,
2.42424242,   2.52525253,   2.62626263,   2.72727273,
2.82828283,   2.92929293,   3.03030303,   3.13131313,
3.23232323,   3.33333333,   3.43434343,   3.53535354,
3.63636364,   3.73737374,   3.83838384,   3.93939394,
4.04040404,   4.14141414,   4.24242424,   4.34343434,
4.44444444,   4.54545455,   4.64646465,   4.74747475,
4.84848485,   4.94949495,   5.05050505,   5.15151515,
5.25252525,   5.35353535,   5.45454545,   5.55555556,
5.65656566,   5.75757576,   5.85858586,   5.95959596,
6.06060606,   6.16161616,   6.26262626,   6.36363636,
6.46464646,   6.56565657,   6.66666667,   6.76767677,
6.86868687,   6.96969697,   7.07070707,   7.17171717,
7.27272727,   7.37373737,   7.47474747,   7.57575758,
7.67676768,   7.77777778,   7.87878788,   7.97979798,
8.08080808,   8.18181818,   8.28282828,   8.38383838,
8.48484848,   8.58585859,   8.68686869,   8.78787879,
8.88888889,   8.98989899,   9.09090909,   9.19191919,
9.29292929,   9.39393939,   9.49494949,   9.5959596 ,
9.6969697 ,   9.7979798 ,   9.8989899 ,  10.        ])



It is useful to define the step size $h$ as:

$$h = t_{i+1} - t_i$$


In :

h = t-t
print("h =", h)




h = 0.10101010101



The numerical solution of an ODE will then be an $M\times N$ array $y_{ij}$ such that:

$$\left[\vec{y}(t_i)\right]_j = y_{ij}$$

In other words, the rows of the array $y_{ij}$ are the state vectors $\vec{y}(t_i)$ at times $t_i$. Here is an array of zeros having the right shape for the values of $N$ and $M$ we are using here:



In :

N = 2 # 2d case
y = np.zeros((M, N))
print("N =", N)
print("M =", M)
print("y.shape =", y.shape)




N = 2
M = 100
y.shape = (100, 2)



A numerical ODE solver takes the ith row of this array y[i,:] and calculates the i+1th row y[i+1,:]. This process starts with the initial condition y[0,:] and continues through all of the times with steps of size $h$. One of the core ideas of numerical ODE solvers is that the error at each step is proportional to $\mathcal{O}(h^n)$ where $n\geq1$. Because $h<1$ you can reduce the error by making $h$ smaller (up to a point) or finding an ODE solver with a larger value of $n$.

Here are some common numerical algorithms for solving ODEs:

1. The Euler method, which has an error of $\mathcal{O}(h)$.
2. The midpoint method, which has an error of $\mathcal{O}(h^2)$.
3. Runga-Kutta methods, the most common (called RK4) of which has an error of $\mathcal{O}(h^4)$. Because Runga-Kutta methods are fast and have a small errors, they are one of the most popular general purpose algorithm for solving ODEs.

There are many other specialized methods and tricks for solving ODEs (see this page). One of the most common tricks is to use an adaptive step size, which changes the value of $h$ at each step to make sure the error stays below a certain threshold.

## Using scipy.integrate.odeint

SciPy provides a general purpose ODE solver, scipy.integrate.odeint, that can handle a wide variety of linear and non-linear multidimensional ODEs.



In :

from scipy.integrate import odeint




In :

odeint?




Signature: odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0)
Docstring:
Integrate a system of ordinary differential equations.

Solve a system of ordinary differential equations using lsoda from the
FORTRAN library odepack.

Solves the initial value problem for stiff or non-stiff systems
of first order ode-s::

dy/dt = func(y,t0,...)

where y can be a vector.

Parameters
----------
func : callable(y, t0, ...)
Computes the derivative of y at t0.
y0 : array
Initial condition on y (can be a vector).
t : array
A sequence of time points for which to solve for y.  The initial
value point should be the first element of this sequence.
args : tuple, optional
Extra arguments to pass to function.
Dfun : callable(y, t0, ...)
Gradient (Jacobian) of func.
col_deriv : bool, optional
True if Dfun defines derivatives down columns (faster),
otherwise Dfun should define derivatives across rows.
full_output : bool, optional
True if to return a dictionary of optional outputs as the second output
printmessg : bool, optional
Whether to print the convergence message

Returns
-------
y : array, shape (len(t), len(y0))
Array containing the value of y for each desired time in t,
with the initial value y0 in the first row.
infodict : dict, only returned if full_output == True

=======  ============================================================
key      meaning
=======  ============================================================
'hu'     vector of step sizes successfully used for each time step.
'tcur'   vector with the value of t reached for each time step.
(will always be at least as large as the input times).
'tolsf'  vector of tolerance scale factors, greater than 1.0,
computed when a request for too much accuracy was detected.
'tsw'    value of t at the time of the last method switch
(given for each time step)
'nst'    cumulative number of time steps
'nfe'    cumulative number of function evaluations for each time step
'nje'    cumulative number of jacobian evaluations for each time step
'nqu'    a vector of method orders for each successful step.
'imxer'  index of the component of largest magnitude in the
weighted local error vector (e / ewt) on an error return, -1
otherwise.
'lenrw'  the length of the double work array required.
'leniw'  the length of integer work array required.
'mused'  a vector of method indicators for each successful time step:
1: adams (nonstiff), 2: bdf (stiff)
=======  ============================================================

Other Parameters
----------------
ml, mu : int, optional
If either of these are not None or non-negative, then the
Jacobian is assumed to be banded.  These give the number of
lower and upper non-zero diagonals in this banded matrix.
For the banded case, Dfun should return a matrix whose
rows contain the non-zero bands (starting with the lowest diagonal).
Thus, the return matrix jac from Dfun should have shape
(ml + mu + 1, len(y0)) when ml >=0 or mu >=0.
The data in jac must be stored such that jac[i - j + mu, j]
holds the derivative of the ith equation with respect to the jth
state variable.  If col_deriv is True, the transpose of this
jac must be returned.
rtol, atol : float, optional
The input parameters rtol and atol determine the error
control performed by the solver.  The solver will control the
vector, e, of estimated local errors in y, according to an
inequality of the form max-norm of (e / ewt) <= 1,
where ewt is a vector of positive error weights computed as
ewt = rtol * abs(y) + atol.
rtol and atol can be either vectors the same length as y or scalars.
Defaults to 1.49012e-8.
tcrit : ndarray, optional
Vector of critical points (e.g. singularities) where integration
care should be taken.
h0 : float, (0: solver-determined), optional
The step size to be attempted on the first step.
hmax : float, (0: solver-determined), optional
The maximum absolute step size allowed.
hmin : float, (0: solver-determined), optional
The minimum absolute step size allowed.
ixpr : bool, optional
Whether to generate extra printing at method switches.
mxstep : int, (0: solver-determined), optional
Maximum number of (internally defined) steps allowed for each
integration point in t.
mxhnil : int, (0: solver-determined), optional
Maximum number of messages printed.
mxordn : int, (0: solver-determined), optional
Maximum order to be allowed for the non-stiff (Adams) method.
mxords : int, (0: solver-determined), optional
Maximum order to be allowed for the stiff (BDF) method.

--------
ode : a more object-oriented integrator based on VODE.
quad : for finding the area under a curve.
File:      /usr/local/lib/python3.4/dist-packages/scipy/integrate/odepack.py
Type:      function



To show how odeint works, we will solve the Lotka–Volterra equations, an example of a predator-prey model:

$$\frac{dx}{dt} = \alpha x - \beta x y$$$$\frac{dy}{dt} = \delta x y - \gamma y$$

where:

• $x(t)$ is the number of prey.
• $y(t)$ is the number of predators.
• $\alpha$ is the natural birth rate of the prey.
• $\gamma$ is the natural death rate of the predators.
• $\beta$ determines the death rate of prey when eaten by predators.
• $\delta$ determines the growth rate of predators when they eat prey.

It is important to note here that $y(t)$ is different from the overall solutions vector $\vec{y}(t)$. In fact, perhaps confusingly, in this case $\vec{y}(t)=[x(t),y(t)]$.

To integrate this system of differential equations, we must define a function derivs that computes the right-hand-side of the differential equation, $\vec{f}(\vec{y}(t), t)$. The signature of this function is set by odeint itself:

def derivs(yvec, t, *args):
...
return dyvec

• yvec will be a 1d NumPy array with $N$ elements that are the values of the solution at the current time, $\vec{y}(t)$.
• t will be the current time.
• *args will be other arguments, typically parameters in the differential equation.

The derivs function must return a 1d NumPy array with elements that are the values of the function $\vec{f}(\vec{y}(t), t)$.



In :

def derivs(yvec, t, alpha, beta, delta, gamma):
x = yvec
y = yvec
dx = alpha*x - beta*x*y
dy = delta*x*y - gamma*y
return np.array([dx, dy])



Here are the parameters and initial condition we will use to solve the differential equation. In this case, our prey variable $x$ is the number of rabbits and the predator variable $y$ is the number of foxes (foxes eat rabbits).



In :

nfoxes = 10
nrabbits = 20
ic = np.array([nrabbits, nfoxes])
maxt = 20.0
alpha = 1.0
beta = 0.1
delta = 0.1
gamma = 1.0



Here we call odeint with our derivs function, initial condition ic, array of times t and the extra parameters:



In :

t = np.linspace(0, maxt, int(100*maxt))
soln = odeint(derivs, # function to compute the derivatives
ic,     # array of initial conditions
t,      # array of times
args=(alpha, beta, delta, gamma), # extra args
atol=1e-9, rtol=1e-8)             # absolute and relative error tolerances



We can plot the componenets of the solution as a function of time as follows:



In :

plt.plot(t, soln[:,0], label='rabbits')
plt.plot(t, soln[:,1], label='foxes')
plt.xlabel('t')
plt.ylabel('count')
plt.legend();






We can also make a parametric plot of $[x(t),y(t)]$:



In :

plt.plot(soln[:,0], soln[:,1])
plt.xlim(0, 25)
plt.ylim(0, 25)
plt.xlabel('rabbits')
plt.ylabel('foxes');