# 10.6: Integrating ODEs with Scipy

- Page ID
- 34862

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

`scipy.integrate.odeint`

SolverIn 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 function`func`

. For example, if`args=(2,3)`

, then`func`

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 parametersy0 = array([1.0, 5.0])# Initial conditions [x, v]t = linspace(0.0, 50.0, 100)# Output timesy = odeint(ydot, y0, t, args=(m, lambd, k))## Plot x versus tplt.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

`scipy.integrate.ode`

SolversApart 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 parametersy0 = array([1.0, 5.0])# Initial conditions [x, v]t = linspace(0.0, 50.0, 100)# Output times## Set up the ODE objectr = ode(ydot)r.set_integrator('dopri5')# A Runge-Kutta solverr.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 tplt.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.