Runge-Kutta Integrators

The rk module provides explicit Runge-Kutta integrators used throughout the project.

Base Classes

class hiten.algorithms.integrators.rk._RungeKuttaBase[source]

Bases: _Integrator

Provide shared functionality of explicit Runge-Kutta schemes.

The class stores a Butcher tableau and provides a single low level helper _rk_embedded_step() that advances one macro time step and, when a second set of weights is available, returns an error estimate suitable for adaptive step-size control.

Parameters:

name (str)

_A

Strictly lower triangular array of stage coefficients a_ij.

Type:

numpy.ndarray of shape (s, s)

_B_HIGH

Weights of the high order solution.

Type:

numpy.ndarray of shape (s,)

_B_LOW

Weights of the lower order solution, optional. When None no error estimate is produced and rk_embedded_step_jit_kernel() falls back to the high order result for both outputs.

Type:

numpy.ndarray or None

_C

Nodes c_i measured in units of the step size.

Type:

numpy.ndarray of shape (s,)

_p

Formal order of accuracy of the high order scheme.

Type:

int

Notes

The class is not intended to be used directly. Concrete subclasses define the specific coefficients and expose a public interface compliant with _Integrator.

property order: int

Return the formal order of accuracy of the method.

Returns:

The order of accuracy of the Runge-Kutta method.

Return type:

int

Fixed-Step Classes

class hiten.algorithms.integrators.rk._FixedStepRK[source]

Bases: _RungeKuttaBase

Implement an explicit fixed-step Runge-Kutta scheme.

Parameters:
  • name (str) – Human readable identifier of the scheme (e.g. "_RK4").

  • A (numpy.ndarray) – Butcher tableau as returned by *.

  • B (numpy.ndarray) – Butcher tableau as returned by *.

  • C (numpy.ndarray) – Butcher tableau as returned by *.

  • order (int) – Formal order of accuracy p of the method.

  • **options – Additional keyword options forwarded to the base _Integrator.

Notes

The step size is assumed to be constant and is inferred from the spacing of the t_vals array supplied to integrate().

integrate(system, y0, t_vals, *, event_fn=None, event_cfg=None, **kwargs)[source]

Integrate with a fixed-step Runge-Kutta method, with events.

Parameters:
  • system (_DynamicalSystemProtocol) – The dynamical system to integrate.

  • y0 (numpy.ndarray) – Initial state of shape (system.dim,). Units follow the provided system.

  • t_vals (numpy.ndarray) – Strictly monotonic time nodes of shape (N,) at which to evaluate the solution. Units follow the provided system.

  • event_fn (Callable[[float, numpy.ndarray], float], optional) – Scalar event function evaluated as g(t, y). A zero crossing may terminate integration or mark an event.

  • event_cfg (_EventConfig | None) – Configuration controlling event directionality, terminal behavior, and tolerances.

  • **kwargs – Additional integration options passed to the implementation.

Returns:

Integration results with times, states, and derivatives when available. Units follow the provided system.

Return type:

_Solution

class hiten.algorithms.integrators.rk._RK4[source]

Bases: _FixedStepRK

Implement the classical 4th-order Runge-Kutta method.

This is the standard 4th-order explicit Runge-Kutta method, also known as RK4 or the “classical” Runge-Kutta method. It uses 4 function evaluations per step and has order 4.

class hiten.algorithms.integrators.rk._RK6[source]

Bases: _FixedStepRK

Implement a 6th-order Runge-Kutta method.

A 6th-order explicit Runge-Kutta method that provides higher accuracy than RK4 at the cost of more function evaluations per step.

class hiten.algorithms.integrators.rk._RK8[source]

Bases: _FixedStepRK

Implement an 8th-order Runge-Kutta method.

An 8th-order explicit Runge-Kutta method that provides very high accuracy for applications requiring precise numerical integration.

Adaptive Classes

class hiten.algorithms.integrators.rk._AdaptiveStepRK[source]

Bases: _RungeKuttaBase

Implement an embedded adaptive Runge-Kutta integrator with PI controller.

The class provides common constants for PI step-size control; concrete methods (e.g. RK45, DOP853) implement the integration drivers.

Parameters:
  • name (str, default "AdaptiveRK") – Identifier passed to the _Integrator base class.

  • rtol (float, optional) – Relative and absolute error tolerances. Defaults are read from TOL.

  • atol (float, optional) – Relative and absolute error tolerances. Defaults are read from TOL.

  • max_step (float, optional) – Upper bound on the step size. infinity by default.

  • min_step (float or None, optional) – Lower bound on the step size. When None the value is derived from machine precision.

SAFETY, MIN_FACTOR, MAX_FACTOR

Magic constants used by the PI controller. They follow SciPy’s implementation and the recommendations by Hairer et al.

Type:

float

Raises:

RuntimeError – If the step size underflows while trying to satisfy the error tolerance.

Parameters:
SAFETY = 0.9
MIN_FACTOR = 0.2
MAX_FACTOR = 10.0
class hiten.algorithms.integrators.rk._RK45[source]

Bases: _AdaptiveStepRK

Implement the Dormand-Prince 5(4) adaptive Runge-Kutta method.

This is the Dormand-Prince 5th-order adaptive Runge-Kutta method with 4th-order error estimation. It provides a good balance between accuracy and computational efficiency for most applications.

integrate(system, y0, t_vals, *, event_fn=None, event_cfg=None, **kwargs)[source]

Integrate with RK45 and dense interpolation, with events.

Parameters:
  • system (_DynamicalSystemProtocol) – The dynamical system to integrate.

  • y0 (numpy.ndarray) – Initial state of shape (system.dim,). Units follow the provided system.

  • t_vals (numpy.ndarray) – Time nodes of shape (N,) at which to evaluate the solution. Units follow the provided system.

  • event_fn (Callable[[float, numpy.ndarray], float], optional) – Scalar event function evaluated as g(t, y).

  • event_cfg (_EventConfig | None) – Configuration controlling directionality, terminal behavior, and tolerances.

  • **kwargs – Additional integration options passed to the implementation.

Returns:

Integration results with times, states, and derivatives when available. Units follow the provided system.

Return type:

_Solution

class hiten.algorithms.integrators.rk._DOP853[source]

Bases: _AdaptiveStepRK

Implement the Dormand-Prince 8(5,3) adaptive Runge-Kutta method.

This is the Dormand-Prince 8th-order adaptive Runge-Kutta method with 5th and 3rd-order error estimation. It provides very high accuracy for applications requiring precise numerical integration.

integrate(system, y0, t_vals, *, event_fn=None, event_cfg=None, **kwargs)[source]

Integrate with DOP853 and dense interpolation, with events.

Parameters:
  • system (_DynamicalSystemProtocol) – The dynamical system to integrate.

  • y0 (numpy.ndarray) – Initial state of shape (system.dim,). Units follow the provided system.

  • t_vals (numpy.ndarray) – Time nodes of shape (N,) at which to evaluate the solution. Units follow the provided system.

  • event_fn (Callable[[float, numpy.ndarray], float], optional) – Scalar event function evaluated as g(t, y).

  • event_cfg (_EventConfig) – Configuration controlling directionality, terminal behavior, and tolerances.

  • **kwargs – Additional integration options passed to the implementation.

Returns:

Integration results with times, states, and derivatives when available. Units follow the provided system.

Return type:

_Solution

Factory Classes

class hiten.algorithms.integrators.rk.FixedRK[source]

Bases: object

Implement a factory class for creating fixed-step Runge-Kutta integrators.

This factory provides convenient access to fixed-step Runge-Kutta methods of different orders. The available orders are 4, 6, and 8.

Examples

>>> rk4 = FixedRK(order=4)
>>> rk6 = FixedRK(order=6)
>>> rk8 = FixedRK(order=8)
static __new__(cls, order=4, **opts)[source]

Create a fixed-step Runge-Kutta integrator of specified order.

Parameters:
  • order (int, default 4) – Order of the Runge-Kutta method. Must be 4, 6, or 8.

  • **opts – Additional options passed to the integrator constructor.

Returns:

A fixed-step Runge-Kutta integrator instance.

Return type:

_FixedStepRK

Raises:

ValueError – If the specified order is not supported.

class hiten.algorithms.integrators.rk.AdaptiveRK[source]

Bases: object

Implement a factory class for creating adaptive step-size Runge-Kutta integrators.

This factory provides convenient access to adaptive step-size Runge-Kutta methods. The available orders are 5 (Dormand-Prince 5(4)) and 8 (Dormand-Prince 8(5,3)).

Examples

>>> rk45 = AdaptiveRK(order=5)
>>> dop853 = AdaptiveRK(order=8)
static __new__(cls, order=5, **opts)[source]

Create an adaptive step-size Runge-Kutta integrator of specified order.

Parameters:
  • order (int, default 5) – Order of the Runge-Kutta method. Must be 5 or 8.

  • **opts – Additional options passed to the integrator constructor.

Returns:

An adaptive step-size Runge-Kutta integrator instance.

Return type:

_AdaptiveStepRK

Raises:

ValueError – If the specified order is not supported.

class hiten.algorithms.integrators.rk.RungeKutta[source]

Bases: object

Implement a factory class for creating Runge-Kutta integrators.

This factory provides convenient access to Runge-Kutta integrators of different orders.

static __new__(cls, order=4, **opts)[source]

Create a Runge-Kutta integrator of specified order.

Parameters:
  • order (int, default 4) – Order of the Runge-Kutta method. Must be 4, 45, 6, 8, or 853.

  • **opts – Additional options passed to the integrator constructor.

Returns:

Kernel Functions

hiten.algorithms.integrators.rk.rk_embedded_step_jit_kernel()[source]

Perform a single step of the RK method (no f closure).

Parameters:
  • f (callable) – The right-hand side of the differential equation.

  • t (float) – The current time.

  • y (numpy.ndarray) – The current state.

  • h (float) – The step size.

  • A (numpy.ndarray) – The Butcher tableau.

  • B_HIGH (numpy.ndarray) – The high order weights.

  • B_LOW (numpy.ndarray) – The low order weights.

  • C (numpy.ndarray) – The nodes.

  • has_b_low (bool) – Whether to use the low order weights.

Returns:

  • numpy.ndarray – The high order solution.

  • numpy.ndarray – The low order solution.

  • numpy.ndarray – The error vector.

hiten.algorithms.integrators.rk.rk_embedded_step_ham_jit_kernel()[source]

Hamiltonian variant of a single RK embedded step (no f closure).

Parameters:
Returns:

  • numpy.ndarray – The high order solution.

  • numpy.ndarray – The low order solution.

  • numpy.ndarray – The error vector.

hiten.algorithms.integrators.rk.rk45_step_jit_kernel()[source]

Perform a single step of the RK45 method.

Parameters:
Returns:

  • numpy.ndarray – The high order solution.

  • numpy.ndarray – The low order solution.

  • numpy.ndarray – The error vector.

hiten.algorithms.integrators.rk.rk45_step_ham_jit_kernel()[source]

Hamiltonian variant of a single RK45 step (no f closure).

Parameters:
Returns:

  • numpy.ndarray – The high order solution.

  • numpy.ndarray – The low order solution.

  • numpy.ndarray – The error vector.

  • numpy.ndarray – The intermediate stages.

hiten.algorithms.integrators.rk.dop853_step_jit_kernel()[source]

Perform a single step of the DOP853 method.

Parameters:
Returns:

  • numpy.ndarray – The high order solution.

  • numpy.ndarray – The low order solution.

  • numpy.ndarray – The error vector.

hiten.algorithms.integrators.rk.dop853_step_ham_jit_kernel()[source]

Hamiltonian variant of a single DOP853 step (no f closure).

Parameters:
Returns:

  • numpy.ndarray – The high order solution.

  • numpy.ndarray – The low order solution.

  • numpy.ndarray – The error vector.

  • numpy.ndarray – The intermediate stages.

Dense Output Functions

hiten.algorithms.integrators.rk._hermite_eval_dense()[source]

Evaluate cubic Hermite interpolant between (t0,y0,f0) and (t1,y1,f1).

Uses standard basis: H00=2x^3-3x^2+1, H10=x^3-2x^2+x, H01=-2x^3+3x^2, H11=x^3-x^2.

Parameters:
Returns:

The interpolated state.

Return type:

numpy.ndarray

hiten.algorithms.integrators.rk._rk45_build_Q_cache()[source]

Build the Q cache for the RK45 method.

Parameters:
Returns:

The Q cache.

Return type:

numpy.ndarray

hiten.algorithms.integrators.rk._rk45_eval_dense()[source]

Evaluate the dense interpolator at time x.

Parameters:
Returns:

The interpolated state.

Return type:

numpy.ndarray

hiten.algorithms.integrators.rk._dop853_build_dense_cache()[source]

Build the dense cache for the DOP853 method.

Parameters:
  • f (callable) – The right-hand side of the differential equation.

  • t_old (float) – The initial time.

  • y_old (numpy.ndarray) – The initial state.

  • f_old (numpy.ndarray) – The initial derivative.

  • y_new (numpy.ndarray) – The final state.

  • f_new (numpy.ndarray) – The final derivative.

  • hseg (float) – The step size.

  • Kseg (numpy.ndarray) – The intermediate stages.

  • A_full (numpy.ndarray) – The full Butcher tableau.

  • C_full (numpy.ndarray) – The full Butcher tableau.

  • D (numpy.ndarray) – The full Butcher tableau.

  • n_stages_extended (int) – The number of stages.

  • interpolator_power (int) – The power of the interpolator.

Returns:

The dense cache.

Return type:

numpy.ndarray

Notes

This function implements the dense cache for the DOP853 method. It uses the dense cache to evaluate the interpolator at the time x. It returns the interpolated state. The interpolated state is returned in the interval [y0, y1]. The time of the event is returned in the interval [t0, t1].

hiten.algorithms.integrators.rk._dop853_eval_dense()[source]

Evaluate the dense interpolator at time x.

Parameters:
  • y_old (numpy.ndarray) – The initial state.

  • F_cache (numpy.ndarray) – The dense cache.

  • interpolator_power (int) – The power of the interpolator.

  • x (float) – The time to evaluate the interpolator at.

Returns:

The interpolated state.

Return type:

numpy.ndarray

Notes

This function implements the dense interpolator. It uses the dense cache to evaluate the interpolator at the time x. It returns the interpolated state. The interpolated state is returned in the interval [y0, y1]. The time of the event is returned in the interval [t0, t1].

Event Functions

hiten.algorithms.integrators.rk._hermite_refine_in_step()[source]

Refine the integration step using bisection on x in [0, 1].

Parameters:
  • event_fn (callable) – The event function.

  • t0 (float) – The initial time.

  • y0 (numpy.ndarray) – The initial state.

  • f0 (numpy.ndarray) – The initial derivative.

  • t1 (float) – The final time.

  • y1 (numpy.ndarray) – The final state.

  • f1 (numpy.ndarray) – The final derivative.

  • h (float) – The step size.

  • direction (int) – The direction of the event.

  • xtol (float) – The tolerance.

Returns:

  • float – The time of the event.

  • numpy.ndarray – The state at the event.

Notes

This function implements the bisection method to find the time of the event. It uses the dense interpolator to evaluate the event function at the time of the event. It returns the time of the event and the state at the event. The time of the event is returned in the interval [t0, t1].

hiten.algorithms.integrators.rk._rk45_refine_in_step()[source]

Refine the integration step using bisection on x in [0, 1].

Parameters:
  • event_fn (callable) – The event function.

  • t0 (float) – The initial time.

  • y0 (numpy.ndarray) – The initial state.

  • t1 (float) – The final time.

  • y1 (numpy.ndarray) – The final state.

  • h (float) – The step size.

  • Kseg (numpy.ndarray) – The intermediate stages.

  • P (numpy.ndarray) – The nodes.

  • direction (int) – The direction of the event.

  • xtol (float) – The tolerance.

Returns:

  • float – The time of the event.

  • numpy.ndarray – The state at the event.

hiten.algorithms.integrators.rk._dop853_refine_in_step()[source]

Refine the integration step using bisection on x in [0, 1].

Parameters:
  • f (callable) – The right-hand side of the differential equation.

  • event_fn (callable) – The event function.

  • t0 (float) – The initial time.

  • y0 (numpy.ndarray) – The initial state.

  • f0 (numpy.ndarray) – The initial derivative.

  • t1 (float) – The final time.

  • y1 (numpy.ndarray) – The final state.

  • f1 (numpy.ndarray) – The final derivative.

  • h (float) – The step size.

  • Kseg (numpy.ndarray) – The intermediate stages.

  • A_full (numpy.ndarray) – The full Butcher tableau.

  • C_full (numpy.ndarray) – The full Butcher tableau.

  • D (numpy.ndarray) – The full Butcher tableau.

  • n_stages_extended (int) – The number of stages.

  • interpolator_power (int) – The power of the interpolator.

  • direction (int) – The direction of the event.

  • xtol (float) – The tolerance.

Returns:

  • float – The time of the event.

  • numpy.ndarray – The state at the event.

Notes

This function implements the bisection method to find the time of the event. It uses the dense cache to evaluate the event function at the endpoints and the intermediate points. It then uses the bisection method to find the time of the event. It returns the time of the event and the state at the event. The time of the event is returned in the interval [t0, t1]. The state at the event is returned in the interval [y0, y1].

hiten.algorithms.integrators.rk._dop853_refine_in_step_ham()[source]

Hamiltonian variant of DOP853 in-step root refinement.

Mirrors _dop853_refine_in_step but builds missing stages via _hamiltonian_rhs.