rkstiff.etd

Exponential Time-Differencing (ETD) Runge-Kutta Methods

Implements Exponential Time-Differencing (ETD) Runge-Kutta methods for solving stiff partial differential equations (PDEs) of the form

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

where \(\mathcal{L}\) is a stiff linear operator and \(\mathcal{N}(\mathbf{U})\) is a nonlinear, typically non-stiff term.

These methods use analytic integration of the linear part via the matrix exponential and handle the nonlinear part with Runge-Kutta-like stages built around the psi-functions:

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

Contents

ETDConfig - configuration for contour integration and cutoff settings.
psi1(), psi2(), psi3() - exponential helper functions.
ETDAS - adaptive-step ETD solver.
ETDCS - constant-step ETD solver.

Notes

ETD methods are particularly efficient when \(\mathcal{N}(\mathbf{U})\) evolves on a slower time scale than the linear component \(\mathcal{L}\mathbf{U}\).

Functions

`psi1`(z)	Compute \(\psi_1(z) = \frac{e^z - 1}{z}\) element-wise.
`psi2`(z)	Compute \(\psi_2(z) = 2*\frac{e^z - 1 - z}{z^2}\) element-wise.
`psi3`(z)	Compute \(\psi_3(z) = 6*\frac{e^z - 1 - z - \frac{z^2}{2}}{z^3}\) element-wise.

Classes

`ETDAS`(lin_op, nl_func[, config, etd_config, ...])	Adaptive-step Exponential Time-Differencing (ETD) Runge–Kutta solver.
`ETDCS`(lin_op, nl_func[, etd_config, loglevel])	Constant-step Exponential Time-Differencing (ETD) Runge–Kutta solver.
`ETDConfig`([modecutoff, contour_points, ...])	Configuration parameters for Exponential Time-Differencing (ETD) solvers.

class rkstiff.etd.ETDConfig(modecutoff=0.01, contour_points=32, contour_radius=1.0)[source]

Bases: object

Configuration parameters for Exponential Time-Differencing (ETD) solvers.

These parameters control the numerical evaluation of the psi-functions and the stability of contour-integral approximations.

Parameters:

modecutoff (float, optional) – Threshold for switching between direct evaluation and contour integration. For \(|z| < \text{modecutoff}\), \(\psi_k(z)\) is computed via Contour integration; otherwise, direct evaluation is used.
contour_points (int, optional) – Number of quadrature nodes used for contour integration.
contour_radius (float, optional) – Radius \(R\) of the circular contour in the complex plane used for the Cauchy integral evaluation of \(\psi_k(z)\).

Notes

For small \(z\), the ETD psi-functions are numerically unstable via direct evaluation and a contour integral representation is preferred.

\[\psi_k(z) = \frac{1}{2\pi i} \oint_\Gamma \frac{e^\xi}{(\xi - z)\xi^{k-1}}\,d\xi, \quad \Gamma : |\xi| = R.\]

For larger \(z\) (not near zero) direct evaluation of the ETD psi-functions works well.

__init__(modecutoff=0.01, contour_points=32, contour_radius=1.0)[source]

property modecutoff: float

Numerical inaccuracy cutoff for psi functions.

Must be between 0.0 and 1.0.

property contour_points: int

Number of points for the contour integral.

Must be an integer greater than 1.

property contour_radius: float

Radius of the circular contour integral.

Must be greater than 0.

rkstiff.etd.psi1(z)[source]

Compute \(\psi_1(z) = \frac{e^z - 1}{z}\) element-wise.

Parameters:: z (np.ndarray) – Real or complex array.
Returns:: Array of the same shape as z.
Return type:: np.ndarray

rkstiff.etd.psi2(z)[source]

Compute \(\psi_2(z) = 2*\frac{e^z - 1 - z}{z^2}\) element-wise.

Parameters:: z (np.ndarray) – Real or complex array.
Returns:: Array of the same shape as z.
Return type:: np.ndarray

rkstiff.etd.psi3(z)[source]

Compute \(\psi_3(z) = 6*\frac{e^z - 1 - z - \frac{z^2}{2}}{z^3}\) element-wise.

Parameters:: z (np.ndarray) – Real or complex array.
Returns:: Array of the same shape as z.
Return type:: np.ndarray

class rkstiff.etd.ETDAS(lin_op, nl_func, config=SolverConfig(), etd_config=ETDConfig(), loglevel='WARNING')[source]

Bases: BaseSolverAS

Adaptive-step Exponential Time-Differencing (ETD) Runge–Kutta solver.

Integrates stiff PDE’s 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 nonlinear.

This variant adjusts the time step \(h\) dynamically based on embedded error estimates, maintaining both accuracy and efficiency.

Parameters:

lin_op (np.ndarray) – Linear operator \(\mathcal{L}\).
nl_func (Callable[[np.ndarray], np.ndarray]) – Nonlinear function \(\mathcal{N}(\mathbf{U})\).
config (SolverConfig, optional) – General solver configuration controlling adaptivity.
etd_config (ETDConfig, optional) – ETD-specific contour and cutoff settings.
loglevel (str or int, optional) – Logging verbosity level.

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

Initialize an adaptive-step solver with validated configuration.

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:

MaxLoopsExceeded – If too many attempts are made to find an acceptable step size.
MinimumStepReached – If the step size drops below the configured minimum minh.

Notes

The algorithm follows this pattern:

Try the proposed step.
Estimate local error and compute scaling factor s.
If s < 1 → reject step and reduce h.
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.

class rkstiff.etd.ETDCS(lin_op, nl_func, etd_config=ETDConfig(), loglevel='WARNING')[source]

Bases: BaseSolverCS

Constant-step Exponential Time-Differencing (ETD) Runge–Kutta solver.

Integrates stiff PDEs of the form

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

using a fixed step size \(h\).

Parameters:

lin_op (np.ndarray) – Linear operator \(\mathcal{L}\).
nl_func (Callable[[np.ndarray], np.ndarray]) – Nonlinear function \(\mathcal{N}(\mathbf{U})\).
etd_config (ETDConfig, optional) – ETD configuration settings.
loglevel (str or int, optional) – Logging level.

Notes

Constant-step ETD methods are efficient when the stiffness and smoothness of the system are known in advance. Each time step advances via

\[\mathbf{u}_{n+1} = e^{h\mathcal{L}}\mathbf{u}_n + h\sum_i b_i\,\mathcal{N}(\mathbf{k}_i),\]

where \(\mathbf{k}_i\) are stage vectors computed using matrix–function coefficients \(\psi_k(h\mathcal{L})\).

evolve(u, t0, tf, h, store_data=True, store_freq=1)

Integrate the system from \(t_0\) to \(t_f\) using fixed step size.

Repeatedly applies step() with constant \(h\) until the final time is reached.

Parameters:

u (np.ndarray) – Initial solution vector at \(t_0\).
t0 (float) – Initial time.
tf (float) – Final time (integration stops when t ≥ tf).
h (float) – Constant step size \(\Delta t\).
store_data (bool, default=True) – Whether to store intermediate results in t and u.
store_freq (int, default=1) – Frequency of storing data; every store_freq steps.

Returns:

Final solution vector at \(t_f\).

Return type:

np.ndarray

Raises:

ValueError – If h exceeds the total time span tf - t0.

Notes

The time grid is uniformly spaced with spacing \(h\).
Stored data can be accessed through t and u.

Example

>>> solver = ETD4(lin_op, nl_func)
>>> u_final = solver.evolve(u0, t0=0.0, tf=10.0, h=0.05)
>>> len(solver.t)
200

reset()

Reset solver state and clear stored time/solution data.

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 constant-step solvers.

Returns:: Always returns SolverType.CONSTANT_STEP.
Return type:: SolverType

Examples

>>> from rkstiff.if4 import IF4
>>> solver = IF4(lin_op, nl_func)
>>> solver.solver_type == SolverType.CONSTANT_STEP
True

step(u, h)

Perform a single constant-step propagation.

Parameters:

u (np.ndarray) – Current solution vector.
h (float) – Constant step size (must be non-negative).

Returns:

Updated solution after one full time step.

Return type:

np.ndarray

Notes

This method simply wraps _update_stages() and performs minimal validation and logging.

__init__(lin_op, nl_func, etd_config=ETDConfig(), loglevel='WARNING')[source]

Initialize a constant-step solver and validate inputs.