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 andrk_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,)
Notes
The class is not intended to be used directly. Concrete subclasses define the specific coefficients and expose a public interface compliant with
_Integrator
.
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 providedsystem
.t_vals (numpy.ndarray) – Strictly monotonic time nodes of shape
(N,)
at which to evaluate the solution. Units follow the providedsystem
.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:
- 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:
- 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 providedsystem
.t_vals (numpy.ndarray) – Time nodes of shape
(N,)
at which to evaluate the solution. Units follow the providedsystem
.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:
- 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 providedsystem
.t_vals (numpy.ndarray) – Time nodes of shape
(N,)
at which to evaluate the solution. Units follow the providedsystem
.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:
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:
- 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:
- 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:
FixedRK
orAdaptiveRK
– A Runge-Kutta integrator instance.
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:
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.
jac_H (numpy.ndarray) – The Jacobian of the Hamiltonian.
clmo_H (numpy.ndarray) – The coefficient-layout mapping objects for the Hamiltonian.
n_dof (int) – The number of degrees of freedom.
- 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:
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.
C (numpy.ndarray) – The nodes.
E (numpy.ndarray) – The error weights.
- 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:
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.
C (numpy.ndarray) – The nodes.
E (numpy.ndarray) – The error weights.
jac_H (numpy.ndarray) – The Jacobian of the Hamiltonian.
clmo_H (numpy.ndarray) – The coefficient-layout mapping objects for the Hamiltonian.
n_dof (int) – The number of degrees of freedom.
- 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:
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.
C (numpy.ndarray) – The nodes.
E5 (numpy.ndarray) – The error weights.
E3 (numpy.ndarray) – The error weights.
- 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:
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.
C (numpy.ndarray) – The nodes.
E5 (numpy.ndarray) – The error weights.
E3 (numpy.ndarray) – The error weights.
jac_H (numpy.ndarray, shape (n_dof, n_dof)) – The Jacobian of the Hamiltonian.
clmo_H (numpy.ndarray) – The coefficient-layout mapping objects for the Hamiltonian.
n_dof (int, shape (n_dof,)) – The number of degrees of freedom.
- 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:
y0 (numpy.ndarray) – The initial state.
f0 (numpy.ndarray) – The initial derivative.
y1 (numpy.ndarray) – The final state.
f1 (numpy.ndarray) – The final derivative.
x (float) – The time to evaluate the interpolator at.
h (float) – The step size.
- Returns:
The interpolated state.
- Return type:
- hiten.algorithms.integrators.rk._rk45_build_Q_cache()[source]
Build the Q cache for the RK45 method.
- Parameters:
Kseg (numpy.ndarray) – The intermediate stages.
P (numpy.ndarray) – The nodes.
dim (int) – The dimension of the state.
- Returns:
The Q cache.
- Return type:
- hiten.algorithms.integrators.rk._rk45_eval_dense()[source]
Evaluate the dense interpolator at time x.
- Parameters:
y_old (numpy.ndarray) – The initial state.
Q_cache (numpy.ndarray) – The Q cache.
P (numpy.ndarray) – The nodes.
x (float) – The time to evaluate the interpolator at.
hseg (float) – The step size.
- Returns:
The interpolated state.
- Return type:
- 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:
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:
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].