rkstiff.etd35

Adaptive-Step Fifth Order (Third Order Embedding) Exponential Time-Differencing Integrator

Exponential Time-Differencing Runge-Kutta Integrator (ETD(3,5))

Implements the ETD(3,5) exponential time-differencing scheme with embedded third-order error estimation and adaptive step control as described in:

Whalen, P., Brio, M., & Moloney, J. V. (2015). Exponential time-differencing with embedded Runge-Kutta adaptive step control. Journal of Computational Physics, 280, 579-601. doi:[10.1016/j.jcp.2014.09.024](https://doi.org/10.1016/j.jcp.2014.09.024)

—

Mathematical Formulation

This solver integrates stiff systems of ODEs or semidiscretized PDEs of the form

\[\frac{\partial \mathbf{U}}{\partial t} = \mathcal{L}\mathbf{U} + \mathcal{N}(\mathbf{U}),\]

where \(\mathcal{L}\) is the linear (possibly stiff) operator and \(\mathcal{N}\) is a nonlinear term evaluated explicitly.

The ETD(3,5) scheme uses exponential Runge–Kutta stages built from the ψ-functions:

\[\psi_r(z) = r \int_0^1 e^{(1-\theta)z}\,\theta^{r-1}\,d\theta, \quad r = 1,2,3,\dots\]

which serve as scaled exponential integrators analogous to the \(\phi_r\) functions in classical ETD formulations (\(\psi_r = r!\,\phi_r\)).

—

Implementation Overview

The module provides specialized strategies for different forms of the linear operator \(\mathcal{L}\):

  • _Etd35Diagonal — diagonal systems (elementwise exponentials)

  • _Etd35Diagonalized — eigen-decomposed systems

  • _Etd35NonDiagonal — full matrix exponential evaluation

  • ETD35 — high-level adaptive solver interface

Adaptive step control follows the embedded-order algorithm of Whalen et al. (2015), balancing efficiency and accuracy through local error estimates derived from the third-order embedding.

Classes

ETD35(lin_op, nl_func[, config, etd_config, ...])

Adaptive fifth-order Exponential Time-Differencing Runge–Kutta solver (ETD(3,5)).

_Etd35Diagonal(lin_op, nl_func, etd_config)

ETD(3, 5) diagonal formulation.

_Etd35Diagonalized(lin_op, nl_func, etd_config)

ETD35 solver for non-diagonal systems via eigenvector diagonalization.

_Etd35NonDiagonal(lin_op, nl_func, etd_config)

ETD35 solver for full (non-diagonal, non-diagonalizable) linear operators.
class rkstiff.etd35.ETD35(lin_op, nl_func, config=SolverConfig(), etd_config=ETDConfig(), diagonalize=False, loglevel='WARNING')[source]

Bases: ETDAS

Adaptive fifth-order Exponential Time-Differencing Runge–Kutta solver (ETD(3,5)).

Solves stiff systems of the form

\[\frac{\partial \mathbf{U}}{\partial t} = \mathcal{L}\mathbf{U} + \mathcal{N}(\mathbf{U}),\]

where \(\mathcal{L}\) is a linear operator and \(\mathcal{N}\) is a nonlinear term.

This implementation follows the ETD(3,5) algorithm developed by Whalen, Brio, and Moloney (2015), which embeds a third-order ETD scheme for adaptive time-step control within a fifth-order integrator.

The ETD(3,5) method advances the solution through exponential Runge-Kutta stages defined by the ψ-functions:

\[\psi_r(z) = r \int_0^1 e^{(1-\theta)z}\,\theta^{r-1}\,d\theta, \quad r = 1,2,3,\dots\]

These appear in the Runge-Kutta coefficients and can be related to the more common \(\phi\)-functions via \(\psi_r(z) = r!\,\phi_r(z)\).

The embedded third-order estimate is used to adapt the time step according to

\[h_{\text{new}} = h_{\text{old}}\,\nu \left( \frac{\varepsilon}{\mathrm{err}} \right)^{1/(q+1)}, \qquad q=4,\]

where \(\varepsilon\) is the tolerance, \(\nu\) is a safety factor, and \(\mathrm{err}\) is the estimated local truncation error.

Supports

  • Diagonal \(\mathcal{L}\) (elementwise exponentials)

  • Diagonalizable \(\mathcal{L}\) (eigenbasis integration)

  • Full \(\mathcal{L}\) (matrix exponentials via contour integration)

References

Whalen, P., Brio, M., & Moloney, J. V. (2015). Exponential time-differencing with embedded Runge–Kutta adaptive step control. Journal of Computational Physics, 280, 579–601.

__init__(lin_op, nl_func, config=SolverConfig(), etd_config=ETDConfig(), diagonalize=False, loglevel='WARNING')[source]

Initialize ETD35 adaptive solver.

Parameters:

  • lin_op (np.ndarray) – Linear operator L in the system. May be a 1D array (diagonal system) or a 2D square matrix (non-diagonal system).

  • nl_func (Callable[[np.ndarray], np.ndarray]) – Nonlinear function N(U).

  • config (SolverConfig, optional) – General solver configuration controlling adaptivity thresholds, safety factors, and other integration parameters.

  • etd_config (ETDConfig, optional) – Configuration for ETD-specific parameters, such as contour integration settings and spectral radius estimation.

  • diagonalize (bool, optional) – If True, attempts eigenvalue decomposition to transform system into diagonal form before solving.

  • loglevel (Union[str, int], optional) – Logging level.

MAX_LOOPS = 50

Maximum retry attempts per adaptive step

MAX_S = 4.0

Maximum allowed step size increase factor

MIN_S = 0.25

Minimum allowed step size reduction factor

exception MaxLoopsExceeded

Bases: SolverError

Raised when too many attempts are made to find a valid adaptive step.

exception MinimumStepReached

Bases: SolverError

Raised when the adaptive step size falls below the minimum allowed value.

exception SolverError

Bases: RuntimeError

Base exception for solver-related runtime errors.

evolve(u, t0, tf, h_init=None, store_data=True, store_freq=1)

Integrate the system from \(t_0\) to \(t_f\) using adaptive time steps.

Repeatedly applies step() to propagate the solution forward while dynamically adjusting the time step size based on local error estimates.

Parameters:

  • u (np.ndarray) – Initial solution vector at \(t_0\).

  • t0 (float) – Initial time.

  • tf (float) – Final time.

  • h_init (float, optional) – Initial step size. Defaults to (tf - t0) / 100 if not provided.

  • store_data (bool, default=True) – Whether to store intermediate results in t and u.

  • store_freq (int, default=1) – Frequency of data storage; store every store_freq accepted steps.

Returns:

Final solution at \(t = t_f\).

Return type:

np.ndarray

Notes

  • The evolution proceeds until \(t \geq t_f\), automatically adjusting step sizes as needed.

  • Stored data is accessible via t and u.

Example

>>> solver = ETD35(lin_op, nl_func)
>>> u_final = solver.evolve(u0, t0=0.0, tf=10.0)
>>> solver.t[-1], np.linalg.norm(solver.u[-1])
(10.0, 0.0134)
reset()

Reset solver and clear adaptive-step state.

Return type:

None

set_loglevel(loglevel)

Adjust the solver’s logging verbosity at runtime.

Parameters:

loglevel (str or int) – New logging level. Accepts standard string levels or numeric constants from logging.

Return type:

None

Examples

>>> solver.set_loglevel("INFO")
>>> solver.set_loglevel(logging.DEBUG)
property solver_type: SolverType

Return the solver type for adaptive-step solvers.

Returns:

Always returns SolverType.ADAPTIVE_STEP.

Return type:

SolverType

Examples

>>> from rkstiff.etd35 import ETD35
>>> solver = ETD35(lin_op, nl_func)
>>> solver.solver_type == SolverType.ADAPTIVE_STEP
True
step(u, h_suggest)

Perform one adaptive integration step.

Attempts to advance the solution by time h_suggest and adjusts the step size automatically based on local error estimates.

Parameters:

  • u (np.ndarray) – Current state vector.

  • h_suggest (float) – Suggested time step size.

Return type:

Tuple[ndarray, float, float]

Returns:

  • unew (np.ndarray) – Updated solution vector after one accepted step.

  • h (float) – Actual step size used.

  • h_suggest (float) – Suggested step size for the next iteration.

Raises:

Notes

The algorithm follows this pattern:

  1. Try the proposed step.

  2. Estimate local error and compute scaling factor s.

  3. If s < 1 → reject step and reduce h.

  4. If s ≥ 1 → accept step and update h_suggest for next step.

Warning

Very small or divergent s values may indicate instability or excessively tight tolerances.

etd_config: ETDConfig