scijit.integrate.solve_ivp

scijit.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, vectorized=False, args=None, rtol=0.001, atol=1e-06, max_step=inf, mxstep=0, npoints=100, terminal=None, direction=None, first_step=None, **options)

Solve an initial value problem for a system of ODEs.

Integrates dy/dt = fun(t, y) from t_span[0] to t_span[1], starting at y0.

Parameters:
fun@njit function f(t, y) or f(t, y, args)

Right-hand side, t-first. y is a 1-D float64 array and the return must be a new 1-D float64 array of the same length. The three-argument form receives args, on every method. A plain @njit function, on every method.

t_span2-tuple or length-2 array of float

(t0, tf). tf < t0 integrates backwards, on every method. Without t_eval, tf == t0 returns the initial state twice; with t_eval it reports nothing.

y0float64 or complex128 array, shape (n,)

Initial state. Must be 1-D. A complex128 y0 integrates in the complex domain on 'RK45', 'RK23' and 'DOP853', and the result’s y and sol are complex128. 'LSODA' raises there. The dtype is read when the call compiles, so the result type never depends on a value.

methodstr, optional

'RK45' (default), 'RK23', 'DOP853' or 'LSODA'. Case sensitive. Inside @njit it must be a string literal, since it selects which solver is compiled in. 'Radau' and 'BDF' raise.

t_evalfloat64 array or None, optional

Times to report at, inside t_span and STRICTLY ordered in the direction of travel, so a repeat is refused. None means the solver’s own steps, on every method. An empty t_eval reports nothing on every method, which is a different input from None. Must be 1-D.

dense_outputbool, optional

Compile-time literal. True adds a callable sol to the result and changes the return type. Available on every method.

events@njit function g(t, y) or g(t, y, args), optional

Event functions, evaluated together and returning a float64 array of ng values, one per event. A root of any of them is reported. ng is learnt by calling the function once at (t_span[0], y0). A single event may return a scalar. The three-argument form receives args, and the arity must agree with args in both directions or the call raises TypeError.

vectorizedbool, optional

Whether the right-hand side is written to take a block of states: True means it is called as fun(t, y[:, None]), shape (n, 1), and its result ravelled. Accepted on every method and on both callback spellings.

Inside @njit it must be a literal, since it selects which adapter is built when the call compiles. A variable raises.

argsfloat64 array or None, optional

Extra parameters, one flat float64 buffer. Reaches fun, jac and events on every method. The callback has to be able to receive it: a three-argument f(t, y, args). A two-argument f(t, y) given args raises TypeError, and a three-argument one given no args raises TypeError, in both directions. An EMPTY args carries no parameters, so it counts as absent in both: f(t, y) runs and f(t, y, args) raises, from Python and from inside @njit.

rtol, atolfloat or float64 array, optional

Tolerances, 1e-3 and 1e-6. Either may be a scalar or an array with one entry per component of y0, on every method. atol is shape-checked and rtol is not, so a length-1 atol array raises where len(y0) > 1 and a length-1 rtol array broadcasts.

An rtol below 100 * eps, 2.220446049250313e-14, is raised to it and a UserWarning says so, on every method. atol = 0 is legal and makes the tolerance purely relative; a negative atol raises. A purely relative tolerance collapses the step size wherever a component passes through zero.

max_stepfloat, optional

Largest step allowed. np.inf by default.

mxstepint, optional

'LSODA' only: step limit per output interval. 0 selects LSODA’s own 500. On the RK methods it is ignored with a UserWarning.

npointsint, optional

'LSODA' with t_eval=None only: how many uniformly spaced output times to produce, INSTEAD of the solver’s own steps. The default, 100, reports the steps; any other value reports the grid. On the RK methods it is ignored with a UserWarning.

terminalint64 array, shape (ng,), optional

n stops the integration at that event’s n-th occurrence and 0 never stops it. None (the default) makes every event non-terminal. A negative or fractional entry raises ValueError.

With two terminal events crossing in one step the earlier root in TIME stops the run, not the lower event index.

directionfloat64 array, shape (ng,), optional

Positive keeps upward crossings only, negative downward only, zero keeps both. None (the default) keeps both for every event.

first_stepfloat or None, optional

Initial step size. None, the default, lets the solver choose one. Reaches all four methods.

**optionsdict, optional

The keyword-only 'LSODA' controls min_step, lband, uband and jac, listed under Other Parameters. Any other keyword raises TypeError.

Returns:
resOdeResult, or OdeResultDense when dense_output=True

res.t shape (m,); res.y shape (n, m), component-major; res.nfev right-hand-side evaluations; res.njev Jacobian evaluations; res.nlu LU decompositions; res.status 0 when the end of t_span was reached and -1 when a step failed; res.message and res.success. res.sol is present only under dense_output=True and is callable as res.sol(t) for a scalar or an array of times. Without events, res.t_events and res.y_events are None; with events both are lists of arrays, one entry per event, shapes (k_i,) and (k_i, n).

Other Parameters:
min_stepfloat, optional

'LSODA' only: smallest step allowed. 0.0, the default, means no floor. On the RK methods it is ignored with a UserWarning.

lband, ubandint or None, optional

'LSODA' only: ignored with a UserWarning on the RK methods. Half-bandwidths of the Jacobian, counting the sub- and super-diagonals and excluding the main diagonal, so d f_i / d y_j is taken as zero outside i - lband <= j <= i + uband. Setting either one selects LSODA’s banded Jacobian and the other becomes 0. None on both, the default, keeps the full Jacobian.

LSODA rebuilds the Jacobian by finite differences unless jac is given, and a full rebuild costs len(y0) right-hand-side evaluations against lband + uband + 1 for a banded one. Measured on an 80-point heat equation by method of lines, whose Jacobian is tridiagonal, rtol=1e-6 atol=1e-9 over [0, 0.5]: 6 rebuilds, nfev 597 without the bands and 135 with lband=uband=1, values agreeing to 0.000e+00, and 2.85x less wall time over 20 solves.

A value outside 0 <= lband, uband <= len(y0) - 1 reaches ODEPACK, which reports “Illegal input detected (internal error)” and returns success=False.

A band NARROWER than the Jacobian is legal and expensive. ODEPACK takes it as an approximation, and the corrector then needs many more iterations. Measured on the same 40-point problem, lband=1 alone against the correct lband=uband=1: nfev 4523 and njev 640 against 265 and 17, landing 1.083e-06 away. This can exhaust mxstep, which ODEPACK counts PER OUTPUT INTERVAL, and the run then reports “Excess work done on this call”. Raise mxstep: at mxstep=100000 every narrow band above succeeds.

jac@njit function j(t, y) or j(t, y, args), optional

'LSODA' only: ignored with a UserWarning on the RK methods. The Jacobian d f / d y. Without it LSODA builds one by finite differences, at the cost above.

Returns (n, n) with jac[i, j] = d f_i / d y_j. With lband/uband it returns the packed banded form instead, (lband + uband + 1, n) with jac_packed[uband + i - j, j] = d f_i / d y_j. A wrong shape raises ValueError, from one evaluation at (t_span[0], y0).

The three-argument form receives args, as fun’s does.

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

Raises:
ValueError

method outside 'RK45', 'RK23', 'DOP853', 'LSODA'; 'Radau' or 'BDF'; a t_span without exactly two elements; a nan in t_span; a y0 that is not 1-D or not finite; a t_eval that is not 1-D, leaves t_span or is not strictly ordered; a non-positive max_step; a non-positive first_step, or one exceeding the span; a negative min_step; a negative atol or one of the wrong shape; an rtol of the wrong length; a terminal entry that is negative or fractional; a terminal or direction whose length is not ng; a jac of the wrong shape; a right-hand side returning the wrong length; a complex y0 on 'LSODA', 'Radau' or 'BDF'; a complex y0 together with events; and, on 'LSODA', a right-hand side that never wrote ydot.

TypeError

A nested or 2-D t_span; an arity mismatch between args and fun or events, in either direction.

RuntimeError

An event root that Brent’s method did not bracket in 100 iterations.

See also

scipy.integrate.solve_ivp

The scipy routine this mirrors.

Notes

Every method is prange-safe.

'RK45', 'RK23' and 'DOP853' hold no module state. 'LSODA' reaches Fortran, and its callback slot is !$omp threadprivate, so each thread reads its own copy.

The solution object. res.sol(t) takes a float or a float64 array. A Python list is not an array: numba types it as a reflected list, which no array guard admits, so res.sol([0.1, 0.2]) raises TypeError. Wrap it, res.sol(np.array([0.1, 0.2])).

A solution object reaches compiled code as an argument, as a module global, or captured in a closure; a registered jitclass instance is a compile-time constant.

Reading a shared solution object from a prange loop is safe and reproduces the serial result exactly. Writing a field from several threads is a data race.

Differences from scipy.

  • sol is absent from the result unless dense_output=True, where scipy always carries the field and sets it to None. So hasattr(res, 'sol') answers what res.sol is None answers in scipy.

  • events is one function returning ng values, where scipy takes a callable or a list of them, and terminal / direction are arrays passed alongside rather than attributes set on the function. A single event may return a scalar, as scipy’s do.

  • nfev after a terminal event covers the WHOLE span. scipy leaves its stepping loop at the event, so its count stops there; this integrates the span and trims the history afterwards.

  • jac, lband, uband, min_step, mxstep and npoints on 'RK45', 'RK23' and 'DOP853' are ignored with one UserWarning per call, “The following arguments have no effect for a chosen solver: mxstep, npoints, min_step, jac, lband, uband.”, which is character for character scipy’s on the same call. The answer is the one the same call gives without them.

    min_step, jac, lband and uband warn when the argument was PASSED, whatever its value, so solve_ivp(..., 'RK45', jac=None) warns. mxstep and npoints warn when they carry a value other than their default, so solve_ivp(..., 'RK45', mxstep=0) does not. Neither has a scipy counterpart.

    The names are listed in signature order. scipy lists them in the order the caller passed them, which a compiled call site does not preserve.

  • An unrecognised keyword argument raises TypeError. scipy absorbs it into **options and names it in the same UserWarning.

  • jac is a plain @njit function taking args as one flat float64 buffer, matching fun, rather than scipy’s callable with unpacked parameters.

  • args reaches the callback as one flat float64 buffer rather than unpacked into separate parameters. scipy calls f(t, y, *args); this calls f(t, y, args) with everything packed into one array, because numba has no *args. Which methods accept it no longer differs: every one does, as in scipy.

  • An arity mismatch between fun and args is refused in both directions, as in scipy: args given to an f(t, y), and args omitted for an f(t, y, args). An empty args is not a mismatch either way. The exception is TypeError, scipy’s, from python, and TypingError from inside @njit. One case differs there: an array’s length is not known when the call compiles, so an empty args reads as given inside @njit and as absent from python.

  • tf == t0 with a t_eval, and an empty t_eval on any span, report nothing, as in scipy. scipy returns t and y as empty Python lists there and this returns empty arrays, shapes (0,) and (n, 0).

  • events and a complex y0 together raise. scipy carries both: it stores the state at each event root, and here that store is float64. Every other combination of a complex y0 with t_eval, dense_output and args is supported.

  • 'Radau' and 'BDF' raise. They are planned.

  • The result is a namedtuple with scipy’s field names, so res.y works and res['y'] does not.

  • 'LSODA' with a t_eval and neither dense_output nor events reports through a run that stops at each requested time, the route scipy’s odeint takes, where scipy’s solve_ivp steps freely and interpolates. This returns scipy’s odeint answer. Asking for dense_output or events takes the stepping route instead, matching scipy’s solve_ivp.

Examples

>>> import numpy as np
>>> from numba import njit
>>> import scijit.integrate as si
>>> @njit
... def rhs(t, y):                       # y'' = -4 y
...     out = np.empty(2)
...     out[0] = y[1]
...     out[1] = -4.0 * y[0]
...     return out
>>> @njit
... def run():
...     return si.solve_ivp(rhs, (0.0, 1.0), np.array([1.0, 0.0]),
...                         'RK45', np.array([0.0, 0.5, 1.0]))
>>> res = run()
>>> res.y.shape                          # (n_states, n_times)
(2, 3)
>>> res.y[0]
array([ 1.        ,  0.54038693, -0.41629169])
>>> res.success
True

dense_output=True adds res.sol, callable on a scalar or an array of times:

>>> res = si.solve_ivp(rhs, (0.0, 1.0), np.array([1.0, 0.0]),
...                    'RK45', None, True)
>>> res.sol(0.5)
array([ 0.54038693, -1.68337531])
>>> res.sol(np.array([0.25, 0.5]))
array([[ 0.8775798 ,  0.54038693],
       [-0.95885456, -1.68337531]])

events is one @njit function returning one value per event. A root of any of them is reported in res.t_events:

>>> @njit
... def hits_zero(t, y):
...     out = np.empty(1)
...     out[0] = y[0]
...     return out
>>> res = si.solve_ivp(rhs, (0.0, 4.0), np.array([1.0, 0.0]), 'RK45',
...                    events=hits_zero)
>>> res.t_events[0]
array([0.78528433, 2.35559604, 3.92597172])

terminal stops the integration at the first root, where scipy sets a .terminal attribute on the event function:

>>> res = si.solve_ivp(rhs, (0.0, 4.0), np.array([1.0, 0.0]), 'RK45',
...                    events=hits_zero, terminal=np.array([1]))
>>> res.t_events[0]
array([0.78528433])
>>> res.status
1
>>> res.message
'A termination event occurred.'