ODE Solvers¶
scikits.odes
contains two main routines for solving ODEs: the simpler
scikits.odes.odeint.odeint()
, and the more configurable
scikits.odes.ode.ode
. Both these routines allow selection of the
solver and solution method used. Additionally, it is also possible to directly
use the low level interface to individual solvers.

scikits.odes.odeint.
odeint
(rhsfun, tout, y0, method='bdf', **options)[source]¶ Integrate a system of ordinary differential equations. odeint is a wrapper around the ode class, as a convenience function to quickly integrate a system of ode. Solves the initial value problem for stiff or nonstiff systems of first order ode’s:
rhs = dy/dt = fun(t, y)
where y can be a vector, then rhsfun must be a function computing rhs with signature:
rhsfun(t, y, rhs)
storing the computed dy/dt in the rhs array passed to the function
 Parameters
rhsfun (callable(t, y, out)) – Computes the derivative at dy/dt in terms of t and y, and stores in out
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.
method (string, solution method to use.) –
Available are all the ode class solvers as well as some convenience shorthands:
Method
Meaning
bdf
This uses the ‘cvode’ solver in default from, which is a variable step, variable coefficient Backward Differentiation Formula solver, good for stiff ODE. Newton iterations are used to solve the nonlinear system.
admo
This uses the ‘cvode’ solver with option lmm_type=’ADAMS’, which is a variable step AdamsMoulton method (linear multistep method), good for nonstiff ODE. Functional iterations are used to solve the nonlinear system.
rk5
This uses the ‘dopri5’ solver, which is a variable step RungeKutta method of order (4)5 (use for nonstiff ODE)
rk8
This uses the ‘dop853’ solver, which is a variable step RungeKutta method of order 8(5,3)
For educational purposes, you can also access following methods:
Method
Meaning
beuler
This is the Implicit Euler (backward Euler) method (order 1), which is obtained via the ‘bdf’ method, setting the order option to 1, setting large tolerances, and fixing the stepsize. Use option ‘step’ to change stepsize, default: step=0.05. Use option ‘rtol’ and ‘atol’ to use more strict tolerances Note: this is not completely the backward Euler method, as the cvode solver has added control options!
trapz
This is the Trapezoidal Rule method (order 2), which is obtained via the ‘admo’ method, setting option order to 2, setting large tolerances and fixing the stepsize. Use option ‘step’ to change stepsize, default: step=0.05. Use option ‘rtol’ and ‘atol’ to use more strict tolerances Note: The cvode solver might change the order to 1 internally in which case this becomes beuler method. Set atol, rtol options as strict as possible.
You can also access the solvers of ode via their names:
Method
Meaning
cvode
This uses the ‘cvode’ solver
dopri5
This uses the ‘dopri5’ solver
dop853
This uses the ‘dop853’ solver
options (extra solver options, optional) – Every solver has it’s own extra options, see the ode class and the details of the solvers available there to know the options possible per solver
 Returns
solution – A single named tuple is returned containing the result of the integration.
Field
Meaning
flag
An integer flag
values
Named tuple with fields t and y
errors
Named tuple with fields t and y
roots
Named tuple with fields t and y
tstop
Named tuple with fields t and y
message
String with message in case of an error
 Return type
named tuple
See also
scikits.odes.ode.ode()
a more objectoriented integrator
scikits.odes.dae.dae()
a solver for differentialalgebraic equations
scipy.integrate.quad()
for finding the area under a curve
Examples
The second order differential equation for the angle theta of a pendulum acted on by gravity with friction can be written:
\[\theta''(t) + b \theta'(t) + c \sin(\theta(t)) = 0\]where b and c are positive constants, and a prime (‘) denotes a derivative. To solve this equation with odeint, we must first convert it to a system of first order equations. By defining the angular velocity
omega(t) = theta'(t)
, we obtain the system:\[\theta'(t) = \omega(t) \omega'(t) = b \omega(t)  c \sin(\theta(t))\]We assume the constants are b = 0.25 and c = 5.0:
>>> b = 0.25 >>> c = 5.0
Let y be the vector [theta, omega]. We implement this system in python as:
>>> def pend(t, y, out): ... theta, omega = y ... out[:] = [omega, b*omega  c*np.sin(theta)] ...
In case you want b and c easily changable, make pend a class method, and consider attributes b and c accessible via self.b and self.c. For initial conditions, we assume the pendulum is nearly vertical with theta(0) = pi  0.1, and it initially at rest, so omega(0) = 0. Then the vector of initial conditions is
>>> y0 = [np.pi  0.1, 0.0]
We generate a solution 101 evenly spaced samples in the interval 0 <= t <= 10. So our array of times is
>>> t = np.linspace(0, 10, 101)
Call odeint to generate the solution.
>>> from scikits.odes.odeint import odeint >>> sol = odeint(pend, t, y0)
The solution is a named tuple sol. sol.values.y is an array with shape (101, 2). The first column is theta(t), and the second is omega(t). The following code plots both components.
>>> import matplotlib.pyplot as plt >>> plt.plot(sol.values.t, sol.values.y[:, 0], 'b', label='theta(t)') >>> plt.plot(sol.values.t, sol.values.y[:, 1], 'g', label='omega(t)') >>> plt.legend(loc='best') >>> plt.xlabel('t') >>> plt.grid() >>> plt.show()

class
scikits.odes.ode.
ode
(integrator_name, eqsrhs, **options)[source]¶ A generic interface class to differential equation solvers.
 Parameters
integrator_name (
'cvode'
,'dopri5'
or'dop853'
) –The solver to use.
'cvode'
selects the CVODE solver from the SUNDIALS package. See py:class:scikits.odes.sundials.cvode.CVODE for solver specific options.'dopri5'
selects the Dormand & Prince RungeKutta order (4)5 solver from scipy.eqsrhs (righthandside function) –
Righthandside of a first order ode. Generally, you can assume the following signature to work:
eqsrhs(x, y, return_rhs)
with
x: independent variable, eg the time, float
y: array of n unknowns in x
return_rhs : array that must be updated with the value of the righthandside, so f(t,y). The dimension is equal to dim(y)
 return value: An integer, 0 for success, 1 for failure.
It is not guaranteed that a solver takes this status into account
Some solvers will allow userdata to be passed to eqsrhs, or optional formats that are more performant.
options (additional options of the solver) – See set_options method of the integrator_name you selected for details. Set option old_api=False to use the new API. In the future, this will become the default!
See also
scikits.odes.odeint.odeint
an ODE integrator with a simpler interface
scipy.integrate
Methods in scipy for ODE integration
Examples
ODE arise in many applications of dynamical systems, as well as in discritisations of PDE (eg moving mesh combined with method of lines). As an easy example, consider the simple oscillator,
>>> from __future__ import print_function >>> from numpy import cos, sin, sqrt >>> k = 4.0 >>> m = 1.0 >>> initx = [1, 0.1] >>> def rhseqn(t, x, xdot): # we create rhs equations for the problem xdot[0] = x[1] xdot[1] =  k/m * x[0]
>>> from scikits.odes import ode >>> solver = ode('cvode', rhseqn, old_api=False) >>> result = solver.solve([0., 1., 2.], initx) >>> print(' t Solution Exact') >>> print('') >>> for t, u in zip(result.values.t, result.values.y): print('%4.2f %15.6g %15.6g' % (t, u[0], initx[0]*cos(sqrt(k/m)*t)+initx[1]*sin(sqrt(k/m)*t)/sqrt(k/m)))
More examples in the Examples directory and IPython worksheets.

get_info
()[source]¶ Return additional information about the state of the integrator.
 Returns
A dictionary filled with internal data as exposed by the integrator.
See the get_info method of your chosen integrator for details.

init_step
(t0, y0)[source]¶ Initializes the solver and allocates memory. It is not needed to call this method if solve is used to compute the solution. In the case step is used, init_step must be called first.
 Parameters
t0 (number) – initial time
y0 (array) – initial condition for y (can be list or numpy array)
 Returns
old_api is False (namedtuple) – namedtuple with the following attributes
Field
Meaning
flag
An integer flag (StatusEnum)
values
Named tuple with entries t and y
errors
Named tuple with entries t and y
roots
Named tuple with entries t and y
tstop
Named tuple with entries t and y
message
String with message in case of an error
old_api is True (tuple) – tuple with the following elements in order
Field
Meaning
flag
boolean status of the computation (successful or error occured)
t_out
inititial time

set_options
(**options)[source]¶ Set specific options for the solver. See the solver documentation for details.
Calling set_options a second time, is only possible for options that can change during runtime.

set_tstop
(tstop)[source]¶ Add a stop time to the integrator past which he is not allowed to integrate.
 Parameters
tstop (float time) – Time point in the future where the integration must stop. You can indicate like this that integration past this point is not allowed, in order to avoid undefined behavior. You can achieve the same result with a call to set_options(tstop=tstop)

solve
(tspan, y0)[source]¶ Runs the solver.
 Parameters
tspan (array (or similar)) – a list of times at which the computed value will be returned. Must contain the start time as first entry.
y0 (array (or similar)) – a list of initial values
 Returns
old_api is False (namedtuple) – namedtuple with the following attributes
Field
Meaning
flag
An integer flag (StatusEnum)
values
Named tuple with entries t and y
errors
Named tuple with entries t and y
roots
Named tuple with entries t and y
tstop
Named tuple with entries t and y
message
String with message in case of an error
old_api is True (tuple) – tuple with the following elements in order
Field
Meaning
flag
indicating return status of the solver
t
numpy array of times at which the computations were successful
y
numpy array of values corresponding to times t (values of y[i, :] ~ t[i])
t_err
float or None  if recoverable error occured (for example reached maximum number of allowed iterations), this is the time at which it happened
y_err
numpy array of values corresponding to time t_err

step
(t, y_retn=None)[source]¶ Method for calling successive next step of the ODE solver to allow more precise control over the solver. The ‘init_step’ method has to be called before the ‘step’ method.
A step is done towards time t, and output at t returned. This time can be higher or lower than the previous time. If option ‘one_step_compute’==True, and the solver supports it, only one internal solver step is done in the direction of t starting at the current step.
If old_api=True, the old behavior is used: if t>0.0 then integration is performed until this time and results at this time are returned in y_retn if t<0.0 only one internal step is perfomed towards time abs(t) and results after this one time step are returned
 Parameters
t (number) –
y_retn (numpy vector (ndim = 1)) – in which the computed value will be stored (needs to be preallocated). If None y_retn is not used.
 Returns
old_api is False (namedtuple) – namedtuple with the following attributes
Field
Meaning
flag
An integer flag (StatusEnum)
values
Named tuple with entries t and y
errors
Named tuple with entries t and y
roots
Named tuple with entries t and y
tstop
Named tuple with entries t and y
message
String with message in case of an error
old_api is True (tuple) – tuple with the following elements in order
Field
Meaning
flag
status of the computation (successful or error occured)
t_out
time, where the solver stopped (when no error occured, t_out == t)