`vp.ode.solve`

Solve an initial value problem y' = f(t, y, p), stiff or not, with the right-hand side in Viv.

Part of the distribution vipyr; it runs from a checkout of the repository. Planned on PyPI as pip install vipyr.

vp.ode.solve is scipy.integrate.solve_ivp with the integrator and the right-hand side compiled together: Viv's Dormand-Prince, Radau IIA and BDF (lib/ode/) call your rhs directly, with no Python in the loop. Dormand-Prince follows SciPy's RK45 step for step. Radau and BDF take SciPy's number of steps to within two and its number of right-hand-side calls to within 2% (python/tests/test_ode.py), so tolerances carry over unchanged.

The right-hand side is Viv because the integrator calls it at every stage of every step, and a call into Python there is what makes solve_ivp slow on small systems. Give it as

a problem of lib/ode/problems.viv by name: "robertson", "vdpol", "hires", "oregonator", "e5", "kepler", "lotka", "ball", "sir", "pk";
a def of your own file: "model.viv:decay", or k.decay after k = vp.load("model.viv");
a Python function in py2viv's numeric subset (tools/py2viv/SUBSET.md), translated to Viv once and cached.

In every form the def is def f(t: Num, y: Arr, p: Arr) -> Arr: the parameters arrive in p (params=), not in variables captured from Python, because a compiled Viv def cannot close over Python state. An events def has the same signature and returns one g value per event. The first call for a problem writes a small driver .viv and compiles it (a few seconds); later calls, in this process or the next, load it from vipyr's cache.

method="auto" starts with Dormand-Prince and switches to Radau IIA when Hairer's stiffness test (15 steps in a row with h |lambda| > 3.25) fires; stiff_at says when.

Signature

vp.ode.solve(rhs, t_span, y0, *, method='auto', rtol=0.001, atol=1e-06, params=None, t_eval=None, events=None, direction=None, terminal=None, jac=None, first_step=None, max_step=None, max_steps=1000000, fast=False)

Parameters

rhs : str | function, required: The right-hand side f(t, y, p), a Viv def: a problem of lib/ode/problems.viv by name, 'file.viv:def', a def of a module from vp.load, or a Python function py2viv translates.
t_span : (float, float), required: The interval (t0, t1); t1 < t0 integrates backward.
y0 : ndarray[float64], required: The initial state, length n.
method : str, default 'auto': 'dopri5' (or 'rk45'), 'radau', 'bdf', or 'auto': Dormand-Prince until its stiffness test fires, then Radau IIA.
rtol : float, default 0.001: Relative tolerance; raised to 100 eps if smaller.
atol : float | ndarray[float64], default 1e-06: Absolute tolerance, one number or one per component.
params : ndarray[float64], default None: The parameter vector p every call of rhs, events and jac receives; one compiled solver serves every p.
t_eval : ndarray[float64], default None: Times to report the solution at (dense output); None reports every step.
events : str | function, default None: A def of the rhs's file returning one g value per event; an event is a zero of g.
direction : ndarray[float64], default None: Per event: +1 rising zeros only, -1 falling only, 0 both.
terminal : ndarray[bool], default None: Per event: stop the integration at its first zero.
jac : str, default None: None: exact, by reverse-mode grad of the rhs; 'fd': forward differences; or the name of a def of the rhs's file returning the n by n Jacobian.
first_step : float, default None: The first step; None chooses it (Hairer's rule).
max_step : float, default None: The largest step allowed; None for no limit.
max_steps : int, default 1000000: Stop with status -2 after this many accepted steps.
fast : bool, default False: Compile with the fast flags (fused multiply-add): faster, last bits may differ.

Returns

A frozen OdeSolution. Iterating it gives the fields below in order; it also has params, tests, log and check().

t : ndarray[float64]: The output times.
y : ndarray[float64]: The solution, shape (n, len(t)) as in SciPy: y[i] is component i.
t_events : list[ndarray]: Per event, the times it fired (None without events).
y_events : list[ndarray]: Per event, the states where it fired.
status : int: 1: reached t1; 2: a terminal event; -1: step size underflow; -2: max_steps; -3: the solution stopped being finite.
message : str: The status in words.
success : bool: Status > 0.
stats : dict: Nfev (rhs evaluations), njev (Jacobians), nlu (LU factorizations), steps (accepted), rejected.
stiff_at : float: With method='auto', the time it changed to Radau IIA (None if it did not).

Raises

ValueError: A method, interval, tolerance or array shape the solver cannot take; vipyr._ode.OdeProblemError (a ValueError) when the rhs, events or jac cannot be made into a Viv problem, or py2viv refuses a Python rhs.
vipyr.VivCompileError: A SyntaxError: the rhs's .viv file does not compile.
TypeError: An argument of the wrong type, or an unknown keyword.
vipyr.VivBuildError: An ImportError: no C compiler, or the program's C does not build.

Notes

Dormand-Prince 5(4) ("dopri5") is a seven-stage explicit Runge-Kutta pair: a fifth-order solution $y_{n + 1} = y_{n} + h \sum_{i} b_{i} k_{i}$ and an embedded fourth-order one whose difference is the error estimate. The last stage is $f (t_{n + 1}, y_{n + 1})$ , the next step's first ("first same as last"), so a step costs six evaluations. A step is accepted when

err = (\frac{1}{n} i \sum (\frac{e _{i}}{atol _{i} + rtol max ( ∣ y _{n, i} ∣ , ∣ y _{n + 1, i} ∣ )})^{2})^{1/2} \leq 1,

and the next step is $h \cdot min (10, 0.9 err^{- 1/5})$ (at least $0.2 h$ after a rejection). The dense output is Hairer's fifth-order continuous extension, built from the stages already computed.

Radau IIA ("radau") is the three-stage collocation method at the right Radau points $c = ((4 - 6) /10, (4 + 6) /10, 1)$ , order 5, L-stable and stiffly accurate. The stage increments $Z$ solve $Z = h (A \otimes I) F (Z)$ by simplified Newton iterations with the Jacobian $J$ at the start of the step. Diagonalizing $A^{- 1} = T Λ T^{- 1}$ splits the $3 n \times 3 n$ Newton matrix into one real system $(γ / h) I - J$ and one complex system $((α + i β) / h) I - J$ , so a step costs one real and one complex LU (written here as a real $2 n \times 2 n$ system). The error estimate is Hairer and Wanner's (IV.8, eq. 8.20), $((γ / h) I - J)^{- 1} (f (t_{n}, y_{n}) + \sum_{i} e_{i} Z_{i} / h)$ , and the step size follows Gustafsson's predictive controller. A Newton iteration that contracts slower than its rate test allows is stopped early; it refreshes a stale Jacobian first and halves the step second.

BDF ("bdf") is variable-order (1 to 5) backward differentiation in Shampine and Reichelt's quasi-constant-step form: the history is kept as backward differences $D^{j} y_{n}$ on an equally spaced grid and rescaled by their matrices $R$ and $U$ when the step changes. The formulas are the NDFs of MATLAB's ode15s, with $κ = (- 0.185, - 1/9, - 0.0823, - 0.0415, 0)$ : each order solves

(I - \frac{h}{α _{k}} J) Δ d = \frac{h}{α _{k}} f (t_{n + 1}, y^{(0)} + d) - ψ - d

for the correction $d$ to the predictor, with one LU of $I - (h / α_{k}) J$ kept while $h$ and $k$ stay put. After $k + 1$ equal steps the errors of orders $k - 1$ , $k$ and $k + 1$ are compared and the order moves to the one that allows the largest step.

Jacobians. By default the driver writes one reverse-mode grad of the rhs per row, so $J$ is exact to rounding and costs about $n$ evaluations of the rhs. jac="fd" takes forward differences instead (n evaluations, an error of order $ϵ$ ); jac="name" calls a def of yours.

Events. After each accepted step every event's $g$ is evaluated at the new point; a sign change (in the direction asked for) is located on the step's dense output by the Illinois variant of regula falsi to four units in the last place of $t$ . A terminal event ends the integration at its root, and the outputs stop there.

Batches. vp.ode.batch runs the same compiled solver over the rows of params and y0 in one call, on one output grid. Each row is an independent solve with its own steps, bit for bit the solve vp.ode.solve makes with that row.

Limits. Dense matrices only: the LU is $O (n^{3})$ , fine up to a few hundred unknowns, wrong for a large sparse system. No mass matrix or DAE form. BDF fails on Van der Pol with $μ = 1 0^{6}$ at its first relaxation jump (as SciPy's does); Radau solves it. E5 (Hairer's test set) is not reproduced by these methods or SciPy's at any rtol and atol tried: the reference's $y_{2}$ and $y_{3}$ differ by $1 0^{- 25}$ at $1 0^{- 22}$ .

References

J. R. Dormand and P. J. Prince, A family of embedded Runge-Kutta formulae, Journal of Computational and Applied Mathematics 6, 19-26 (1980). doi:10.1016/0771-050X(80)90013-3
E. Hairer, S. P. Norsett and G. Wanner, Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed., Springer (1993). doi:10.1007/978-3-540-78862-1
E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd ed., Springer (1996). doi:10.1007/978-3-642-05221-7
L. F. Shampine and M. W. Reichelt, The MATLAB ODE Suite, SIAM J. Sci. Comput. 18, 1-22 (1997). doi:10.1137/S1064827594276424
F. Mazzia and C. Magherini, Test Set for Initial Value Problem Solvers, release 2.4, University of Bari (2008): the reference solutions of Robertson, Van der Pol, HIRES and the Oregonator. https://archimede.uniba.it/~testset/
P. Virtanen et al., SciPy 1.0: fundamental algorithms for scientific computing in Python, Nat. Methods 17, 261-272 (2020): scipy.integrate.solve_ivp, whose step logic lib/ode follows. doi:10.1038/s41592-019-0686-2

Examples

import vipyr as vp
# Robertson's kinetics, stiff, to t = 1e5 by Radau IIA with the exact Jacobian
r = vp.ode.solve("robertson", (0.0, 1e5), [1.0, 0.0, 0.0], method="radau", rtol=1e-8, atol=1e-12,
                 params=[0.04, 1e4, 3e7], t_eval=[1.0, 1e3, 1e5])
print(r.y[:, -1], r.stats["steps"])
assert r.check()

[1.78659211e-02 7.27475147e-08 9.82134006e-01] 572

A bouncing ball's landing as a terminal event, located on the dense output.

import numpy as np
import vipyr as vp
r = vp.ode.solve("ball", (0.0, 10.0), [10.0, 0.0], method="dopri5", params=[9.81],
                 events="ball_hits", direction=[-1], terminal=[True])
print(r.message, r.t_events[0], np.sqrt(2 * 10.0 / 9.81))

A termination event occurred. [1.42784312] 1.4278431229270645

Van der Pol with mu = 1000 (Hairer's scaled form, eps = 10⁻⁶) by auto: it starts explicit and changes to Radau IIA.

import vipyr as vp
r = vp.ode.solve("vdpol", (0.0, 2.0), [2.0, 0.0], method="auto", rtol=1e-6, atol=1e-6, params=[1e-6])
print(r.y[:, -1], r.stiff_at is not None, r.stats)

[ 1.70616772 -0.89280971] True {'nfev': 7442, 'njev': 205, 'nlu': 578, 'steps': 895, 'rejected': 26}

Each example ran as written when these pages were generated, and the output shown is what it printed.

Checked against

Reference: The stiff test set's published reference solutions (Mazzia and Magherini, Bari, 2008): Robertson to 10¹¹, Van der Pol (mu = 1000), HIRES and the Oregonator, by Radau IIA (largest relative error 5.2 × 10⁻¹⁰ at rtol 10⁻⁸) and BDF (1.6 × 10⁻⁶, as SciPy's BDF); Kepler, Lotka-Volterra, events, dense output and a batch against closed forms (lib/ode/stiff_reference.viv, lib/ode/nonstiff_events.viv); step counts against SciPy's solve_ivp (python/tests/test_ode.py).

Source

Viv: lib/ode/solve.viv
Viv: lib/ode/dopri5.viv
Viv: lib/ode/radau.viv
Viv: lib/ode/bdf.viv
Viv: lib/ode/core.viv
Viv: lib/ode/lu.viv
The defs it reaches and the calls between them: map of the code

vp.ode.solve