diff --git a/v3.0.0a3/_modules/index.html b/v3.0.0a3/_modules/index.html new file mode 100644 index 0000000..03e300a --- /dev/null +++ b/v3.0.0a3/_modules/index.html @@ -0,0 +1,112 @@ + + + +
+ + +
+# -*- coding: utf-8 -*-
+# Authors: B. Malengier based on ode.py
+r"""
+First-order DAE solver
+======================
+
+User-friendly interface to various numerical integrators for solving an
+algebraic system of first order ODEs with prescribed initial conditions:
+
+.. math::
+ A \frac{dy(t)}{dt} = f(t,y(t)),
+
+ y(t=0)[i] = y0[i],
+
+ \frac{d y(t=0)}{dt}[i] = yprime0[i],
+
+where :math:`i = 0, ..., len(y0) - 1`; :math:`A` is a (possibly singular) matrix
+of size :math:`i × i`; and :math:`f(t,y)` is a vector of size :math:`i` or more generally, equations of the form
+
+.. math::
+ G(t,y,y') = 0
+
+"""
+__all__ = ['dae']
+__docformat__ = "restructuredtext en"
+import re
+import sys
+
+from scikits_odes_core import DaeBase
+
+
+#------------------------------------------------------------------------------
+# User interface
+#------------------------------------------------------------------------------
+
+
+[docs]
+class dae(object):
+ """
+ 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)
+ """
+
+ LOADED = False
+
+ def __init__(self, integrator_name, eqsres, **options):
+ integrator = find_dae_integrator(integrator_name)
+ if integrator is None:
+ raise ValueError('No integrator name match with %s or is not available.'\
+ %(repr(integrator_name)))
+ else:
+ self._integrator = integrator(eqsres, **options)
+
+
+[docs]
+ def set_options(self, **options):
+ """
+ Set specific options for the solver.
+ See the solver documentation for details.
+
+ Calling set_options a second time, normally resets the solver.
+ """
+ return self._integrator.set_options(**options)
+
+
+
+[docs]
+ def solve(self, tspan, y0, yp0):
+ """
+ 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
+ ========== ==========================================
+
+ """
+ return self._integrator.solve(tspan, y0, yp0)
+
+
+
+[docs]
+ def init_step(self, t0, y0, yp0, y_ic0_retn = None, yp_ic0_retn = None):
+ """
+ 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)
+ =========== ==========================================
+
+ """
+ return self._integrator.init_step(t0, y0, yp0, y_ic0_retn, yp_ic0_retn)
+
+
+
+[docs]
+ def step(self, t, y_retn=None, yp_retn=None):
+ """
+ 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)
+ =========== ==========================================
+
+ """
+ return self._integrator.step(t, y_retn, yp_retn)
+
+
+ def __del__(self):
+ """
+ Clean up what is needed
+ """
+ if hasattr(self, '_integrator'):
+ del self._integrator
+
+
+#------------------------------------------------------------------------------
+# DAE integrators
+#------------------------------------------------------------------------------
+
+def find_dae_integrator(name):
+ if not dae.LOADED:
+ ## ida
+ try:
+ from scikits_odes_sundials import ida
+ DaeBase.integrator_classes.append(ida.IDA)
+ except ValueError as msg:
+ print('Could not load IDA solver', msg)
+ except ImportError:
+ print(sys.exc_info()[1])
+
+ ## idas
+ try:
+ from scikits_odes_sundials import idas
+ DaeBase.integrator_classes.append(idas.IDAS)
+ except ValueError as msg:
+ print('Could not load IDAS solver', msg)
+ except ImportError:
+ print(sys.exc_info()[1])
+
+ ## ddaspk
+ try:
+ from scikits_odes_daepack.ddaspkint import ddaspk
+ except ImportError:
+ print(sys.exc_info()[1])
+
+ ## lsodi
+ try:
+ from scikits_odes_daepack.lsodiint import lsodi
+ except ImportError:
+ print(sys.exc_info()[1])
+
+ dae.LOADED = True
+
+ for cl in DaeBase.integrator_classes:
+ if re.match(name, cl.__name__, re.I):
+ return cl
+ elif hasattr(cl, name) and re.match(name, cl.name, re.I):
+ return cl
+ raise ValueError('Integrator name %s does not exsist' % name)
+
+# -*- coding: utf-8 -*-
+# Created on Thu Jan 28 14:59:25 2016
+# @author: benny
+"""
+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:
+
+ E. 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”.
+
+"""
+from collections import namedtuple
+from enum import IntEnum
+from warnings import warn
+import sys
+
+import numpy as np
+from scipy.integrate import ode as _runner
+
+from scikits_odes_core import OdeBase
+
+SolverReturn = namedtuple(
+ "SolverReturn", [
+ "flag", "values", "errors", "roots",
+ "tstop", "message"
+ ]
+)
+
+SolverVariables = namedtuple("SolverVariables", ["t", "y"])
+
+
+
+[docs]
+class StatusEnumDOP(IntEnum):
+ SUCCESS = 1
+ SOLOUT = 2
+
+ INPUT_FAIL = -1
+ NMAX_FAIL = -2
+ STEPSIZE_FAIL = -3
+ STIFFNESS_FAIL = -4
+
+ UNEXPECTED_IDID = -10
+
+
+STATUS_MESSAGE = {
+ StatusEnumDOP.SUCCESS: 'computation successful',
+ StatusEnumDOP.SOLOUT: 'comput. successful (interrupted by solout)',
+ StatusEnumDOP.INPUT_FAIL: 'input is not consistent',
+ StatusEnumDOP.NMAX_FAIL: 'larger nmax is needed',
+ StatusEnumDOP.STEPSIZE_FAIL: 'step size becomes too small',
+ StatusEnumDOP.STIFFNESS_FAIL: 'problem is probably stiff (interrupted)',
+ StatusEnumDOP.UNEXPECTED_IDID: 'Unexpected idid, check warnings for info',
+}
+
+
+[docs]
+class DOPSolveException(Exception):
+ """Base class for exceptions raised by `DOP.validate_flags`."""
+ def __init__(self, soln):
+ self.soln = soln
+ self.args = (self._message.format(soln),)
+
+
+
+[docs]
+class DOPSolveFailed(DOPSolveException):
+ """`DOP.solve` failed to reach endpoint"""
+ _message = (
+ "Solver failed with flag {0.flag} and finished at {0.errors.t}"
+ "with values {0.errors.y}."
+ )
+
+
+
+[docs]
+class dopri5(OdeBase):
+ name = 'dopri5'
+ _runner = _runner
+ default_values = {
+ 'rtol': 1e-6,
+ 'atol': 1e-12,
+ 'nsteps': 500,
+ 'max_step': 0.0,
+ 'first_step': 0.0, # determined by solver
+ 'safety': 0.9,
+ 'ifactor': 10.0,
+ 'dfactor': 0.2,
+ 'beta': 0.0,
+ 'verbosity': -1, # no messages if negative
+ }
+
+ def __init__(self, Rfn, **options):
+ """
+ Initialize the ODE Solver and it's default values
+
+ Parameters
+ ----------
+ Rfn - right-hand-side function
+ options - additional options for initialization
+ """
+ self.options = self.default_values
+ self.set_options(rfn=Rfn, **options)
+ self._validate_flags = None
+ self.solver = None
+ self.initialized = False
+
+
+[docs]
+ def set_options(self, **options):
+ """
+ Set specific options for the solver.
+
+ Calling set_options a second time, normally resets the solver.
+ """
+ for (key, value) in options.items():
+ self.options[key.lower()] = value
+ self.initialized = False
+
+
+ def _init_data(self):
+ self.rfn = self.options['rfn']
+ self.N = len(self.y0)
+
+ def _wrap_Rfn(self, t, y):
+ yout = np.empty(self.N, float)
+ self.rfn(t, y, yout)
+ return yout
+
+
+[docs]
+ def init_step(self, t0, y0):
+ """
+ Initializes the solver and allocates memory.
+
+ Parameters
+ ----------
+ t0 - initial time
+ y0 - initial condition for y (can be list or numpy array)
+
+ 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
+
+ """
+ self.y0 = y0
+ self.t = t0
+ self._init_data()
+
+ self.solver = self._runner(self._wrap_Rfn)\
+ .set_integrator(self.name,
+ rtol = self.options['rtol'],
+ atol = self.options['atol'],
+ nsteps = self.options['nsteps'],
+ max_step = self.options['max_step'],
+ first_step = self.options['first_step'], # determined by solver
+ safety = self.options['safety'],
+ ifactor = self.options['ifactor'],
+ dfactor = self.options['dfactor'],
+ beta = self.options['beta'],
+ verbosity = self.options['verbosity'])\
+ .set_initial_value(y0, t0)
+ self.initialized = True
+ y_retn = np.empty(len(y0), float)
+ y_retn[:] = y0[:]
+ soln = SolverReturn(
+ flag=StatusEnumDOP.SUCCESS,
+ values=SolverVariables(t=t0, y=y_retn),
+ errors=SolverVariables(t=None, y=None),
+ roots=SolverVariables(t=None, y=None),
+ tstop=SolverVariables(t=None, y=None),
+ message=STATUS_MESSAGE[StatusEnumDOP.SUCCESS]
+ )
+ if self._validate_flags:
+ return self.validate_flags(soln)
+ return soln
+
+
+
+[docs]
+ def step(self, t, y_retn=None):
+ """
+ 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 - 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
+ 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
+ -------
+ 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
+
+ """
+ if not self.initialized == False:
+ raise ValueError("DOPRI:step: init_step must be run before running the step method.")
+ if t <= self.t:
+ raise ValueError("Integration must be forward! t must be > previous timestep")
+ y = self.solver.integrate(t)
+ y_err = None
+ t_err = None
+ if not self.solver.successful():
+ flag = StatusEnumDOP.UNEXPECTED_IDID
+ y_err = y
+ t_err = t
+ else:
+ flag = StatusEnumDOP.SUCCESS
+ self.y = y
+ self.t = self.solver.t
+
+ return SolverReturn(
+ flag=flag,
+ values=SolverVariables(t=self.t, y=y),
+ errors=SolverVariables(t=t_err, y=y_err),
+ roots=SolverVariables(t=None, y=None),
+ tstop=SolverVariables(t=None, y=None),
+ message=STATUS_MESSAGE[flag]
+ )
+
+
+
+[docs]
+ def solve(self, tspan, y0):
+ """
+ Runs the solver.
+
+ Parameters
+ ----------
+ tspan - an list/array of times at which the computed value will be
+ returned. Must contain the start time as first entry..
+ y0 - list/numpy array of initial values
+
+ 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
+ """
+ self.initialized = False
+ soln = self.init_step(tspan[0], y0)
+ nrt = len(tspan)
+ t_retn = np.empty(np.shape(tspan), float)
+ y_retn = np.empty([len(tspan), len(y0)], float)
+
+ t_retn[0] = tspan[0]
+ y_retn[0, :] = y0
+
+ idx = 1
+ while self.solver.successful() and idx<nrt:
+ time = tspan[idx]
+ y = self.solver.integrate(time)
+ t_retn[idx] = time
+ y_retn[idx, :] = y
+ idx += 1
+
+ y_err = None
+ t_err = None
+ if not self.solver.successful():
+ flag = StatusEnumDOP.UNEXPECTED_IDID
+ y_err = y_retn[idx-1]
+ t_err = t_retn[idx-1]
+ # return values computed so far
+ t_retn = t_retn[0:idx-1]
+ y_retn = y_retn[0:idx-1, :]
+ else:
+ flag = StatusEnumDOP.SUCCESS
+
+ soln = SolverReturn(
+ flag=flag,
+ values=SolverVariables(t=t_retn, y=y_retn),
+ errors=SolverVariables(t=t_err, y=y_err),
+ roots=SolverVariables(t=None, y=None),
+ tstop=SolverVariables(t=None, y=None),
+ message=STATUS_MESSAGE[flag]
+ )
+ if self._validate_flags:
+ return self.validate_flags(soln)
+ return soln
+
+
+
+[docs]
+ def validate_flags(self, soln):
+ """
+ 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`.
+
+ """
+ if soln.flag == StatusEnumDOP.SUCCESS:
+ return soln
+ if soln.flag < 0:
+ raise DOPSolveFailed(soln)
+ warn(WARNING_STR.format(soln.flag, *soln.err_values))
+ return soln
+
+
+
+
+
+[docs]
+class dop853(dopri5):
+ name = 'dop853'
+ default_values = {
+ 'rtol': 1e-6,
+ 'atol': 1e-12,
+ 'nsteps': 500,
+ 'max_step': 0.0,
+ 'first_step': 0.0, # determined by solver
+ 'safety': 0.9,
+ 'ifactor': 6.0,
+ 'dfactor': 0.3,
+ 'beta': 0.0,
+ 'verbosity': -1, # no messages if negative
+ }
+
+
+
+OdeBase.integrator_classes.append(dopri5)
+OdeBase.integrator_classes.append(dop853)
+
+# -*- coding: utf-8 -*-
+# Authors: B. Malengier based on ode.py
+r"""
+First-order ODE solver
+======================
+
+User-friendly interface to various numerical integrators for solving a
+system of first order ODEs with prescribed initial conditions:
+
+.. math::
+
+ \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
+
+:math:`f(t,y)` is the right hand side function and returns a vector of size
+:math:`\mathrm{len}(y_0)`.
+
+"""
+__all__ = ['ode']
+__docformat__ = "restructuredtext en"
+
+import re
+import sys
+
+
+from scikits_odes_core import OdeBase
+
+#------------------------------------------------------------------------------
+# User interface
+#------------------------------------------------------------------------------
+
+
+[docs]
+class ode(object):
+ """
+ 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.
+
+ .. _Examples: https://github.com/bmcage/odes/tree/master/docs/src/examples
+ .. _IPython: https://github.com/bmcage/odes/tree/master/docs/ipython
+ """
+ LOADED = False
+
+ def __init__(self, integrator_name, eqsrhs, **options):
+ integrator = find_ode_integrator(integrator_name)
+ if integrator is None:
+ raise ValueError('No integrator name match with %s or is not available.'\
+ %(repr(integrator_name)))
+ else:
+ self._integrator = integrator(eqsrhs, **options)
+
+
+[docs]
+ def set_options(self, **options):
+ """
+ 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.
+ """
+ return self._integrator.set_options(**options)
+
+
+
+[docs]
+ def solve(self, tspan, y0):
+ """
+ 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
+ ========== ==========================================
+
+ """
+ return self._integrator.solve(tspan, y0)
+
+
+
+[docs]
+ def init_step(self, t0, y0):
+ """
+ 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
+ ========== ==========================================
+
+ """
+ return self._integrator.init_step(t0, y0)
+
+
+
+[docs]
+ def step(self, t, y_retn=None):
+ """
+ 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)
+ ========== ==========================================
+
+ """
+ return self._integrator.step(t, y_retn)
+
+
+
+[docs]
+ def set_tstop(self, tstop):
+ """
+ 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)`
+ """
+ if hasattr(self._integrator, 'set_tstop'):
+ self._integrator.set_tstop(tcrit)
+ else:
+ self._integrator.set_options(tstop=tstop)
+
+
+
+[docs]
+ def get_info(self):
+ """
+ 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.
+
+ """
+ if hasattr(self._integrator, 'get_info'):
+ return self._integrator.get_info()
+ else:
+ return {}
+
+
+
+#------------------------------------------------------------------------------
+# ODE integrators
+#------------------------------------------------------------------------------
+
+def find_ode_integrator(name):
+ if not ode.LOADED:
+ ## cvode
+ try:
+ from scikits_odes_sundials import cvode
+ OdeBase.integrator_classes.append(cvode.CVODE)
+ except ValueError as msg:
+ print('Could not load CVODE solver', msg)
+ except ImportError:
+ print(sys.exc_info()[1])
+ ## cvodes
+ try:
+ from scikits_odes_sundials import cvodes
+ OdeBase.integrator_classes.append(cvodes.CVODES)
+ except ValueError as msg:
+ print('Could not load CVODES solver', msg)
+ except ImportError:
+ print(sys.exc_info()[1])
+
+ ## dopri5 and dop853
+ try:
+ from .dopri5 import dopri5, dop853
+ except ImportError:
+ print(sys.exc_info()[1])
+
+ ode.LOADED = True
+
+ for cl in OdeBase.integrator_classes:
+ if re.match(name, cl.__name__, re.I):
+ return cl
+ elif hasattr(cl, name) and re.match(name, cl.name, re.I):
+ return cl
+ raise ValueError('Integrator name %s does not exist' % name)
+
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Feb 1 13:38:24 2016
+
+@author: benny
+"""
+from copy import copy
+from .ode import ode
+
+
+
+[docs]
+def odeint(rhsfun, tout, y0, method='bdf', **options):
+ """
+ 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 non-stiff 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 Adams-Moulton method (linear
+ multistep method), good for non-stiff ODE. Functional
+ iterations are used to solve the nonlinear system.
+ rk5 This uses the 'dopri5' solver, which is a variable step
+ Runge-Kutta method of order (4)5 (use for non-stiff ODE)
+ rk8 This uses the 'dop853' solver, which is a variable step
+ Runge-Kutta 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 : named tuple
+ 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
+ ======== ==========================================
+
+ See Also
+ --------
+ scikits.odes.ode.ode : a more object-oriented integrator
+ scikits.odes.dae.dae : a solver for differential-algebraic 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:
+
+ .. math:: \\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:
+
+ .. math:: \\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()
+ """
+ t = copy(tout)
+ y0 = copy(y0)
+
+ int_name = 'cvode'
+ #always use old api
+ options['old_api'] = False
+ if (method == 'bdf'):
+ int_name = 'cvode'
+ elif (method == 'admo'):
+ options['lmm_type'] = 'ADAMS'
+ options['nonlinsolver'] = 'fixedpoint'
+ int_name = 'cvode'
+ elif (method == 'rk5'):
+ int_name = 'dopri5'
+ elif (method == 'rk8'):
+ int_name = 'dop853'
+ elif (method == 'beuler' or method == 'trapz'):
+ int_name = 'cvode'
+ options['order'] = 1
+ if not options.has_key('step'):
+ options['step'] = 0.05
+ options['min_step_size'] = options['step']
+ options['max_step_size'] = options['step']
+ #no error control possible, try to disable it:
+ if not options.has_key('atol'):
+ options['atol'] = 1e3
+ if not options.has_key('rtol'):
+ options['rtol'] = 1e3
+ del options['step']
+ if (method == 'trapz'):
+ options['lmm_type'] = 'ADAMS'
+ options['nonlinsolver'] = 'fixedpoint'
+ options['order'] = 2
+ else:
+ int_name = method
+
+ solver = ode(int_name, rhsfun, **options)
+ return solver.solve(t, y0)
+
+
+from . import _version
+__version__ = _version.get_versions()['version']
+
+
+
+[docs]
+class DaeBase:
+ """
+ The interface which DAE solvers must implement.
+
+ Parameters
+ ----------
+ Rfn : function
+ residual function
+ options : mapping
+ Additional options for initialization, solver dependent
+ """
+
+ integrator_classes = []
+
+ def __init__(self, Rfn, **options):
+ raise NotImplementedError('all DAE solvers must implement this')
+
+
+[docs]
+ def set_options(self, **options):
+ """
+ Set specific options for the solver.
+
+ Calling set_options a second time, normally resets the solver.
+ """
+ raise NotImplementedError('all DAE solvers must implement this')
+
+
+
+[docs]
+ def solve(self, tspan, y0, yp0):
+ """
+ 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 of initial values
+ yp0 : list/array
+ List 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
+ ========== ==========================================
+
+ """
+ raise NotImplementedError('all DAE solvers must implement this')
+
+
+
+[docs]
+ def init_step(self, t0, y0, yp0, y_ic0_retn = None, yp_ic0_retn = None):
+ """
+ Initializes the solver and allocates memory.
+
+ 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)
+ =========== ==========================================
+
+ """
+ raise NotImplementedError('all DAE solvers must implement this')
+
+
+
+[docs]
+ def step(self, t, y_retn=None, yp_retn=None):
+ """
+ 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)
+ =========== ==========================================
+
+ """
+ raise NotImplementedError('all DAE solvers must implement this')
+
+
+
+
+
+[docs]
+class OdeBase:
+ """
+ The interface which ODE solvers must implement.
+
+ Parameters
+ ----------
+ Rfn : function
+ A function which computes the required derivatives. The signature
+ should be ``func(t, y, y_dot, *args, **kwargs)``. Note that ``*args``
+ and ``**kwargs`` handling are solver dependent.
+
+ options : mapping
+ Additional options for initialization, solver dependent
+ """
+
+ integrator_classes = []
+
+ def __init__(self, Rfn, **options):
+ raise NotImplementedError('all ODE solvers must implement this')
+
+
+[docs]
+ def set_options(self, **options):
+ """
+ Set specific options for the solver.
+
+ Calling set_options a second time, normally resets the solver.
+ """
+ raise NotImplementedError('all ODE solvers must implement this')
+
+
+
+[docs]
+ def solve(self, tspan, y0):
+ """
+ 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 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
+ ========== ==========================================
+ """
+ raise NotImplementedError('all ODE solvers must implement this')
+
+
+
+[docs]
+ def init_step(self, t0, y0):
+ """
+ Initializes the solver and allocates memory.
+
+ 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 occurred)
+ ``t_out`` initial time
+ ========== ==========================================
+ """
+ raise NotImplementedError('all ODE solvers must implement this')
+
+
+
+[docs]
+ def step(self, t, y_retn=None):
+ """
+ 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 performed 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 occurred)
+ ``t_out`` time, where the solver stopped (when no error occurred, t_out == t)
+ ========== ==========================================
+ """
+ raise NotImplementedError('all ODE solvers must implement this')
+
+
+
+"""
+SUNDIALS wrapper
+"""
+
+import inspect
+from . import _version
+__version__ = _version.get_versions()['version']
+
+
+
+[docs]
+class CVODESolveException(Exception):
+ """Base class for exceptions raised by ``CVODE.validate_flags``."""
+ def __init__(self, soln):
+ self.soln = soln
+ self.args = (self._message.format(soln),)
+
+
+
+[docs]
+class CVODESolveFailed(CVODESolveException):
+ """``CVODE.solve`` failed to reach endpoint"""
+ _message = (
+ "Solver failed with flag {0.flag} and finished at {0.errors.t}"
+ "with values {0.errors.y}."
+ )
+
+
+
+[docs]
+class CVODESolveFoundRoot(CVODESolveException):
+ """``CVODE.solve`` found a root"""
+ _message = "Solver found a root at {0.roots.t[0]}."
+
+
+
+[docs]
+class CVODESolveReachedTSTOP(CVODESolveException):
+ """``CVODE.solve`` reached the endpoint specified by tstop."""
+ _message = "Solver reached tstop at {0.tstop.t[0]}."
+
+
+
+[docs]
+class IDASolveException(Exception):
+ """Base class for exceptions raised by ``IDA.validate_flags``."""
+ def __init__(self, soln):
+ self.soln = soln
+ self.args = (self._message.format(soln),)
+
+
+
+[docs]
+class IDASolveFailed(IDASolveException):
+ """``IDA.solve`` failed to reach endpoint"""
+ _message = (
+ "Solver failed with flag {0.flag} and finished at {0.errors.t}"
+ "with values {0.errors.y} and derivatives {0.errors.ydot}."
+ )
+
+
+
+[docs]
+class IDASolveFoundRoot(IDASolveException):
+ """``IDA.solve`` found a root"""
+ _message = "Solver found a root at {0.roots.t[0]}."
+
+
+
+[docs]
+class IDASolveReachedTSTOP(IDASolveException):
+ """``IDA.solve`` reached the endpoint specified by tstop."""
+ _message = "Solver reached tstop at {0.tstop.t[0]}."
+
+
+
+def _get_num_args(func):
+ """
+ Python 2/3 compatible method of getting number of args that `func` accepts
+ """
+ if hasattr(inspect, "getfullargspec"):
+ argspec = inspect.getfullargspec(func)
+ else:
+ argspec = inspect.getargspec(func)
+ arg_cnt = 0
+ for arg in argspec.args:
+ if arg not in ("self", "cls"):
+ arg_cnt += 1
+ return arg_cnt
+
Note
" + + "" + msg + "
"; + var parent = document.querySelector('div.body') + || document.querySelector('div.document') + || document.body; + parent.insertBefore(warning, parent.firstChild); + } + + +} + +function addVersionsMenu() { + // We assume that we can load versions.json from + // https://' + + '' + + _("Hide Search Matches") + + "
" + ) + ); + }, + + /** + * helper function to hide the search marks again + */ + hideSearchWords: () => { + document + .querySelectorAll("#searchbox .highlight-link") + .forEach((el) => el.remove()); + document + .querySelectorAll("span.highlighted") + .forEach((el) => el.classList.remove("highlighted")); + localStorage.removeItem("sphinx_highlight_terms") + }, + + initEscapeListener: () => { + // only install a listener if it is really needed + if (!DOCUMENTATION_OPTIONS.ENABLE_SEARCH_SHORTCUTS) return; + + document.addEventListener("keydown", (event) => { + // bail for input elements + if (BLACKLISTED_KEY_CONTROL_ELEMENTS.has(document.activeElement.tagName)) return; + // bail with special keys + if (event.shiftKey || event.altKey || event.ctrlKey || event.metaKey) return; + if (DOCUMENTATION_OPTIONS.ENABLE_SEARCH_SHORTCUTS && (event.key === "Escape")) { + SphinxHighlight.hideSearchWords(); + event.preventDefault(); + } + }); + }, +}; + +_ready(() => { + /* Do not call highlightSearchWords() when we are on the search page. + * It will highlight words from the *previous* search query. + */ + if (typeof Search === "undefined") SphinxHighlight.highlightSearchWords(); + SphinxHighlight.initEscapeListener(); +}); diff --git a/v3.0.0a3/_static/theme_overrides.css b/v3.0.0a3/_static/theme_overrides.css new file mode 100644 index 0000000..eb66180 --- /dev/null +++ b/v3.0.0a3/_static/theme_overrides.css @@ -0,0 +1,14 @@ +/* from https://rackerlabs.github.io/docs-rackspace/tools/rtd-tables.html */ +/* override table width restrictions */ +@media screen and (min-width: 767px) { + + .wy-table-responsive table td { + /* !important prevents the common CSS stylesheets from overriding + this as on RTD they are loaded after this stylesheet */ + white-space: normal !important; + } + + .wy-table-responsive { + overflow: visible !important; + } +} diff --git a/v3.0.0a3/api/compat.html b/v3.0.0a3/api/compat.html new file mode 100644 index 0000000..5c45a6b --- /dev/null +++ b/v3.0.0a3/api/compat.html @@ -0,0 +1,173 @@ + + + + + + + +scikits.odes is a scikit offering a cython wrapper around some extra ode/dae +solvers, so they can mature outside of scipy.
+CVODE
IDA
ddaspk <http://www.netlib.org/ode/ddaspk.f> (included)
lsodi <http://www.netlib.org/ode/lsodi.f> (included)
The interface which DAE solvers must implement.
+Rfn (function) – residual function
options (mapping) – Additional options for initialization, solver dependent
Initializes the solver and allocates memory.
+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.
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnumXXX) |
+
|
+Named tuple with entries t and y and ydot. y will correspond to y_retn value and ydot to yp_retn! |
+
|
+Named tuple with entries t and y and ydot |
+
|
+Named tuple with entries t and y and ydot |
+
|
+Named tuple with entries t and y and ydot |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+status of the computation (successful or error occurred) |
+
|
+time, where the solver stopped (when no error occurred, t_out == t) |
+
Set specific options for the solver.
+Calling set_options a second time, normally resets the solver.
+Runs the solver.
+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 of initial values
yp0 (list/array) – List of initial values of derivatives
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnumXXX) |
+
|
+Named tuple with entries array t and array y and array ydot. y will correspond to y_retn value and ydot to yp_retn! |
+
|
+Named tuple with entries t and y and ydot of error |
+
|
+Named tuple with entries array t and array y and array ydot |
+
|
+Named tuple with entries array t and array y and array ydot |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+indicating return status of the solver |
+
|
+numpy array of times at which the computations were successful |
+
|
+numpy array of values corresponding to times t (values of y[i, :] ~ t[i]) |
+
|
+numpy array of derivatives corresponding to times t (values of yp[i, :] ~ t[i]) |
+
|
+float or None - if recoverable error occurred (for example reached maximum number of allowed iterations), this is the time at which it happened |
+
|
+numpy array of values corresponding to time t_err |
+
|
+numpy array of derivatives corresponding to time t_err |
+
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.
+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.
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnumXXX) |
+
|
+Named tuple with entries t and y and ydot. y will correspond to y_retn value and ydot to yp_retn! |
+
|
+Named tuple with entries t and y and ydot |
+
|
+Named tuple with entries t and y and ydot |
+
|
+Named tuple with entries t and y and ydot |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+status of the computation (successful or error occurred) |
+
|
+time, where the solver stopped (when no error occurred, t_out == t) |
+
The interface which ODE solvers must implement.
+Rfn (function) – A function which computes the required derivatives. The signature
+should be func(t, y, y_dot, *args, **kwargs)
. Note that *args
+and **kwargs
handling are solver dependent.
options (mapping) – Additional options for initialization, solver dependent
Initializes the solver and allocates memory.
+t0 (number) – initial time
y0 (array) – initial condition for y (can be list or numpy array)
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnum) |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+boolean status of the computation (successful or error occurred) |
+
|
+initial time |
+
Set specific options for the solver.
+Calling set_options a second time, normally resets the solver.
+Runs the solver.
+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
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnum) |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+indicating return status of the solver |
+
|
+numpy array of times at which the computations were successful |
+
|
+numpy array of values corresponding to times t (values of y[i, :] ~ t[i]) |
+
|
+float or None - if recoverable error occurred (for example reached maximum number of allowed iterations), this is the time at which it happened |
+
|
+numpy array of values corresponding to time t_err |
+
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 performed towards time abs(t) +and results after this one time step are returned
+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.
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnum) |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+status of the computation (successful or error occurred) |
+
|
+time, where the solver stopped (when no error occurred, t_out == t) |
+
scikits-odes is a scikit offering a cython wrapper around some extra ode/dae +solvers, so they can mature outside of scipy.
+CVODE
IDA
ddaspk <http://www.netlib.org/ode/ddaspk.f> (included)
lsodi <http://www.netlib.org/ode/lsodi.f> (included)
User-friendly interface to various numerical integrators for solving a +system of first order ODEs with prescribed initial conditions:
+\(f(t,y)\) is the right hand side function and returns a vector of size +\(\mathrm{len}(y_0)\).
+A generic interface class to differential equation solvers.
+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)
+
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.
+Return additional information about the state of the integrator.
+A dictionary filled with internal data as exposed by the integrator.
See the get_info method of your chosen integrator for details.
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.
+t0 (number) – initial time
y0 (array) – initial condition for y (can be list or numpy array)
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnum) |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+boolean status of the computation (successful or error occured) |
+
|
+inititial time |
+
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.
+Add a stop time to the integrator past which he is not allowed to +integrate.
+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)
+Runs the solver.
+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
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnum) |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+indicating return status of the solver |
+
|
+numpy array of times at which the computations were successful |
+
|
+numpy array of values corresponding to times t (values of y[i, :] ~ t[i]) |
+
|
+float or None - if recoverable error occured (for example reached maximum number of allowed iterations), this is the time at which it happened |
+
|
+numpy array of values corresponding to time t_err |
+
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
+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.
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnum) |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+status of the computation (successful or error occured) |
+
|
+time, where the solver stopped (when no error occured, t_out == t) |
+
User-friendly interface to various numerical integrators for solving an +algebraic system of first order ODEs with prescribed initial conditions:
+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
+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.
+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)
+
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.
+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.
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnumXXX) |
+
|
+Named tuple with entries t and y and ydot. y will correspond to y_retn value and ydot to yp_retn! |
+
|
+Named tuple with entries t and y and ydot |
+
|
+Named tuple with entries t and y and ydot |
+
|
+Named tuple with entries t and y and ydot |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+status of the computation (successful or error occurred) |
+
|
+time, where the solver stopped (when no error occurred, t_out == t) |
+
Set specific options for the solver. +See the solver documentation for details.
+Calling set_options a second time, normally resets the solver.
+Runs the solver.
+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
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnumXXX) |
+
|
+Named tuple with entries array t and array y and array ydot. y will correspond to y_retn value and ydot to yp_retn! |
+
|
+Named tuple with entries t and y and ydot of error |
+
|
+Named tuple with entries array t and array y and array ydot |
+
|
+Named tuple with entries array t and array y and array ydot |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+indicating return status of the solver |
+
|
+numpy array of times at which the computations were successful |
+
|
+numpy array of values corresponding to times t (values of y[i, :] ~ t[i]) |
+
|
+numpy array of derivatives corresponding to times t (values of yp[i, :] ~ t[i]) |
+
|
+float or None - if recoverable error occurred (for example reached maximum number of allowed iterations), this is the time at which it happened |
+
|
+numpy array of values corresponding to time t_err |
+
|
+numpy array of derivatives corresponding to time t_err |
+
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.
+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.
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnumXXX) |
+
|
+Named tuple with entries t and y and ydot. y will correspond to y_retn value and ydot to yp_retn! |
+
|
+Named tuple with entries t and y and ydot |
+
|
+Named tuple with entries t and y and ydot |
+
|
+Named tuple with entries t and y and ydot |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+status of the computation (successful or error occurred) |
+
|
+time, where the solver stopped (when no error occurred, t_out == t) |
+
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)
+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”.
+Base class for exceptions raised by DOP.validate_flags.
+DOP.solve failed to reach endpoint
+Alias for field number 2
+Alias for field number 0
+Alias for field number 5
+Alias for field number 3
+Alias for field number 4
+Alias for field number 1
+Alias for field number 0
+Alias for field number 1
+Initializes the solver and allocates memory.
+time (t0 - initial)
array) (y0 - initial condition for y (can be list or numpy)
if old_api – not supported
if old_api False –
+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 specific options for the solver.
+Calling set_options a second time, normally resets the solver.
+Runs the solver.
+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)
if old_api – Not supported
if old_api False –
+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
+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.
+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 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 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.
if old_api – not supported
if old_api False –
+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
+SUNDIALS wrapper
+Base class for exceptions raised by CVODE.validate_flags
.
CVODE.solve
failed to reach endpoint
CVODE.solve
reached the endpoint specified by tstop.
Base class for exceptions raised by IDA.validate_flags
.
Simple wrapper for functions called when ROOT or TSTOP are returned.
+Prototype for jacobian function.
+Note that evaluate must return a integer, 0 for success, positive for +recoverable failure, negative for unrecoverable failure (as per CVODE +documentation).
+Returns the Jacobi matrix of the right hand side function, as:
+++d(rhs)/d y
+
(for dense the full matrix, for band only bands). Result has to be +stored in the variable J, which is preallocated to the corresponding +size.
+This is a generic class, you should subclass is for the problem specific +purposes.
+Prototype for jacobian times setup function.
+Note that evaluate must return a integer, 0 for success, non-zero for error +(as per CVODE documentation), with >0 a recoverable error (step is retried).
+This function calculates the product of the Jacobian with a given vector v. +Use the userdata object to expose Jacobian related data to the solve function.
+This is a generic class, you should subclass it for the problem specific +purposes.
+Prototype for jacobian times vector function.
+Note that evaluate must return a integer, 0 for success, non-zero for error +(as per CVODE documentation).
+This function calculates the product of the Jacobian with a given vector v. +Use the userdata object to expose Jacobian related data to the solve function.
+This is a generic class, you should subclass it for the problem specific +purposes.
+Prototype for preconditioning setup function.
+Note that evaluate must return a integer, 0 for success, positive for +recoverable failure, negative for unrecoverable failure (as per CVODE +documentation).
+This function preprocesses and/or evaluates Jacobian-related data +needed by the preconditioner. Use the userdata object to expose +the preprocessed data to the solve function.
+This is a generic class, you should subclass it for the problem specific +purposes.
+Prototype for precondititioning solution function.
+Note that evaluate must return a integer, 0 for success, positive for +recoverable failure, negative for unrecoverable failure (as per CVODE +documentation).
+This function solves the preconditioned system P*z = r, where P may be +either a left or right preconditioner matrix. Here P should approximate +(at least crudely) the Newton matrix M = I − gamma*J, where J is the +Jacobian of the system. If preconditioning is done on both sides, +the product of the two preconditioner matrices should approximate M.
+This is a generic class, you should subclass it for the problem specific +purposes.
+Prototype for rhs function.
+Note that evaluate must return a integer, 0 for success, positive for +recoverable failure, negative for unrecoverable failure (as per CVODE +documentation).
+Prototype for root function.
+Note that evaluate must return a integer, 0 for success, non-zero if error +(as per CVODE documentation).
+Returns the Jacobi matrix (for dense the full matrix, for band only +bands. Result has to be stored in the variable J, which is preallocated +to the corresponding size.
+Set some jacobian equations as a JacRhsFunction executable class.
+Set some CV_JacTimesSetupFn executable class.
+Set some CV_JacTimesVecFn executable class.
+set a precondititioning setup method as a CV_PrecSetupFunction +executable class
+set a precondititioning solve method as a CV_PrecSolveFunction +executable class
+set some rhs equations as a RhsFunction executable class
+‘int’
+with_userdata
+set root-ing condition(equations) as a RootFunction executable class
+Simple wrapper for functions called when ROOT or TSTOP are returned.
+Prototype for jacobian function.
+Note that evaluate must return a integer, 0 for success, positive for +recoverable failure, negative for unrecoverable failure (as per IDA +documentation).
+Returns the Jacobi matrix of the residual function, as:
+++d(res)/d y + cj d(res)/d ydot
+
(for dense the full matrix, for band only bands). Result has to be +stored in the variable J, which is preallocated to the corresponding +size.
+This is a generic class, you should subclass is for the problem specific +purposes.”
+Prototype for jacobian times setup function.
+Note that evaluate must return a integer, 0 for success, non-zero for error +(as per CVODE documentation), with >0 a recoverable error (step is retried).
+This function calculates the product of the Jacobian with a given vector v. +Use the userdata object to expose Jacobian related data to the solve function.
+This is a generic class, you should subclass it for the problem specific +purposes.
+Prototype for jacobian times vector function.
+Note that evaluate must return a integer, 0 for success, non-zero for error +(as per IDA documentation).
+This function calculates the product of the Jacobian with a given vector v. +Use the userdata object to expose Jacobian related data to the solve function.
+This is a generic class, you should subclass it for the problem specific +purposes.
+Prototype for preconditioning setup function.
+Note that evaluate must return a integer, 0 for success, positive for +recoverable failure, negative for unrecoverable failure (as per CVODE +documentation).
+This function preprocesses and/or evaluates Jacobian-related data +needed by the preconditioner. Use the userdata object to expose +the preprocessed data to the solve function.
+This is a generic class, you should subclass it for the problem specific +purposes.
+Prototype for precondititioning solution function.
+Note that evaluate must return a integer, 0 for success, positive for +recoverable failure, negative for unrecoverable failure (as per CVODE +documentation).
+This function solves the preconditioned system P*z = r, where P may be +either a left or right preconditioner matrix. Here P should approximate +(at least crudely) the Newton matrix M = I − gamma*J, where J is the +Jacobian of the system. If preconditioning is done on both sides, +the product of the two preconditioner matrices should approximate M.
+This is a generic class, you should subclass it for the problem specific +purposes.
+Prototype for rhs function.
+Note that evaluate must return a integer, 0 for success, positive for +recoverable failure, negative for unrecoverable failure (as per IDA +documentation).
+Prototype for root function.
+Note that evaluate must return a integer, 0 for success, non-zero for error +(as per IDA documentation).
+Returns the Jacobi matrix (for dense the full matrix, for band only +bands. Result has to be stored in the variable J, which is preallocated +to the corresponding size.
+Set some jacobian equations as a JacResFunction executable class.
+Set some IDA_JacTimesSetupFn executable class.
+Set some IDA_JacTimesVecFn executable class.
+set a precondititioning setup method as a IDA_PrecSetupFunction +executable class
+set a precondititioning solve method as a IDA_PrecSolveFunction +executable class
+set some residual equations as a ResFunction executable class
+set root-ing condition(equations) as a RootFunction executable class
+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.
+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)
+
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.
+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.
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnumXXX) |
+
|
+Named tuple with entries t and y and ydot. y will correspond to y_retn value and ydot to yp_retn! |
+
|
+Named tuple with entries t and y and ydot |
+
|
+Named tuple with entries t and y and ydot |
+
|
+Named tuple with entries t and y and ydot |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+status of the computation (successful or error occurred) |
+
|
+time, where the solver stopped (when no error occurred, t_out == t) |
+
Set specific options for the solver. +See the solver documentation for details.
+Calling set_options a second time, normally resets the solver.
+Runs the solver.
+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
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnumXXX) |
+
|
+Named tuple with entries array t and array y and array ydot. y will correspond to y_retn value and ydot to yp_retn! |
+
|
+Named tuple with entries t and y and ydot of error |
+
|
+Named tuple with entries array t and array y and array ydot |
+
|
+Named tuple with entries array t and array y and array ydot |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+indicating return status of the solver |
+
|
+numpy array of times at which the computations were successful |
+
|
+numpy array of values corresponding to times t (values of y[i, :] ~ t[i]) |
+
|
+numpy array of derivatives corresponding to times t (values of yp[i, :] ~ t[i]) |
+
|
+float or None - if recoverable error occurred (for example reached maximum number of allowed iterations), this is the time at which it happened |
+
|
+numpy array of values corresponding to time t_err |
+
|
+numpy array of derivatives corresponding to time t_err |
+
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.
+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.
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnumXXX) |
+
|
+Named tuple with entries t and y and ydot. y will correspond to y_retn value and ydot to yp_retn! |
+
|
+Named tuple with entries t and y and ydot |
+
|
+Named tuple with entries t and y and ydot |
+
|
+Named tuple with entries t and y and ydot |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+status of the computation (successful or error occurred) |
+
|
+time, where the solver stopped (when no error occurred, t_out == t) |
+
+ | + |
+ |
+ |
|
+
+ |
|
+ + |
+ |
+ | + |
+ | + |
+ |
+ |
The ODES scikit provides access to Ordinary Differential Equation (ODE) solvers and Differential Algebraic Equation (DAE) solvers.
+This is the generated API documentation for version 3.0.0a3, see https://scikits-odes.readthedocs.io/en/latest/ for the main documentation and user guide.
+For other versions of the API documentation, see https://bmcage.github.io/odes.
+Contents:
+dae
+CV_ContinuationFunction
CV_JacRhsFunction
CV_JacTimesSetupFunction
CV_JacTimesVecFunction
CV_PrecSetupFunction
CV_PrecSolveFunction
CV_RhsFunction
CV_RootFunction
CV_WrapJacRhsFunction
CV_WrapJacTimesSetupFunction
CV_WrapJacTimesVecFunction
CV_WrapPrecSetupFunction
CV_WrapPrecSolveFunction
CV_WrapRhsFunction
CV_WrapRootFunction
StatusEnum
no_continue_fn()
IDA_ContinuationFunction
IDA_JacRhsFunction
IDA_JacTimesSetupFunction
IDA_JacTimesVecFunction
IDA_PrecSetupFunction
IDA_PrecSolveFunction
IDA_RhsFunction
IDA_RootFunction
IDA_WrapJacRhsFunction
IDA_WrapJacTimesSetupFunction
IDA_WrapJacTimesVecFunction
IDA_WrapPrecSetupFunction
IDA_WrapPrecSolveFunction
IDA_WrapRhsFunction
IDA_WrapRootFunction
StatusEnumIDA
no_continue_fn()
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.
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 non-stiff 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
+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 Adams-Moulton method (linear +multistep method), good for non-stiff ODE. Functional +iterations are used to solve the nonlinear system. |
+
rk5 |
+This uses the ‘dopri5’ solver, which is a variable step +Runge-Kutta method of order (4)5 (use for non-stiff ODE) |
+
rk8 |
+This uses the ‘dop853’ solver, which is a variable step +Runge-Kutta 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
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 |
+
named tuple
+See also
+scikits.odes.ode.ode
a more object-oriented integrator
+scikits.odes.dae.dae
a solver for differential-algebraic 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:
+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:
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()
+
A generic interface class to differential equation solvers.
+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)
+
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.
+Return additional information about the state of the integrator.
+A dictionary filled with internal data as exposed by the integrator.
See the get_info method of your chosen integrator for details.
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.
+t0 (number) – initial time
y0 (array) – initial condition for y (can be list or numpy array)
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnum) |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+boolean status of the computation (successful or error occured) |
+
|
+inititial time |
+
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.
+Add a stop time to the integrator past which he is not allowed to +integrate.
+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)
+Runs the solver.
+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
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnum) |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+indicating return status of the solver |
+
|
+numpy array of times at which the computations were successful |
+
|
+numpy array of values corresponding to times t (values of y[i, :] ~ t[i]) |
+
|
+float or None - if recoverable error occured (for example reached maximum number of allowed iterations), this is the time at which it happened |
+
|
+numpy array of values corresponding to time t_err |
+
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
+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.
old_api is False (namedtuple) – namedtuple with the following attributes
+Field |
+Meaning |
+
---|---|
|
+An integer flag (StatusEnum) |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+Named tuple with entries t and y |
+
|
+String with message in case of an error |
+
old_api is True (tuple) – tuple with the following elements in order
+Field |
+Meaning |
+
---|---|
|
+status of the computation (successful or error occured) |
+
|
+time, where the solver stopped (when no error occured, t_out == t) |
+
+ s | ||
+ |
+ scikits | + |
+ |
+ scikits.odes | + |
+ |
+ scikits.odes.dae | + |
+ |
+ scikits.odes.ddaspkint | + |
+ |
+ scikits.odes.dopri5 | + |
+ |
+ scikits.odes.lsodiint | + |
+ |
+ scikits.odes.ode | + |
+ |
+ scikits.odes.sundials | + |
+ |
+ scikits.odes.sundials.cvode | + |
+ |
+ scikits.odes.sundials.ida | + |
+ |
+ scikits_odes | + |
+ |
+ scikits_odes.dae | + |
+ |
+ scikits_odes.dopri5 | + |
+ |
+ scikits_odes.ode | + |
+ |
+ scikits_odes_core | + |
+ |
+ scikits_odes_sundials | + |
+ |
+ scikits_odes_sundials.cvode | + |
+ |
+ scikits_odes_sundials.ida | + |
+ Searching for multiple words only shows matches that contain + all words. +
+ + + + + + + + +Sundials can be built with different precisions (currently ‘single’ which maps +to C ‘float’; ‘double’ (the default) which maps to C ‘double’; and +‘extended’ which maps to C ‘long double’), and scikits-odes-sundials +supports using whichever precision Sundials is built with. To take advantage of +this, you must first have sundials built and installed with the desired +precision setting.
+Once you have done this, build scikits-odes-sundials against this particular
+version of sundials (see the main documentation for how to do this).
+scikits-odes-sundials will automatically detect the precision, and store this in
+a variable called
+DTYPE
. DTYPE
should be accessed from
+scikits_odes_sundials
(other modules may have DTYPE
+defined, but scikits_odes_sundials
should be preferred). Additionally
+scikits_odes_sundials.precision
contains the precision setting found
+by scikits.odes.
To use DTYPE
, treat it as a numpy dtype; use it whenever you need to
+create an array:
np.array([1,2,3], dtype=DTYPE)
+
or when using scalars:
+DTYPE(0.1) * DTYPE(0.1)
+