scijit.integrate.odeint

scijit.integrate.odeint(func, y0, t, args=(), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=None, atol=None, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0, tfirst=False, success_out=False)

Integrate a system of ODEs with LSODA over a time grid.

LSODA switches between Adams (non-stiff) and BDF (stiff) methods automatically. The result is time-major: row i is the state at t[i].

Parameters:

func@njit function

The right-hand side, f(y, t, *args), RETURNING a new 1-D float64 array of neq derivatives. One parameter per entry of args. Under tfirst=True the order is f(t, y, *args). A plain @njit function.

y0float, sequence or ndarray

Initial state at t[0]. A scalar becomes length 1. Rank >= 2 raises ValueError.

tsequence or ndarray

Times to report at, monotonic; repeated values are allowed. t[0] is the initial time. Rank >= 2 raises ValueError. An empty t returns a (0, neq) array.

argstuple, optional

Extra parameters, passed to the callback as separate arguments. An entry may be a float, an int, a bool or an array of any rank, and arrives with its own type and shape. Anything other than a tuple raises ODEpackError("Extra arguments must be in a tuple."). A callback whose arity does not fit raises TypeError naming both counts.

Dfun@njit function or None, optional

The Jacobian d f / d y, j(y, t, *args), or j(t, y, *args) under tfirst=True. None (default) leaves LSODA to build one by finite differences, at neq extra right-hand-side evaluations per rebuild. A plain @njit function, as func is.

Shape (neq, neq) with jac[i, j] = d f_i / d y_j, or, with ml/mu set, the packed (ml + mu + 1, neq) form with jac[mu + i - j, j] = d f_i / d y_j. col_deriv transposes both. A wrong shape raises RuntimeError.

Measured on Robertson, 3 states, rtol=1e-8 atol=1e-10: nfe 502 without and 424 with, nje 26 either way, values agreeing to 1.369777e-11. The 78 evaluations removed are exactly the 26 * 3 the finite-difference rebuilds cost.

col_derivint, optional

Non-zero means Dfun returns the transpose: (neq, neq) read as jac[j, i] = d f_i / d y_j, and with ml/mu the packed array transposed to (neq, ml + mu + 1). Inert with Dfun=None. Must be a compile-time constant inside @njit.

full_outputbool, optional

False (default) returns y. True returns (y, infodict). Must be a compile-time constant inside @njit.

ml, muint or None, optional

Half-bandwidths of the Jacobian, counting the sub- and super-diagonals and excluding the main diagonal. Setting either selects LSODA’s banded Jacobian, jt=5 with Dfun=None and jt=4 with a Dfun; the unset one becomes 0. None (default) on both keeps the full Jacobian. A negative value means the same as None.

Cost, on an 80-point heat equation by method of lines whose Jacobian is tridiagonal, rtol=1e-6 atol=1e-9: LSODA rebuilds the Jacobian 6 times and a full finite-difference rebuild spends neq = 80 right-hand-side evaluations each, against ml + mu + 1 = 3 for the banded one. nfe 597 against 135.

rtol, atolfloat, sequence or None, optional

Relative and absolute tolerances. None (the default) means 1.49012e-8, sqrt(finfo(float64).eps). Either may be a scalar or a vector of length neq; the four combinations map to LSODA’s itol 1, 2, 3 and 4. A sequence of any other length raises ODEpackError, and rank >= 2 raises ODEpackError("Error converting relative tolerance.").

tcritfloat, sequence or None, optional

Critical times the integrator must not step past, in the order it will meet them. The first is in force from the first leg; the integrator advances to the next each time an output time passes the current one, and steps freely once the sequence is exhausted.

h0float, optional

First step size. 0.0 (default) lets LSODA choose.

hmax, hminfloat, optional

Step-size bounds. 0.0 (default) means no bound.

ixprint, optional

Print a message on each method switch.

mxstepint, optional

Maximum internal steps between two output times. 0 (default) means LSODA’s internal 500.

mxhnilint, optional

Maximum “step size too small” messages. 0 (default) means LSODA’s internal 10.

mxordn, mxordsint, optional

Maximum order for the Adams (12) and BDF (5) methods.

printmessgint, optional

Non-zero asks for the success message as an ODEintWarning. An abnormal exit warns whatever this is set to.

tfirstbool, optional

False (default) selects the f(y, t, *args) callback order. True selects f(t, y, *args). Must be a compile-time constant inside @njit.

success_outbool, optional

True appends a success bool to the return. False is the default. Must be a compile-time constant inside @njit.

Returns:

yfloat64 ndarray, shape (len(t), len(y0))

Solution at each requested time, time-major. y[0] is y0. On an abnormal exit the rows up to and including the failing leg hold real values and the rest are 0.0.

infodictInfoDict

Only with full_output=True. The 13 keys reached as attributes: hu, imxer, leniw, lenrw, message, mused, nfe, nje, nqu, nst, tcur, tolsf, tsw. The nine vector fields have length len(t) - 1, one entry per integration leg. Entries past a failure are 0.0. tolsf is written by LSODA only on an “excess accuracy” exit, so on a successful run it reads 0.0.

successbool

Only with success_out=True. istate == 2.

Raises:

TypeError

func or Dfun is not a plain @njit function (an integer or other non-@njit value); or its arity does not fit args.

ODEpackError

args is not a tuple, or a tolerance has the wrong length or rank.

RuntimeError

The right-hand side returned the wrong number of derivatives, or Dfun returned the wrong shape.

ValueError

y0 or t has rank >= 2, or t is not monotonic.

Warns:

ODEintWarning

On every abnormal exit, and on a normal one under printmessg.

See also

scipy.integrate.odeint

The scipy routine this mirrors.

Notes

success_out has no scipy counterpart.

ODEpackError derives straight from Exception, where scipy raises its private scipy.integrate._odepack.error, also derived from Exception.

On an abnormal exit the tail of y past the failing leg is 0.0, and the infodict vector entries past that leg are 0.0; scipy leaves both uninitialized. tolsf reads 0.0 on a successful run, where scipy returns whatever the heap held.

The right-hand side is called once before the run and the length of what it returns is checked, which is what scipy checks inside its callback on every call. Dfun’s shape is checked the same way. A callback whose shape CHANGES mid-run is caught by scipy on the offending call and is not caught here.

Dfun’s arity is checked before the run as well. scipy reaches a wrong arity only when LSODA calls the Jacobian, which on a non-stiff problem never happens.

Two refusals arrive earlier from inside @njit than the check order scipy runs. A func or Dfun that is not a plain @njit function, and an args entry that is neither a real number nor an array of them, are both settled when the call compiles, so a call carrying one of them and a non-monotonic t reports the refusal where scipy reports the t error. Called from python, both follow scipy’s order.

infodict is a namedtuple where scipy returns a dict, so info['nfe'] is spelled info.nfe.

ixpr=1 prints LSODA’s method-switch messages on Fortran unit 6. scipy accepts the argument and prints nothing. Unit 6 is not reachable through contextlib.redirect_stdout, and concurrent solves under prange interleave on it, so ixpr=0, the default, writes nothing.

Prange-safe. The callback address, the args pointer, the Jacobian address and the tfirst flag travel through Fortran module variables, and all slots are !$omp threadprivate, so each OS thread reads its own copy.

Examples

>>> import numpy as np
>>> from numba import njit
>>> import scijit.integrate as si
>>> @njit
... def rhs(y, t, k):                  # scipy's order, y first
...     out = np.empty(2)              # y'' = -k y
...     out[0] = y[1]
...     out[1] = -k * y[0]
...     return out
>>> @njit
... def run():
...     return si.odeint(rhs, np.array([1.0, 0.0]),
...                      np.array([0.0, 0.5, 1.0]), (4.0,))
>>> run()
array([[ 1.        ,  0.        ],
       [ 0.54030231, -1.68294195],
       [-0.4161468 , -1.81859487]])

full_output=True adds scipy’s diagnostics as a namedtuple, so info['nfe'] is spelled info.nfe:

>>> y, info = si.odeint(rhs, np.array([1.0, 0.0]),
...                     np.array([0.0, 0.5, 1.0]), (4.0,),
...                     full_output=True)
>>> info.message
'Integration successful.'