mathematicskit.integrators#

Shared ODE integrators used by ode_dynamics and pde: fixed-step explicit Euler and RK4, Adams-Bashforth multistep methods, symplectic leapfrog and Yoshida, adaptive Dormand-Prince, hand-rolled backward Euler and scipy’s Radau/BDF for stiff systems, collocation for two-point boundary-value problems, and each method’s linear stability region. The “Numerical integration” and “Stiffness and boundary-value problems” sections of the ode_dynamics gallery demonstrate each one, and Breakthroughs in Dynamical Systems covers their history.

Numba-accelerated (when Numba is installed) ODE integrators shared across mathematicskit subpackages.

Right-hand-side / force callbacks use the f(state_or_pos, t, params) -> ndarray convention throughout, so a single compiled callback can be reused across systems with different parameter values without recompiling (see mathematicskit.integrators.fixed_step).

stability_function(), is_absolutely_stable() and adams_bashforth_boundary_locus() describe each explicit and implicit method’s linear stability region in the z = lam * dt plane (see mathematicskit.integrators.stability).

mathematicskit.ode_dynamics builds its phase portraits, bifurcation diagrams, and Poincare sections on top of this shared module rather than reimplementing per-domain integrators, exactly as physicskit’s classical/core/integrators.py wraps its own shared integrators for that subpackage’s njit-callback calling convention.

The integrator kernels are compiled per process and deliberately not cached to disk (no cache=True). Each one takes the callback itself as an argument, so a cache entry is keyed on that callback’s dispatcher. Numba cannot re-save an index holding keys whose dispatchers belonged to an earlier process, so the second script to use a new callback would crash with ReferenceError: underlying object has vanished. The kernels that take only arrays and scalars (e.g. the rhs functions in mathematicskit.ode_dynamics and mathematicskit.pde) stay cached.

njit() is re-exported here for writing those callbacks: it is numba.njit() when Numba is installed (pip install mathematicskit[fast]) and a no-op decorator otherwise, with HAS_NUMBA reporting which.

class mathematicskit.integrators.BVPResult(x, y, sol, success, message, max_rms_residual)[source]#

Bases: object

Solution of a two-point boundary-value problem (see collocation_bvp()).

Parameters:

  • x (NDArray[float64])

  • y (NDArray[float64])

  • sol (Callable[[ArrayLike], NDArray[float64]])

  • success (bool)

  • message (str)

  • max_rms_residual (float)

max_rms_residual: float#

Largest per-interval RMS relative residual of the collocation spline.

Type:

float

message: str#

Solver termination message.

Type:

str

sol: Callable[[ArrayLike], NDArray[float64]]#

sol(x) -> ndarray (len(x), dim), the continuous cubic-spline solution evaluated anywhere in the interval.

Type:

callable

success: bool#

Whether the requested tolerance was met.

Type:

bool

x: NDArray[float64]#

Final (refined) mesh.

Type:

ndarray, shape (n_nodes,)

y: NDArray[float64]#

Solution at the mesh nodes.

Type:

ndarray, shape (n_nodes, dim)

mathematicskit.integrators.adams_bashforth_boundary_locus(order, n_points=721)[source]#

Boundary locus \(z(\theta) = \rho(e^{i\theta})/\sigma(e^{i\theta})\) of Adams-Bashforth.

Every point where the characteristic polynomial has a root on the unit circle lies on this curve, so the stability region’s boundary is part of it (for orders 1 and 2 it is all of it).

Parameters:

  • order (int) – 1 to 4 (order 1 is explicit Euler).

  • n_points (int) – Number of \(\theta\) samples on \([0, 2\pi]\).

Return type:

ndarray

Returns:

ndarray of complex, shape (n_points,)

Examples

>>> import numpy as np
>>> locus = adams_bashforth_boundary_locus(1, n_points=5)
>>> np.allclose(locus, np.exp(1j * np.linspace(0, 2 * np.pi, 5)) - 1.0)
True
>>> # AB2's locus crosses the negative real axis at z = -1 (theta = pi).
>>> round(float(adams_bashforth_boundary_locus(2, n_points=3)[1].real), 12)
-1.0
mathematicskit.integrators.adams_bashforth_integrate(rhs, state0, t0, dt, n_steps, params, order=4)[source]#

Integrate n_steps of the order-order Adams-Bashforth method.

The first order - 1 steps are RK4 steps, which supply the past slopes the multistep formula needs. Order 1 is explicit Euler.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n_steps + 1,)) – Time at each step, starting at t0.

  • states (ndarray of float, shape (n_steps + 1, dim)) – State at each step, starting at state0.

Raises:

ValueError – If order is not 1, 2, 3 or 4.

mathematicskit.integrators.collocation_bvp(rhs, bc, x, state_guess, params, tol=1e-06, max_nodes=10000)[source]#

Solve state' = rhs(state, x, params) subject to bc(state(a), state(b)) = 0.

Parameters:
Return type:

BVPResult

Returns:

BVPResult – The solution; check success before trusting it.

Examples

y'' = -y with y(0) = 0, y(pi/2) = 1 has solution sin x:

>>> import numpy as np
>>> def rhs(state, x, params):
...     return np.array([state[1], -state[0]])
>>> def bc(ya, yb, params):
...     return np.array([ya[0], yb[0] - 1.0])
>>> x = np.linspace(0.0, np.pi / 2, 11)
>>> res = collocation_bvp(rhs, bc, x, np.zeros((11, 2)), np.zeros(0))
>>> res.success, bool(np.allclose(res.sol([0.5])[0, 0], np.sin(0.5), atol=1e-6))
(True, True)
mathematicskit.integrators.dopri5_integrate(rhs, state0, t0, t_end, dt0, params, rtol=1e-06, atol=1e-09, dt_min=1e-12, dt_max=1000000.0, safety=0.9, max_steps=100000)[source]#

Integrate forward from t0 to t_end with adaptive step-size control.

Steps are accepted or rejected from the embedded RK5(4) error estimate against an atol + rtol * |state| tolerance (the same convention as scipy.integrate.solve_ivp()), with the step size adjusted after every attempt. Only forward integration (t_end > t0, dt0 > 0) is supported.

Parameters:

  • rhs (Callable[[NDArray[double], float, NDArray[double]], NDArray[double]]) – Numba-jitted right-hand-side function rhs(state, t, params) -> ndarray.

  • state0 (NDArray[double]) – Initial state vector.

  • t0 (float) – Integration interval, with t_end > t0.

  • t_end (float) – Integration interval, with t_end > t0.

  • dt0 (float) – Initial step size to attempt.

  • params (NDArray[double]) – Parameter vector passed through to rhs.

  • rtol (float) – Relative and absolute error tolerances.

  • atol (float) – Relative and absolute error tolerances.

  • dt_min (float) – Step-size bounds; a step is accepted once dt shrinks to dt_min regardless of its error estimate, to guarantee progress.

  • dt_max (float) – Step-size bounds; a step is accepted once dt shrinks to dt_min regardless of its error estimate, to guarantee progress.

  • safety (float) – Safety factor applied to the step-size update.

  • max_steps (int) – Upper bound on the number of attempted steps (accepted or rejected), to guarantee termination.

Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n_accepted + 1,)) – Times of the accepted steps, starting at t0 and ending at t_end.

  • states (ndarray of float, shape (n_accepted + 1, dim)) – State at each accepted step.

mathematicskit.integrators.dopri5_step(rhs, state, t, dt, params)[source]#

Single embedded Dormand-Prince RK5(4) step, with an error estimate.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • state_new (ndarray of float, shape (dim,)) – The 5th-order state estimate at t + dt.

  • error (ndarray of float, shape (dim,)) – Difference between the 5th- and embedded 4th-order estimates, for step-size control (see dopri5_integrate()).

mathematicskit.integrators.euler_integrate(rhs, state0, t0, dt, n_steps, params)[source]#

Integrate n_steps of explicit Euler starting from state0.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n_steps + 1,)) – Time at each step, starting at t0.

  • states (ndarray of float, shape (n_steps + 1, dim)) – State at each step, starting at state0.

mathematicskit.integrators.euler_step(rhs, state, t, dt, params)[source]#

Single explicit (forward) Euler step: state + dt * rhs(state, t).

Follows the tangent line for one step (Euler 1768). The local error is \(O(dt^2)\), so the global error over a fixed interval is \(O(dt)\). For y' = lam * y one step multiplies y by \(R(z) = 1 + z\) with \(z = lam \cdot dt\), so the method is stable only inside the disc \(|1 + z| \le 1\) (see mathematicskit.integrators.stability).

Parameters:
Return type:

NDArray[double]

Returns:

ndarray of float, shape (dim,) – The state advanced by one step of size dt.

mathematicskit.integrators.implicit_euler_integrate(rhs, state0, t0, dt, n_steps, params, tol=1e-10, max_iter=50)[source]#

Integrate n_steps of backward Euler starting from state0.

First-order accurate but A-stable: for any dt > 0 a decaying mode stays decaying, so stiff systems can be stepped at a dt far beyond the explicit stability limit (compare rk4_integrate()).

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n_steps + 1,)) – Time at each step, starting at t0.

  • states (ndarray of float, shape (n_steps + 1, dim)) – State at each step, starting at state0.

mathematicskit.integrators.implicit_euler_step(rhs, state, t, dt, params, tol=1e-10, max_iter=50)[source]#

Single backward (implicit) Euler step, solved by Newton’s method.

Solves G(y) = y - state - dt * rhs(y, t + dt, params) = 0 with Newton iterations (I - dt * J) delta = -G(y), where J is a forward-difference Jacobian of rhs re-evaluated every iteration.

Parameters:
Return type:

NDArray[double]

Returns:

ndarray of float, shape (dim,) – The state at t + dt.

Raises:

RuntimeError – If Newton’s method does not converge within max_iter iterations.

mathematicskit.integrators.is_absolutely_stable(method, z)[source]#

Whether method is absolutely stable at each \(z = \lambda h\).

Parameters:

  • method (str) – One of STABILITY_METHODS: "euler", "implicit_euler", "rk4", or "adams_bashforth2" to "adams_bashforth4".

  • z (complex or array-like of complex)

Return type:

ndarray

Returns:

ndarray of bool, same shape as z – True inside the stability region (boundary included).

Examples

>>> is_absolutely_stable("euler", [-1.0, -2.5]).tolist()
[True, False]
>>> # AB2's real stability interval is [-1, 0], half of Euler's.
>>> is_absolutely_stable("adams_bashforth2", [-0.9, -1.1]).tolist()
[True, False]
>>> bool(is_absolutely_stable("implicit_euler", -1e6))
True
mathematicskit.integrators.leapfrog_integrate(force, pos0, vel0, t0, dt, n_steps, params)[source]#

Integrate n_steps of the symplectic leapfrog scheme.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n_steps + 1,)) – Time at each step, starting at t0.

  • positions (ndarray of float, shape (n_steps + 1, dim)) – Position at each step, starting at pos0.

  • velocities (ndarray of float, shape (n_steps + 1, dim)) – Velocity at each step, starting at vel0.

mathematicskit.integrators.leapfrog_step(force, pos, vel, t, dt, params)[source]#

Single velocity-Verlet (symplectic) step for pos'' = force(pos, t).

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • pos_new (ndarray of float, shape (dim,)) – Position advanced by one step of size dt.

  • vel_new (ndarray of float, shape (dim,)) – Velocity advanced by one step of size dt.

mathematicskit.integrators.njit(*args, **kws)[source]#

Legacy decorator that is equivalent to the preferred API: jit().

See documentation for jit function/decorator for full description.

Return type:

_JITWrapper

mathematicskit.integrators.rk4_integrate(rhs, state0, t0, dt, n_steps, params)[source]#

Integrate n_steps of RK4 starting from state0.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n_steps + 1,)) – Time at each step, starting at t0.

  • states (ndarray of float, shape (n_steps + 1, dim)) – State at each step, starting at state0.

mathematicskit.integrators.rk4_step(rhs, state, t, dt, params)[source]#

Single classical 4th-order Runge-Kutta step.

Parameters:
Return type:

NDArray[double]

Returns:

ndarray of float, shape (dim,) – The state advanced by one step of size dt.

mathematicskit.integrators.stability_function(method, z)[source]#

Stability function \(R(z)\) of a one-step method.

Parameters:
Return type:

ndarray

Returns:

ndarray of complex – \(R(z)\), the factor one step multiplies \(y\) by on \(y' = \lambda y\).

Examples

>>> import numpy as np
>>> complex(stability_function("euler", -0.5))
(0.5+0j)
>>> # RK4 reproduces e^z to fourth order.
>>> bool(abs(stability_function("rk4", -0.1) - np.exp(-0.1)) < 1e-7)
True
mathematicskit.integrators.stiff_integrate(rhs, state0, t0, t_end, params, method='Radau', rtol=1e-06, atol=1e-09, t_eval=None, jac=None)[source]#

Integrate a stiff system with scipy.integrate.solve_ivp().

Adapts the rhs(state, t, params) convention to scipy’s fun(t, y) and returns the same (times, states) layout as dopri5_integrate(), so the two are drop-in comparable.

Parameters:

  • rhs (Callable[[NDArray[double], float, NDArray[double]], NDArray[double]]) – Right-hand side rhs(state, t, params) -> ndarray; an @njit function or a plain Python one.

  • state0 (ArrayLike) – Initial state vector.

  • t0 (float) – Integration interval.

  • t_end (float) – Integration interval.

  • params (NDArray[double]) – Parameter vector passed through to rhs (and jac).

  • method (str) – Implicit Radau IIA (order 5), variable-order BDF (orders 1-5), or LSODA (switches automatically between Adams and BDF).

  • rtol (float) – Relative and absolute error tolerances.

  • atol (float) – Relative and absolute error tolerances.

  • t_eval (ArrayLike | None) – Times at which to report the solution; by default the solver’s own accepted steps are returned.

  • jac (Callable[[NDArray[double], float, NDArray[double]], NDArray[double]] | None) – Analytic Jacobian jac(state, t, params) -> ndarray (dim, dim); estimated by finite differences if omitted.

Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n,)) – Output times, from t0 to t_end.

  • states (ndarray of float, shape (n, dim)) – State at each output time.

Raises:

  • ValueError – If method is not one of the supported implicit methods.

  • RuntimeError – If the solver fails.

Examples

>>> import numpy as np
>>> def decay(state, t, params):
...     return -params[0] * state
>>> ts, ys = stiff_integrate(decay, [1.0], 0.0, 1.0, np.array([1e4]), method="BDF")
>>> len(ts) < 200, bool(abs(ys[-1, 0]) < 1e-8)
(True, True)
mathematicskit.integrators.velocity_verlet_integrate(force, pos0, vel0, t0, dt, n_steps, params)#

Integrate n_steps of the symplectic leapfrog scheme.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n_steps + 1,)) – Time at each step, starting at t0.

  • positions (ndarray of float, shape (n_steps + 1, dim)) – Position at each step, starting at pos0.

  • velocities (ndarray of float, shape (n_steps + 1, dim)) – Velocity at each step, starting at vel0.

mathematicskit.integrators.velocity_verlet_step(force, pos, vel, t, dt, params)#

Single velocity-Verlet (symplectic) step for pos'' = force(pos, t).

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • pos_new (ndarray of float, shape (dim,)) – Position advanced by one step of size dt.

  • vel_new (ndarray of float, shape (dim,)) – Velocity advanced by one step of size dt.

mathematicskit.integrators.yoshida4_integrate(force, pos0, vel0, t0, dt, n_steps, params)[source]#

Integrate n_steps of the 4th-order symplectic Yoshida scheme.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double], NDArray[double]]

Returns:

  • times (ndarray of float, shape (n_steps + 1,)) – Time at each step, starting at t0.

  • positions (ndarray of float, shape (n_steps + 1, dim)) – Position at each step, starting at pos0.

  • velocities (ndarray of float, shape (n_steps + 1, dim)) – Velocity at each step, starting at vel0.

mathematicskit.integrators.yoshida4_step(force, pos, vel, t, dt, params)[source]#

Single 4th-order symplectic (Yoshida) step for pos'' = force(pos, t).

Composes three leapfrog_step() sub-steps with Yoshida’s (1990) coefficients to raise the (still symplectic) accuracy from 2nd to 4th order, at roughly 3x the cost per step of plain leapfrog. Because it remains symplectic, it – like leapfrog – conserves invariants far better than RK4 over long integrations, now with much smaller local truncation error too.

Parameters:
Return type:

tuple[NDArray[double], NDArray[double]]

Returns:

  • pos_new (ndarray of float, shape (dim,)) – Position advanced by one step of size dt.

  • vel_new (ndarray of float, shape (dim,)) – Velocity advanced by one step of size dt.