scikits-odes¶

scikits_odes¶

scikits-odes is a scikit offering a cython wrapper around some extra ode/dae solvers, so they can mature outside of scipy.

It offers wrappers around the following solvers from SUNDIALS

CVODE
IDA

It additionally offers wrappers around

ddaspk <http://www.netlib.org/ode/ddaspk.f> (included)
lsodi <http://www.netlib.org/ode/lsodi.f> (included)

scikits_odes.ode¶

First-order ODE solver¶

User-friendly interface to various numerical integrators for solving a system of first order ODEs with prescribed initial conditions:

\[ \begin{align}\begin{aligned}\frac{dy(t)}{dt} = f(t, y(t)), \quad y(t_0)=y_0\\y(t_0)[i] = y_0[i], i = 0, ..., \mathrm{len}(y_0) - 1\end{aligned}\end{align} \]

\(f(t,y)\) is the right hand side function and returns a vector of size \(\mathrm{len}(y_0)\).

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 Runge-Kutta order (4)5 solver from scipy.
eqsrhs (right-hand-side function) –
Right-hand-side 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 right-hand-side, 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)

scikits_odes.dae¶

First-order DAE solver¶

User-friendly interface to various numerical integrators for solving an algebraic system of first order ODEs with prescribed initial conditions:

\[ \begin{align}\begin{aligned}A \frac{dy(t)}{dt} = f(t,y(t)),\\y(t=0)[i] = y0[i],\\\frac{d y(t=0)}{dt}[i] = yprime0[i],\end{aligned}\end{align} \]

where \(i = 0, ..., len(y0) - 1\); \(A\) is a (possibly singular) matrix of size \(i × i\); and \(f(t,y)\) is a vector of size \(i\) or more generally, equations of the form

\[G(t,y,y') = 0\]

class scikits_odes.dae.dae(integrator_name, eqsres, **options)[source]¶

A generic interface class to differential algebraic equations.

Define equation res = G(t,y,y’) which can eg be G = f(y,t) - A y’ when solving A y’ = f(y,t), and where (optional) jac is the jacobian matrix of the nonlinear system see fortran source code), so d res/dy + scaling * d res/dy’ or d res/dy depending on the backend.

Parameters:

integrator_name ('ida', 'ddaspk' or 'lsodi') – The integrator solver to use.
eqsres (residual function) –
Residual of the DAE. The signature of this function depends on the solver used, see the solver documentation for details. Generally however, you can assume the following signature to work:

eqsres(x, y, yprime, return_residual)

with x : independent variable, eg the time, float y : array of n unknowns in x yprime : dy/dx array of n unknowns in x, dimension = dim(y) return_residual: array that must be updated with the value of the residuals, so G(t,y,y’). The dimension is equal to dim(y) return value: integer, 0 for success. It is not guaranteed that a solver takes this status into account

Some solvers will allow userdata to be passed to eqsres, or optional formats that are more performant.
options (mapping) – Additional options for initialization, solver dependent See set_options method of the integrator_name you selected for details.

See also

odeint: an ODE integrator with a simpler interface based on lsoda from ODEPACK
ode: class around vode ODE integrator

Notes

Possible future solvers

ddaskr: Not included, starting hints: http://osdir.com/ml/python.f2py.user/2005-07/msg00014.html Modified Extended Backward Differentiation Formulae (MEBDF): Not included. Fortran codes: http://www.ma.ic.ac.uk/~jcash/IVP_software/readme.html

Examples

DAE 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, which we write as G(y,y’,t) = 0 instead of the normal ode, and solve as a DAE.

>>> from __future__ import print_function
>>> from numpy import cos, sin, sqrt
>>> k = 4.0
>>> m = 1.0
>>> initx = [1, 0.1]
>>> initxp = [initx[1], -k/m*initx[0]]
>>> def reseqn(t, x, xdot, result):
    ... # we create residual equations for the problem
    ... result[0] = m*xdot[1] + k*x[0]
    ... result[1] = xdot[0] - x[1]
>>> from scikits.odes import dae
>>> solver = dae('ida', reseqn)
>>> result = solver.solve([0., 1., 2.], initx, initxp)

init_step(t0, y0, yp0, y_ic0_retn=None, yp_ic0_retn=None)[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 (list/array) – initial condition for y
yp0 (list/array) – initial condition for yp
y_ic0 (numpy array) – (optional) returns the calculated consistent initial condition for y
yp_ic0 (numpy array) – (optional) returns the calculated consistent initial condition for y derivated.

Returns:

old_api is False (namedtuple) – namedtuple with the following attributes

Field	Meaning
`flag`	An integer flag (StatusEnumXXX)
`values`	Named tuple with entries t and y and ydot. y will correspond to y_retn value and ydot to yp_retn!
`errors`	Named tuple with entries t and y and ydot
`roots`	Named tuple with entries t and y and ydot
`tstop`	Named tuple with entries t and y and ydot
`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 occurred)

t_out

time, where the solver stopped (when no error occurred, t_out == t)

set_options(**options)[source]¶

Set specific options for the solver. See the solver documentation for details.

Calling set_options a second time, normally resets the solver.

solve(tspan, y0, yp0)[source]¶

Runs the solver.

Parameters:

tspan (list/array) – A list of times at which the computed value will be returned. Must contain the start time as first entry.
y0 (list/array) – list array of initial values
yp0 (list/array) – list array of initial values of derivatives

Returns:

old_api is False (namedtuple) – namedtuple with the following attributes

Field	Meaning
`flag`	An integer flag (StatusEnumXXX)
`values`	Named tuple with entries array t and array y and array ydot. y will correspond to y_retn value and ydot to yp_retn!
`errors`	Named tuple with entries t and y and ydot of error
`roots`	Named tuple with entries array t and array y and array ydot
`tstop`	Named tuple with entries array t and array y and array ydot
`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])
`yp`	numpy array of derivatives corresponding to times t (values of yp[i, :] ~ t[i])
`t_err`	float or None - if recoverable error occurred (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
`yp_err`	numpy array of derivatives corresponding to time t_err

step(t, y_retn=None, yp_retn=None)[source]¶

Method for calling successive next step of the IDA solver to allow more precise control over the IDA 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; else if if t<0.0 only one internal step is performed towards time abs(t) and results after this one time step are returned.

Parameters:

t (number)
y_retn (numpy array (ndim = 1) or None.) – (Needs to be preallocated) If not None, will be filled with y at time t. If None y_retn is not used.
yp_retn (numpy array (ndim = 1) or None.) – (Needs to be preallocated) If not None, will be filled with derivatives of y at time t. If None yp_retn is not used.

Returns:

old_api is False (namedtuple) – namedtuple with the following attributes

Field	Meaning
`flag`	An integer flag (StatusEnumXXX)
`values`	Named tuple with entries t and y and ydot. y will correspond to y_retn value and ydot to yp_retn!
`errors`	Named tuple with entries t and y and ydot
`roots`	Named tuple with entries t and y and ydot
`tstop`	Named tuple with entries t and y and ydot
`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 occurred)

t_out

time, where the solver stopped (when no error occurred, t_out == t)

scikits_odes.dopri5¶

Making scipy ode solvers available via the ode API¶

dopri5¶

This is an explicit runge-kutta method of order (4)5 due to Dormand & Prince (with stepsize control and dense output). The API of this solver is as the other scikit.odes ODE solvers

Authors:

Hairer and G. Wanner Universite de Geneve, Dept. de Mathematiques CH-1211 Geneve 24, Switzerland e-mail: ernst.hairer@math.unige.ch, gerhard.wanner@math.unige.ch

This code is described in [HNW93].

This integrator accepts the following options:

atol : float or sequence absolute tolerance for solution, default 1e-12 rtol : float or sequence relative tolerance for solution, default 1e-6 nsteps : int Maximum number of (internally defined) steps allowed during one call to the solver. Default=500 first_step : float max_step : float safety : float Safety factor on new step selection (default 0.9) ifactor : float dfactor : float Maximum factor to increase/decrease step size by in one step beta : float Beta parameter for stabilised step size control. verbosity : int Switch for printing messages (< 0 for no messages).

References [HNW93] (1, 2) E. Hairer, S.P. Norsett and G. Wanner, Solving Ordinary Differential Equations i. Nonstiff Problems. 2nd edition. Springer Series in Computational Mathematics, Springer-Verlag (1993)

dop853¶

This is an explicit runge-kutta method of order 8(5,3) due to Dormand & Prince (with stepsize control and dense output). Options and references the same as “dopri5”.

exception scikits_odes.dopri5.DOPSolveException(soln)[source]¶: Base class for exceptions raised by DOP.validate_flags.

exception scikits_odes.dopri5.DOPSolveFailed(soln)[source]¶: DOP.solve failed to reach endpoint

class scikits_odes.dopri5.SolverReturn(flag, values, errors, roots, tstop, message)¶

errors¶: Alias for field number 2

flag¶: Alias for field number 0

message¶: Alias for field number 5

roots¶: Alias for field number 3

tstop¶: Alias for field number 4

values¶: Alias for field number 1

class scikits_odes.dopri5.SolverVariables(t, y)¶

t¶: Alias for field number 0

y¶: Alias for field number 1

class scikits_odes.dopri5.StatusEnumDOP(value, names=None, *, module=None, qualname=None, type=None, start=1, boundary=None)[source]¶

class scikits_odes.dopri5.dop853(Rfn, **options)[source]¶

class scikits_odes.dopri5.dopri5(Rfn, **options)[source]¶

init_step(t0, y0)[source]¶

Initializes the solver and allocates memory.

Parameters:

time (t0 - initial)
array) (y0 - initial condition for y (can be list or numpy)

Returns:

if old_api – not supported
if old_api False –

A named tuple, with entries:
flag = An integer flag (StatusEnumDop) values = Named tuple with entries t and y and ydot. y will correspond to y_retn value and ydot to yp_retn! errors = Named tuple with entries t_err and y_err roots = Named tuple with entries t_roots and y_roots tstop = Named tuple with entries t_stop and y_tstop message= String with message in case of an error

set_options(**options)[source]¶

Set specific options for the solver.

Calling set_options a second time, normally resets the solver.

solve(tspan, y0)[source]¶

Runs the solver.

Parameters:

be (tspan - an list/array of times at which the computed value will) – returned. Must contain the start time as first entry..
values (y0 - list/numpy array of initial)

Returns:

if old_api – Not supported
if old_api False –

A named tuple, with entries:
flag = An integer flag 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

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.

Parameters:

t (t - A step is done towards time) –
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
returned. (and output at t) –
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
computed (y_retn - numpy vector (ndim = 1) in which the) – value will be stored (needs to be preallocated). If None y_retn is not used.

Returns:

if old_api – not supported
if old_api False –

A named tuple, with entries:
flag = An integer flag (StatusEnumDOP) values = Named tuple with entries t and y. y will correspond to y_retn value errors = Named tuple with entries t_err and y_err roots = Named tuple with entries t_roots and y_roots tstop = Named tuple with entries t_stop and y_tstop message= String with message in case of an error

validate_flags(soln)[source]¶

Validates the flag returned by dopri.solve.

Validation happens using the following scheme:

failures (flag < 0) raise DOPSolveFailed or a subclass of it;
otherwise, return an instance of SolverReturn.

scikits-odes¶

scikits_odes¶

scikits_odes.ode¶

First-order ODE solver¶

scikits_odes.dae¶

First-order DAE solver¶

scikits_odes.dopri5¶

Making scipy ode solvers available via the ode API¶

dopri5¶

dop853¶

Odes

Navigation

Related Topics

Field	Meaning
`flag`	boolean status of the computation (successful or error occured)
`t_out`	inititial time

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)