Ordinary Differential Equations

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

Imports


In [1]:
%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 [2]:
tmax = 10.0 # The max time
M = 100     # Use 100 times between [0,tmax]
t = np.linspace(0,tmax,M)
t


Out[2]:
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 [3]:
h = t[1]-t[0]
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 [61]:
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 [16]:
from scipy.integrate import odeint

In [17]:
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
    Dictionary containing additional output information

    =======  ============================================================
    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 `i`th equation with respect to the `j`th
    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.

See Also
--------
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 [62]:
def derivs(yvec, t, alpha, beta, delta, gamma):
    x = yvec[0]
    y = yvec[1]
    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 [57]:
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 [60]:
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 [59]:
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 [67]:
plt.plot(soln[:,0], soln[:,1])
plt.xlim(0, 25)
plt.ylim(0, 25)
plt.xlabel('rabbits')
plt.ylabel('foxes');