pyodesys package

Straightforward numerical integration of ODE systems from SymPy.

The pyodesys package for enables intuitive solving of systems of first order differential equations numerically. The system may be symbolically defined.

pyodesys ties computer algebra systems (like SymPy and symengine), and numerical solvers (such as ODEPACK in SciPy, CVode in sundials, GSL or odeint) together.

Submodules

pyodesys.core module

Core functionality from OdeSys.

Note that it is possible to use new custom ODE integrators with pyodesys by providing a module with two functions named integrate_adaptive and integrate_predefined. See the pyodesys.integrators module for examples.

class pyodesys.core.OdeSys(f, jac=None, dfdx=None, roots=None, nroots=None, band=None, names=None, pre_processors=None, post_processors=None)[source]

Bases: object

Object representing an ODE system.

OdeSys provides unified interface to:

  • scipy.integarte.ode
  • pygslodeiv2
  • pyodeint
  • pycvodes

The numerical integration can be performed either in an adaptive() or predefined() mode. Where locations to report the solution is chosen by the stepper or the user respectively. For convenience in user code one may use integrate() which automatically chooses between the two based on the length of xout provided by the user.

Parameters:

f: callback

first derivatives of dependent variables (y) with respect to dependent variable (x). Signature is any of:

  • rhs(x, y[:]) –> f[:]
  • rhs(x, y[:], p[:]) –> f[:]
  • rhs(x, y[:], p[:], backend=math) –> f[:]

jac: callback

Jacobian matrix (dfdy). Required for implicit methods.

dfdx: callback

Signature dfdx(x, y[:], p[:]) -> out[:] (used by e.g. GSL)

band: tuple of 2 integers or None (default: None)

If jacobian is banded: number of sub- and super-diagonals

names: iterable of strings (default: None)

names of variables, e.g. used for plotting

pre_processors: iterable of callables (optional)

signature: f(x1[:], y1[:], params1[:]) -> x2[:], y2[:], params2[:]. When modifying: insert at beginning.

post_processors: iterable of callables (optional)

signature: f(x2[:], y2[:, :], params2[:]) -> x1[:], y1[:, :], params1[:] When modifying: insert at end.

Notes

banded jacobians are supported by “scipy” and “cvode” integrators

Examples

>>> odesys = OdeSys(lambda x, y, p: p[0]*x + p[1]*y[0]*y[0])
>>> yout, info = odesys.predefined([1], [0, .2, .5], [2, 1])
>>> print(info['success'])
True

Attributes

f_cb (callback) for evaluating the vector of derivatives
j_cb (callback) for evaluating the Jacobian matrix of f
roots_cb (callback)
nroots (int)
names (iterable of strings)
internal_xout (1D array of floats) internal values of dependent variable before post-processing
internal_yout (2D (or higher) array of floats) internal values of dependent variable before post-processing
internal_params (1D array of floats) internal parameter values before post-processing

Methods

adaptive(y0, x0, xend[, params]) Integrate with integrator chosen output.
integrate(x, y0[, params]) Integrate the system of ODE’s.
plot_phase_plane([indices]) Plots a phase portrait from last integration.
plot_result(**kwargs) Plots the integrated dependent variables from last integration.
post_process(xout, yout, params) Transforms internal values to output, used inernally.
pre_process(xout, y0[, params]) Transforms input to internal values, used inernally.
predefined(y0, xout[, params]) Integrate with user chosen output.
stiffness([xyp, eigenvals_cb]) Running stiffness ratio from last integration.
adaptive(y0, x0, xend, params=(), **kwargs)[source]

Integrate with integrator chosen output.

Parameters:

integrator: str

y0: array_like

x0: float

initial value of the independent variable

xend: float

final value of the independent variable

params: array_like

**kwargs:

Returns:

Same as integrate()

integrate(x, y0, params=(), **kwargs)[source]

Integrate the system of ODE’s.

Parameters:

x: array_like or pair (start and final time) or float

if float:

make it a pair: (0, x)

if pair or length-2 array:

initial and final value of the independent variable

if array_like:

values of independent variable report at

y0: array_like

Initial values at x[0] for the dependent variables.

params: array_like (default: tuple())

Value of parameters passed to user-supplied callbacks.

integrator: str or None

Name of integrator, one of:
  • ‘scipy’: _integrate_scipy()
  • ‘gsl’: _integrate_gsl()
  • ‘odeint’: _integrate_odeint()
  • ‘cvode’: _integrate_cvode()

See respective method for more information. If None: os.environ.get('PYODESYS_INTEGRATOR', 'scipy')

atol: float

Absolute tolerance

rtol: float

Relative tolerance

with_jacobian: bool or None (default)

Whether to use the jacobian. When None the choice is done automatically (only used when required). This matters when jacobian is derived at runtime (high computational cost).

force_predefined: bool (default: False)

override behaviour of len(x) == 2 => adaptive()

**kwargs:

Additional keyword arguments for _integrate_$(integrator).

Returns:

Length 3 tuple: (x, yout, info)

x: array of values of the independent variable

yout: array of the dependent variable(s) for the different values of x

info: dict (‘nfev’ is guaranteed to be a key)

plot_phase_plane(indices=None, **kwargs)[source]

Plots a phase portrait from last integration.

See pyodesys.plotting.plot_phase_plane()

plot_result(**kwargs)[source]

Plots the integrated dependent variables from last integration.

See pyodesys.plotting.plot_result()

post_process(xout, yout, params)[source]

Transforms internal values to output, used inernally.

pre_process(xout, y0, params=())[source]

Transforms input to internal values, used inernally.

predefined(y0, xout, params=(), **kwargs)[source]

Integrate with user chosen output.

Parameters:

integrator: str

y0: array_like

xout: array_like

params: array_like

**kwargs:

Returns:

Length 2 tuple: (yout, info)

stiffness(xyp=None, eigenvals_cb=None)[source]

Running stiffness ratio from last integration.

Calculate sittness ratio, i.e. the ratio between the largest and smallest absolute eigenvalue of the jacobian matrix. The user may supply their own routine for calculating the eigenvalues, or they will be calculated from the SVD (singular value decomposition). Note that calculating the SVD for any but the smallest Jacobians may prove to be prohibitively expensive.

Parameters:

xyp: length 3 tuple (default: None)

internal_xout, internal_yout, internal_params, taken from last integration if not specified.

eigenvals_cb: callback (optional)

signature (x, y, p) (internal variables), when not provided an internal routine will use self.j_cb and scipy.linalg.svd.

pyodesys.integrators module

This module is for demonstration purposes only and the integrators here are not meant for production use. Consider them provisional, i.e., API here may break without prior deprecation.

class pyodesys.integrators.RK4_example_integartor[source]

This is an example of how to implement a custom integrator. It uses fixed step size and is usually not useful for real problems.

Methods

integrate_adaptive(rhs, jac, y0, x0, xend, ...)
integrate_predefined(rhs, jac, y0, xout, ...)
static integrate_adaptive(rhs, jac, y0, x0, xend, dx0, **kwargs)[source]
static integrate_predefined(rhs, jac, y0, xout, **kwargs)[source]
with_jacobian = False

pyodesys.plotting module

pyodesys.plotting.plot_phase_plane(x, y, params=(), indices=None, post_processors=(), plot=None, names=None, **kwargs)[source]

Plot the phase portrait of two dependent variables

Parameters:

x: array_like

Values of the independent variable

y: array_like

Values of the dependent variables

params: array_like

parameters

indices: pair of integers (default: None)

what dependent variable to plot for (None => (0, 1))

post_processors: iterable of callbles

see :class:OdeSystem

plot: callable (default: None)

Uses matplotlib.pyplot.plot if None

names: iterable of strings

labels for x and y axis

**kwargs:

keyword arguemtns passed to plot()

pyodesys.plotting.plot_result(x, y, params=(), indices=None, plot=None, plot_kwargs_cb=None, ls=('-', '--', ':', '-.'), c=('k', 'r', 'g', 'b', 'c', 'm', 'y'), m=('o', 'v', '8', 's', 'p', 'x', '+', 'd', 's'), m_lim=-1, lines=None, interpolate=None, interp_from_deriv=None, names=None, post_processors=(), xlabel=None, ylabel=None, xscale=None, yscale=None, latex_attr='latex_name')[source]

Plot the depepndent variables vs. the independent variable

Parameters:

x : array_like

Values of the independent variable.

y : array_like

Values of the independent variable. This must hold y.shape[0] == len(x), plot_results will draw y.shape[1] lines. If interpolate != None y is expected two be three dimensional, otherwise two dimensional.

params : array_like

Parameters used.

indices : iterable of integers

What indices to plot (default: None => all).

plot : callback (default: None)

If None, use matplotlib.pyplot.plot.

plot_kwargs_cb : callback(int) -> dict

Keyword arguments for plot for each index (0:len(y)-1).

ls : iterable

Linestyles to cycle through (only used if plot and plot_kwargs_cb are both None).

c : iterable

Colors to cycle through (only used if plot and plot_kwargs_cb are both None).

m : iterable

Markers to cycle through (only used if plot and plot_kwargs_cb are both None and m_lim > 0).

m_lim : int (default: -1)

Upper limit (exclusive, number of points) for using markers instead of lines.

lines : None

default: draw between markers unless we are interpolating as well.

interpolate : bool or int (default: None)

Density-multiplier for grid of independent variable when interpolating if True => 20. negative integer signifies log-spaced grid.

interp_from_deriv : callback (default: None)

When None: scipy.interpolate.BPoly.from_derivatives

post_processors : iterable of callback (default: tuple())

pyodesys.symbolic module

class pyodesys.symbolic.PartiallySolvedSystem(original_system, analytic_factory, **kwargs)[source]

Bases: pyodesys.symbolic.SymbolicSys

Use analytic expressions for some dependent variables

Parameters:

original_system: SymbolicSys

analytic_factory: callable

signature: solved(x0, y0, p0) -> dict, where dict maps independent variables as analytic expressions in remaining variables

Examples

>>> import sympy as sp
>>> odesys = SymbolicSys.from_callback(
...     lambda x, y, p: [
...         -p[0]*y[0],
...         p[0]*y[0] - p[1]*y[1]
...     ], 2, 2)
>>> dep0 = odesys.dep[0]
>>> partsys = PartiallySolvedSystem(odesys, lambda x0, y0, p0: {
...         dep0: y0[0]*sp.exp(-p0[0]*(odesys.indep-x0))
...     })
>>> print(partsys.exprs)  
(_Dummy_29*p_0*exp(-p_0*(-_Dummy_28 + x)) - p_1*y_1,)
>>> y0, k = [3, 2], [3.5, 2.5]
>>> xout, yout, info = partsys.integrate([0, 1], y0, k, integrator='scipy')
>>> info['success'], yout.shape[1]
(True, 2)

Attributes

ny Number of dependent variables in the system.

Methods

analytic_stiffness([xyp]) Running stiffness ratio from last integration.
from_callback(cb, ny[, nparams]) Create an instance from a callback.
from_other(ori, **kwargs)
get_dfdx() Calculates 2nd derivatives of self.exprs
get_dfdx_callback() Generate a callback for evaluating derivative of self.exprs
get_f_ty_callback() Generates a callback for evaluating self.exprs.
get_j_ty_callback() Generates a callback for evaluating the jacobian.
get_jac() Derives the jacobian from self.exprs and self.dep.
get_roots_callback() Generate a callback for evaluating self.roots
class pyodesys.symbolic.ScaledSys(dep_exprs, indep=None, dep_scaling=1, indep_scaling=1, params=(), **kwargs)[source]

Bases: pyodesys.symbolic.TransformedSys

Transformed system where the variables have been scaled linearly.

Parameters:

dep_exprs: iterable of (symbol, expression)-pairs

indep: Symbol

dep_scaling: number (>0) or iterable of numbers

scaling of the dependent variables (default: 1)

indep_scaling: number (>0)

scaling of the independent variable (default: 1)

params:

**kwargs:

keyword arguments passed onto TransformedSys

Examples

>>> from sympy import Symbol
>>> x = Symbol('x')
>>> scaled = ScaledSys([(x, x*x)], dep_scaling=1000)
>>> scaled.exprs
(x**2/1000,)

Attributes

ny Number of dependent variables in the system.

Methods

analytic_stiffness([xyp]) Running stiffness ratio from last integration.
from_callback(cb, ny[, nparams, ...]) Create an instance from a callback.
from_other(ori, **kwargs)
get_dfdx() Calculates 2nd derivatives of self.exprs
get_dfdx_callback() Generate a callback for evaluating derivative of self.exprs
get_f_ty_callback() Generates a callback for evaluating self.exprs.
get_j_ty_callback() Generates a callback for evaluating the jacobian.
get_jac() Derives the jacobian from self.exprs and self.dep.
get_roots_callback() Generate a callback for evaluating self.roots
classmethod from_callback(cb, ny, nparams=0, dep_scaling=1, indep_scaling=1, **kwargs)[source]

Create an instance from a callback.

Analogous to SymbolicSys.from_callback().

Parameters:

cb: callable

Signature rhs(x, y[:], p[:]) -> f[:]

ny: int

length of y

nparams: int

length of p

dep_scaling: number (>0) or iterable of numbers

scaling of the dependent variables (default: 1)

indep_scaling: number (>0)

scaling of the independent variable (default: 1)

**kwargs:

keyword arguments passed onto ScaledSys

Examples

>>> def f(x, y, p):
...     return [p[0]*y[0]**2]
>>> odesys = ScaledSys.from_callback(f, 1, 1, dep_scaling=10)
>>> odesys.exprs
(p_0*y_0**2/10,)
class pyodesys.symbolic.SymbolicSys(dep_exprs, indep=None, params=(), jac=True, dfdx=True, roots=None, backend=None, **kwargs)[source]

Bases: pyodesys.core.OdeSys

ODE System from symbolic expressions

Creates a OdeSys instance from symbolic expressions. Jacboian and second derivatives are derived when needed.

Parameters:

dep_exprs : iterable of (symbol, expression)-pairs

indep : Symbol

Independent variable (default: None => autonomous system).

params : iterable of symbols

Problem parameters.

jac : ImmutableMatrix or bool (default: True)

if True:

calculate jacobian from exprs

if False:

do not compute jacobian (use explicit steppers)

if instance of ImmutableMatrix:

user provided expressions for the jacobian

roots : iterable of expressions

Equations to look for root’s for during integration (currently available through cvode).

backend : str or sym.Backend

See documentation of sym.Backend.

**kwargs:

See OdeSys

Attributes

ny Number of dependent variables in the system.
dep (iterable of symbols) Dependent variables.
indep (Symbol or None) Independent variable (None indicates autonomous system).
params (iterable of symbols) Problem parameters.
roots (iterable of expressions or None) Roots to report for during integration.
be (module) Symbolic backend.

Methods

analytic_stiffness([xyp]) Running stiffness ratio from last integration.
from_callback(cb, ny[, nparams]) Create an instance from a callback.
from_other(ori, **kwargs)
get_dfdx() Calculates 2nd derivatives of self.exprs
get_dfdx_callback() Generate a callback for evaluating derivative of self.exprs
get_f_ty_callback() Generates a callback for evaluating self.exprs.
get_j_ty_callback() Generates a callback for evaluating the jacobian.
get_jac() Derives the jacobian from self.exprs and self.dep.
get_roots_callback() Generate a callback for evaluating self.roots
analytic_stiffness(xyp=None)[source]

Running stiffness ratio from last integration.

Calculate sittness ratio, i.e. the ratio between the largest and smallest absolute eigenvalue of the (analytic) jacobian matrix.

See OdeSys.stiffness() for more info.

classmethod from_callback(cb, ny, nparams=0, **kwargs)[source]

Create an instance from a callback.

Parameters:

cb : callbable

Signature rhs(x, y[:], p[:]) -> f[:]

ny : int

length of y

nparams : int (default: 0)

length of p

**kwargs :

keyword arguments passed onto SymbolicSys

Returns:

An instance of SymbolicSys.

classmethod from_other(ori, **kwargs)[source]
get_dfdx()[source]

Calculates 2nd derivatives of self.exprs

get_dfdx_callback()[source]

Generate a callback for evaluating derivative of self.exprs

get_f_ty_callback()[source]

Generates a callback for evaluating self.exprs.

get_j_ty_callback()[source]

Generates a callback for evaluating the jacobian.

get_jac()[source]

Derives the jacobian from self.exprs and self.dep.

get_roots_callback()[source]

Generate a callback for evaluating self.roots

ny

Number of dependent variables in the system.

class pyodesys.symbolic.TransformedSys(dep_exprs, indep=None, dep_transf=None, indep_transf=None, params=(), exprs_process_cb=None, **kwargs)[source]

Bases: pyodesys.symbolic.SymbolicSys

SymbolicSys with abstracted variable transformations.

Parameters:

dep_exprs : iterable of pairs

indep : Symbol

dep_transf : iterable of (expression, expression) pairs

pairs of (forward, backward) transformations for the dependents variables

indep_transf : pair of expressions

forward and backward transformation of the independent variable

params :

exprs_process_cb : callbable

signatrue f(exprs) -> exprs post processing of the expressions for the derivatives of the dependent variables after transformation have been applied.

**kwargs :

keyword arguments passed onto SymbolicSys

Attributes

ny Number of dependent variables in the system.

Methods

analytic_stiffness([xyp]) Running stiffness ratio from last integration.
from_callback(cb, ny[, nparams, ...]) Create an instance from a callback.
from_other(ori, **kwargs)
get_dfdx() Calculates 2nd derivatives of self.exprs
get_dfdx_callback() Generate a callback for evaluating derivative of self.exprs
get_f_ty_callback() Generates a callback for evaluating self.exprs.
get_j_ty_callback() Generates a callback for evaluating the jacobian.
get_jac() Derives the jacobian from self.exprs and self.dep.
get_roots_callback() Generate a callback for evaluating self.roots
classmethod from_callback(cb, ny, nparams=0, dep_transf_cbs=None, indep_transf_cbs=None, **kwargs)[source]

Create an instance from a callback.

Analogous to SymbolicSys.from_callback().

Parameters:

cb : callable

Signature rhs(x, y[:], p[:]) -> f[:]

ny : int

length of y

nparams : int

length of p

dep_transf_cbs : iterable of pairs callables

callables should have the signature f(yi) -> expression in yi

indep_transf_cbs : pair of callbacks

callables should have the signature f(x) -> expression in x

**kwargs :

keyword arguments passed onto TransformedSys

pyodesys.symbolic.symmetricsys(dep_tr=None, indep_tr=None, **kwargs)[source]

A factory function for creating symmetrically transformed systems.

Parameters:

dep_tr : pair of callables (default: None)

Forward and backward transformation callbacks to be applied to the dependent variables.

indep_tr : pair of callables (default: None)

Forward and backward transformation to be applied to the independent variable.

**kwargs :

Default keyword arguments for the TransformedSys subclass.

Returns:

Subclass of TransformedSys.

Examples

>>> import sympy
>>> logexp = (sympy.log, sympy.exp)
>>> LogLogSys = symmetricsys(
...     logexp, logexp, exprs_process_cb=lambda exprs: [
...         sympy.powsimp(expr.expand(), force=True) for expr in exprs])

pyodesys.util module

pyodesys.util.check_transforms(fw, bw, symbs)[source]

Verify validity of a pair of forward and backward transformations

Parameters:

fw: expression

forward transformation

bw: expression

backward transformation

symbs: iterable of symbols

the variables that are transformed

pyodesys.util.stack_1d_on_left(x, y)[source]

Stack a 1D array on the left side of a 2D array

Parameters:

x: 1D array

y: 2D array

Requirement: shape[0] == x.size

pyodesys.util.transform_exprs_dep(fw, bw, dep_exprs, check=True)[source]

Transform y[:] in dydx

Parameters:

fw: expression

forward transformation

bw: expression

backward transformation

dep_exprs: iterable of (symbol, expression) pairs

pairs of (dependent variable, derivative expressions), i.e. (y, dydx) pairs

check: bool (default: True)

whether to verification of the analytic correctness should be performed

Returns:

List of transformed expressions for dydx

pyodesys.util.transform_exprs_indep(fw, bw, dep_exprs, indep, check=True)[source]

Transform x in dydx

Parameters:

fw: expression

forward transformation

bw: expression

backward transformation

dep_exprs: iterable of (symbol, expression) pairs

pairs of (dependent variable, derivative expressions)

check: bool (default: True)

whether to verification of the analytic correctness should be performed

Returns:

List of transformed expressions for dydx