10.6: Integrating ODEs with Scipy
Except for educational purposes, it is almost always a bad idea to implement your own ODE solver; instead, you should use a pre-written solver.
10.6.1 The
scipy.integrate.odeint
Solver
In Scipy, the simplest ODE solver to use is the
scipy.integrate.odeint
function, which is in the
scipy.integrate
module. This is actually a wrapper around a low-level numerical library known as LSODE (the
L
ivermore
S
olver for
ODE
s"), which is part of a widely-used ODE solver library known as
ODEPACK
. The most important feature of this solver is that it is "adaptive": it can automatically figure out (i) which integration scheme to use (choosing between either a high-order Adams-Moulton method, or another implicit method known as the
Backward Differentiation Formula
which we haven't described), and (ii) the size of the discrete time steps, based on the behavior of the solutions as they are being worked out. In other words, the user only needs to specify the derivative function, the initial state, and the desired output times, without having to worry about the internal details of the solution method.
The function takes several inputs, of which the most important ones are:
-
func
, a function corresponding to the derivative function \(\vec{F}(\vec{y}, t)\). -
y0
, either a number or 1D array, corresponding to the initial state \(\vec{y}(t_0)\). -
t
, an array of times at which to output the ODE solution. The first element corresponding to the initial time \(t_{0}\). Note that these are the "output" times only—they do not specify the actual time steps which the solver uses for finding the solutions; those are automatically determined by the solver. -
(optional)
args
, a tuple of extra inputs to pass to the derivative functionfunc
. For example, ifargs=(2,3)
, thenfunc
should accept four inputs, and it will be passed 2 and 3 as the last two inputs.
The function then returns an array
y
, where
y[n]
contains the solution at time
t[n]
. Note that
y[0]
will be exactly the same as the input
y0
, the initial state which you specified.
Here is an example of using
odeint
to solve the damped harmonic oscillator problem \(m \ddot{x} = - \lambda \dot{x} - k x(t)\), using the previously-mentioned vectorization trick to cast it into a first-order ODE:
from scipy import * import matplotlib.pyplot as plt from scipy.integrate import odeint def ydot(y, t, m, lambd, k): x, v = y[0], y[1] return array([v, -(lambd/m) * v - k * x / m]) m, lambd, k = 1.0, 0.1, 1.0 # Oscillator parameters y0 = array([1.0, 5.0]) # Initial conditions [x, v] t = linspace(0.0, 50.0, 100) # Output times y = odeint(ydot, y0, t, args=(m, lambd, k)) ## Plot x versus t plt.plot(t, y[:,0], 'b-') plt.xlabel('t') plt.ylabel('x') plt.show()
There is an important limitation of
odeint
: it
does not handle complex ODEs
, and always assumes that \(\vec{y}\) and \(\vec{F}\) are real. However, this is not a problem in practice, because you can always convert a complex first-order ODE into a real one, by replacing the complex vectors \(\vec{y}\) and \(\vec{F}\) with double-length real vectors:
\[\vec{y}' \equiv \begin{bmatrix}\mathrm{Re}(\vec{y})\\ \mathrm{Im}(\vec{y})\end{bmatrix}, \;\; \vec{F}' \equiv \begin{bmatrix}\mathrm{Re}(\vec{F})\\ \mathrm{Im}(\vec{F})\end{bmatrix}.\]
10.6.2 The
scipy.integrate.ode
Solvers
Apart from
odeint
, Scipy provides a more general interface to a variety of ODE solvers, in the form of the
scipy.integrate.ode
class. This is a much more low-level interface; instead of calling a single function, you have to create an ODE "object", then use the methods of this object to specify the type of ODE solver to use, the initial conditions, etc.; then you have to repeatedly call the ODE object's
integrate
method, to integrate the solution up to each desired output time step.
The is an extremely aggravating inconsistency between the
odeint
function and this
ode
class: the expected order of inputs for the derivative functions are reversed! The
odeint
function assumes the derivative function has the form
F(y,t)
, but the
ode
class assumes it has the form
F(t,y)
. Watch out for this!
Here is an example of using
ode
class with the damped harmonic oscillator problem \(m \ddot{x} = - \lambda \dot{x} - k x(t)\), using a Runge-Kutta solver:
from scipy import * import matplotlib.pyplot as plt from scipy.integrate import ode ## Note the order of inputs (different from odeint)! def ydot(t, y, m, lambd, k): x, v = y[0], y[1] return array([v, -(lambd/m) * v - k * x / m]) m, lambd, k = 1.0, 0.1, 1.0 # Oscillator parameters y0 = array([1.0, 5.0]) # Initial conditions [x, v] t = linspace(0.0, 50.0, 100) # Output times ## Set up the ODE object r = ode(ydot) r.set_integrator('dopri5') # A Runge-Kutta solver r.set_initial_value(y0) r.set_f_params(m, lambd, k) ## Perform the integration. Note that the "integrate" method only integrates ## up to one single final time point, rather than an array of times. x = zeros(len(t)) x[0] = y0[0] for n in range(1,len(t)): r.integrate(t[n]) assert r.successful() x[n] = (r.y)[0] ## Plot x versus t plt.plot(t, x, 'b-') plt.xlabel('t') plt.ylabel('x') plt.show()
See the documentation for a more detailed list of options, including the list of ODE solvers that you can choose from.