Solvers API Reference

ODE, PDE, and SDE solvers.

ODE Solvers

Fractional Ordinary Differential Equation Solvers

This module provides comprehensive solvers for fractional ODEs including various numerical methods, adaptive step size control, and error estimation.

Performance Note: - Uses FFT-based convolution for O(N log N) history summation instead of O(N²) - Intelligent backend selection for optimal performance (v2.1.0)

hpfracc.solvers.ode_solvers._get_gamma_function()

Get gamma function through adapter system.

hpfracc.solvers.ode_solvers._get_intelligent_selector()

Get intelligent backend selector instance.

hpfracc.solvers.ode_solvers._select_fft_backend(data_size)

Select optimal FFT backend based on data size.

Parameters:

data_size (int) – Number of elements to process

Returns:

“numpy”, “jax”, or “scipy”

Return type:

Backend name

hpfracc.solvers.ode_solvers._fft_convolution(coeffs, values, axis=0)

Fast convolution using FFT for O(N log N) performance with intelligent backend selection.

This replaces the O(N²) direct summation loop for history terms in fractional ODEs. Now uses intelligent backend selection to choose optimal FFT implementation.

Parameters:

• coeffs (ndarray) – Coefficient array (1D)

• values (ndarray) – Value array (can be 1D or 2D with axis specifying which dimension to convolve)

• axis (int) – Axis along which to perform convolution (default 0)

Returns:

Convolution result (same shape as values)

Return type:

ndarray

Mathematical basis:

conv(C, Y) = ifft(fft(C) * fft(Y))

Performance:

• Direct summation: O(N²)

• FFT-based: O(N log N)

• Backend selection: < 0.001 ms overhead

hpfracc.solvers.ode_solvers._fast_history_sum(coeffs, f_hist, reverse=True, verbose=False)

Compute weighted sum of history using FFT-based convolution.

This is the key optimization for fractional ODE solvers, replacing:

sum_{j=0}^{n} coeffs[n-j] * f_hist[j]

with an O(N log N) FFT-based approach instead of O(N²) direct summation.

Parameters:

• coeffs (ndarray) – Coefficient array of length N

• f_hist (ndarray) – History array of shape (N,) or (N, m) where m is state dimension

• reverse (bool) – If True, apply coefficients in reverse order (typical for convolution)

• verbose (bool)

Returns:

Weighted sum result (scalar or array of length m)

Return type:

ndarray

class hpfracc.solvers.ode_solvers.FixedStepODESolver(derivative_type='caputo', method='predictor_corrector', adaptive=True, tol=1e-06, max_iter=1000, *, fractional_order=None, rtol=None, atol=None, order=1, min_h=None, max_h=None, min_step=None, max_step=None)

Bases: object

Base class for fixed-step fractional ODE solvers.

Provides common functionality for solving fractional ordinary differential equations of the form:

D^α y(t) = f(t, y(t))

where D^α is a fractional derivative operator.

Parameters:

• derivative_type (str)

• method (str)

• adaptive (bool)

• tol (float)

• max_iter (int)

• fractional_order (float | FractionalOrder | None)

• rtol (float | None)

• atol (float | None)

• order (int)

• min_h (float | None)

• max_h (float | None)

• min_step (float | None)

• max_step (float | None)

__init__(derivative_type='caputo', method='predictor_corrector', adaptive=True, tol=1e-06, max_iter=1000, *, fractional_order=None, rtol=None, atol=None, order=1, min_h=None, max_h=None, min_step=None, max_step=None)

Initialize fractional ODE solver.

Parameters:

• derivative_type (str) – Type of fractional derivative (“caputo”, “riemann_liouville”, “grunwald_letnikov”)

• method (str) – Numerical method (“predictor_corrector”, “adams_bashforth”, “runge_kutta”)

• adaptive (bool) – Use adaptive step size control

• tol (float) – Tolerance for convergence

• max_iter (int) – Maximum number of iterations

• fractional_order (float | FractionalOrder | None)

• rtol (float | None)

• atol (float | None)

• order (int)

• min_h (float | None)

• max_h (float | None)

• min_step (float | None)

• max_step (float | None)

_validate_alpha(alpha)

Validate the fractional order.

Parameters:

alpha (float | FractionalOrder)

solve(f, t_span, y0, alpha, h=None, **kwargs)

Solve fractional ODE.

Parameters:

• f (Callable) – Right-hand side function f(t, y)

• t_span (Tuple[float, float]) – Time interval (t0, tf)

• y0 (float | ndarray) – Initial condition(s)

• alpha (float | FractionalOrder) – Fractional order

• h (float | None) – Step size (None for adaptive)

• **kwargs – Additional solver parameters

Returns:

Tuple of (t_values, y_values)

Return type:

Tuple[ndarray, ndarray]

_solve_predictor_corrector(f, t0, tf, y0, alpha, h, **kwargs)

Solve using Adams-Bashforth-Moulton predictor-corrector for Caputo fractional ODEs.

Based on the Volterra integral formulation: y(t) = y_0 + (1/Γ(α)) ∫_0^t (t-τ)^(α-1) f(τ, y(τ)) dτ

Reference: Diethelm, K. (2010). The Analysis of Fractional Differential Equations.

_solve_adams_bashforth(f, t0, tf, y0, alpha, h, **kwargs)

Solve using Adams-Bashforth method.

Parameters:

• f (Callable) – Right-hand side function

• t0 (float) – Initial time

• tf (float) – Final time

• y0 (float | ndarray) – Initial condition

• alpha (float | FractionalOrder) – Fractional order

• h (float) – Step size

• **kwargs – Additional parameters

Returns:

Tuple of (t_values, y_values)

Return type:

Tuple[ndarray, ndarray]

_solve_runge_kutta(f, t0, tf, y0, alpha, h, **kwargs)

Solve using fractional Runge-Kutta method.

Parameters:

• f (Callable) – Right-hand side function

• t0 (float) – Initial time

• tf (float) – Final time

• y0 (float | ndarray) – Initial condition

• alpha (float | FractionalOrder) – Fractional order

• h (float) – Step size

• **kwargs – Additional parameters

Returns:

Tuple of (t_values, y_values)

Return type:

Tuple[ndarray, ndarray]

_solve_euler(f, t0, tf, y0, alpha, h, **kwargs)

Solve using fractional Euler method.

Parameters:

• f (Callable) – Right-hand side function

• t0 (float) – Initial time

• tf (float) – Final time

• y0 (float | ndarray) – Initial condition

• alpha (float | FractionalOrder) – Fractional order

• h (float) – Step size

• **kwargs – Additional parameters

Returns:

Tuple of (t_values, y_values)

Return type:

Tuple[ndarray, ndarray]

_euler_step(f, t_values, y_values, n, alpha, h)

Minimal fractional Euler update for 0<α≤1.

y_n = y_{n-1} + h^α / Γ(α+1) * f(t_{n-1}, y_{n-1})

Parameters:

• f (Callable)

• t_values (ndarray)

• y_values (ndarray)

• n (int)

• alpha (float | FractionalOrder)

• h (float)

Return type:

ndarray

_adams_bashforth_step(f, t_values, y_values, n, alpha, coeffs, h)

Compute one fractional Adams-Bashforth-like explicit step.

Uses available history with precomputed coefficients. For the first step, it falls back to Euler behavior.

Parameters:

• f (Callable)

• t_values (ndarray)

• y_values (ndarray)

• n (int)

• alpha (float | FractionalOrder)

• coeffs (ndarray)

• h (float)

Return type:

ndarray

_runge_kutta_step(f, t_values, y_values, n, alpha, h)

Compute one fractional RK2-style step.

Uses midpoint evaluation with fractional scaling factor.

Parameters:

• f (Callable)

• t_values (ndarray)

• y_values (ndarray)

• n (int)

• alpha (float | FractionalOrder)

• h (float)

Return type:

ndarray

_compute_fractional_coefficients(alpha, N)

Compute fractional derivative coefficients.

Parameters:

• alpha (float | FractionalOrder) – Fractional order

• N (int) – Number of time steps

Returns:

Array of coefficients

Return type:

ndarray

hpfracc.solvers.ode_solvers.solve_fractional_ode(f, t_span, y0, alpha, derivative_type='caputo', method='predictor_corrector', adaptive=False, h=None, **kwargs)

Solve fractional ODE.

Parameters:

• f (Callable)

• t_span (Tuple[float, float])

• y0 (float | ndarray)

• alpha (float | FractionalOrder)

• derivative_type (str)

• method (str)

• adaptive (bool)

• h (float | None)

Return type:

Tuple[ndarray, ndarray]

hpfracc.solvers.ode_solvers.solve_fractional_system(f, t_span, y0, alpha, derivative_type='caputo', method='predictor_corrector', **kwargs)

Solve system of fractional ODEs.

Parameters:

• f (Callable) – Right-hand side function f(t, y)

• t_span (Tuple[float, float]) – Time interval (t0, tf)

• y0 (ndarray) – Initial conditions

• alpha (float | ndarray) – Fractional order (scalar only). Per-component orders are not implemented.

• derivative_type (str) – Type of fractional derivative

• method (str) – Numerical method

• **kwargs – Additional solver parameters

Returns:

Tuple of (t_values, y_values)

Return type:

Tuple[ndarray, ndarray]

PDE Solvers

Fractional Partial Differential Equation Solvers

This module provides comprehensive solvers for fractional PDEs including finite difference methods, spectral methods, and adaptive mesh refinement.

Performance Note (v2.1.0): - Intelligent backend selection for sparse matrix operations - Optimal array operations based on problem size - Grünwald–Letnikov interior spatial operator is assembled directly in CSR (no dense fill)

hpfracc.solvers.pde_solvers._get_intelligent_selector()

Get intelligent backend selector instance for PDE solvers.

hpfracc.solvers.pde_solvers._select_array_backend(data_size, operation_type='element_wise')

Select optimal backend for array operations in PDE solving.

Parameters:

• data_size (int) – Number of elements to process

• operation_type (str) – Type of operation (e.g., “element_wise”, “matrix_multiply”)

Returns:

“numpy”, “numba”, or “jax”

Return type:

Backend name

hpfracc.solvers.pde_solvers._select_fft_backend(data_size)

Select optimal FFT backend.

Parameters:

data_size (int)

Return type:

str

hpfracc.solvers.pde_solvers._fft_convolution(coeffs, values, axis=0)

Fast convolution with backend selection.

Parameters:

• coeffs (ndarray)

• values (ndarray)

• axis (int)

Return type:

ndarray

hpfracc.solvers.pde_solvers._grunwald_letnikov_coeffs_pde(order, n_points)

Grünwald–Letnikov coefficients c_k for spatial discretization (recursive form).

Parameters:

• order (float)

• n_points (int)

Return type:

ndarray

hpfracc.solvers.pde_solvers._grunwald_letnikov_interior_sparse(nx, beta, dx)

Sparse interior Grünwald–Letnikov spatial operator (Toeplitz lower triangular).

Interior unknowns correspond to physical indices 1 .. nx-2; entry (i, j) uses offset k = i - j >= 0 with GL coefficient c_k, scaled by dx**(-beta).

Parameters:

• nx (int)

• beta (float | FractionalOrder)

• dx (float)

Return type:

csr_matrix

class hpfracc.solvers.pde_solvers.FractionalPDESolver(pde_type='diffusion', method='finite_difference', spatial_order=2, temporal_order=1, adaptive=False, *, fractional_order=None, boundary_conditions=None)

Bases: object

Base class for fractional PDE solvers.

Provides common functionality for solving fractional partial differential equations of various types.

Parameters:

• pde_type (str)

• method (str)

• spatial_order (int)

• temporal_order (int)

• adaptive (bool)

• fractional_order (float | FractionalOrder | None)

• boundary_conditions (str | None)

__init__(pde_type='diffusion', method='finite_difference', spatial_order=2, temporal_order=1, adaptive=False, *, fractional_order=None, boundary_conditions=None)

Initialize fractional PDE solver.

Parameters:

• pde_type (str) – Type of PDE (“diffusion”, “advection”, “reaction_diffusion”)

• method (str) – Numerical method (“finite_difference”, “spectral”, “finite_element”)

• spatial_order (int) – Order of spatial discretization

• temporal_order (int) – Order of temporal discretization

• adaptive (bool) – Use adaptive mesh refinement

• fractional_order (float | FractionalOrder | None)

• boundary_conditions (str | None)

_validate_orders(alpha, beta)

Validate fractional orders for PDE solver.

Parameters:

• alpha (float)

• beta (float)

_grunwald_letnikov_coeffs(order, n_points)

Compute Grünwald–Letnikov coefficients (shared by FD spatial operators).

Parameters:

• order (float)

• n_points (int)

Return type:

ndarray

_get_spatial_operator(nx, beta, dx)

Interior Grünwald–Letnikov operator as CSR (no dense fill).

Parameters:

• nx (int)

• beta (float | FractionalOrder)

• dx (float)

Return type:

csr_matrix

_build_spatial_matrix(nx, beta, diffusion_coeff, dx, dt, alpha, side='lhs')

Build spatial discretization matrix for interior unknowns (sparse assembly).

Combines identity and scaled Grünwald–Letnikov operator without forming a dense intermediate.

Parameters:

• nx (int)

• beta (float | FractionalOrder)

• diffusion_coeff (float)

• dx (float)

• dt (float)

• alpha (float | FractionalOrder)

• side (str)

Return type:

spmatrix

solve(t_span, x_span, initial_condition, boundary_conditions, alpha, beta, **kwargs)

Solve the fractional PDE.

Parameters:

• t_span (Tuple[float, float])

• x_span (Tuple[float, float])

• initial_condition (Callable)

• boundary_conditions (Dict)

• alpha (float)

• beta (float)

Return type:

Dict[str, ndarray]

class hpfracc.solvers.pde_solvers.FractionalDiffusionSolver(method='finite_difference', spatial_order=2, temporal_order=1, derivative_type='caputo')

Bases: FractionalPDESolver

Solver for fractional diffusion equations.

Solves equations of the form: ∂^α u/∂t^α = D ∂^β u/∂x^β + f(x, t, u)

where α and β are fractional orders.

Parameters:

• method (str)

• spatial_order (int)

• temporal_order (int)

• derivative_type (str)

__init__(method='finite_difference', spatial_order=2, temporal_order=1, derivative_type='caputo')

Initialize fractional diffusion solver.

Parameters:

• method (str) – Numerical method

• spatial_order (int) – Order of spatial discretization

• temporal_order (int) – Order of temporal discretization

• derivative_type (str) – Type of fractional derivative

solve(x_span, t_span, initial_condition, boundary_conditions, alpha, beta, diffusion_coeff=1.0, source_term=None, nx=100, nt=100, **kwargs)

Solve fractional diffusion equation.

Parameters:

• x_span (Tuple[float, float]) – Spatial interval (x0, xf)

• t_span (Tuple[float, float]) – Time interval (t0, tf)

• initial_condition (Callable) – Initial condition u(x, 0)

• boundary_conditions (Tuple[Callable, Callable]) – Boundary conditions (left_bc, right_bc)

• alpha (float | FractionalOrder) – Temporal fractional order

• beta (float | FractionalOrder) – Spatial fractional order

• diffusion_coeff (float) – Diffusion coefficient

• source_term (Callable | None) – Source term f(x, t, u)

• nx (int) – Number of spatial grid points

• nt (int) – Number of temporal grid points

• **kwargs – Additional solver parameters

Returns:

Tuple of (t_values, x_values, solution)

Return type:

Tuple[ndarray, ndarray, ndarray]

_solve_spatial_step(solution, n, x_values, t_values, alpha, beta, source_term, nx, nt, pde_params)

Solve spatial problem at current time step.

Parameters:

• solution (ndarray) – Solution array

• n (int) – Current time step

• x_values (ndarray) – Spatial grid

• t_values (ndarray) – Time grid

• alpha (float | FractionalOrder) – Temporal fractional order

• beta (float | FractionalOrder) – Spatial fractional order

• source_term (Callable | None) – Source term

• nx (int) – Number of spatial grid points

• nt (int) – Number of temporal grid points

• pde_params (dict) – Dictionary of parameters (diffusion_coeff, dx, dt)

Returns:

Solution at interior points

Return type:

ndarray

_finite_difference_step(solution, n, x_values, t_values, alpha, beta, source_term, nx, nt, pde_params)

Finite difference step.

Parameters:

• solution (ndarray) – Solution array

• n (int) – Current time step

• x_values (ndarray) – Spatial grid

• t_values (ndarray) – Time grid

• alpha (float | FractionalOrder) – Temporal fractional order

• beta (float | FractionalOrder) – Spatial fractional order

• source_term (Callable | None) – Source term

• nx (int) – Number of spatial grid points

• nt (int) – Number of temporal grid points

• pde_params (dict) – Dictionary of parameters (diffusion_coeff, dx, dt)

Returns:

Solution at interior points

Return type:

ndarray

_spectral_step(solution, n, x_values, t_values, alpha, beta, source_term, nx, nt, pde_params)

Correct, efficient, diagonal spectral solve for periodic BCs.

_compute_temporal_derivative(solution, n, alpha, dt, **kwargs)

Computes the RHS vector of the system, based on previous time steps.

Parameters:

• solution (ndarray)

• n (int)

• alpha (float | FractionalOrder)

• dt (float)

Return type:

ndarray

_compute_spatial_derivative(u, beta, dx)

Compute spatial fractional derivative.

Parameters:

• u (ndarray) – Solution at current time

• beta (float | FractionalOrder) – Spatial fractional order

• dx (float) – Spatial step size

Returns:

Spatial derivative

Return type:

ndarray

class hpfracc.solvers.pde_solvers.FractionalAdvectionSolver(method='finite_difference', spatial_order=2, temporal_order=1, derivative_type='caputo')

Bases: FractionalPDESolver

Solver for fractional advection-diffusion equations.

Solves equations of the form:

D_t^α u = v * D_x^β u

Currently, only integer orders (alpha=1, beta=1) are supported.

Parameters:

• method (str)

• spatial_order (int)

• temporal_order (int)

• derivative_type (str)

__init__(method='finite_difference', spatial_order=2, temporal_order=1, derivative_type='caputo')

Initialize fractional PDE solver.

Parameters:

• pde_type – Type of PDE (“diffusion”, “advection”, “reaction_diffusion”)

• method (str) – Numerical method (“finite_difference”, “spectral”, “finite_element”)

• spatial_order (int) – Order of spatial discretization

• temporal_order (int) – Order of temporal discretization

• adaptive – Use adaptive mesh refinement

• derivative_type (str)

solve(t_span, x_span, initial_condition, boundary_conditions, alpha, beta, **kwargs)

Solve the fractional advection equation.

Parameters:

• t_span (Tuple[float, float])

• x_span (Tuple[float, float])

• initial_condition (Callable)

• boundary_conditions (Dict)

• alpha (float)

• beta (float)

Return type:

Tuple[ndarray, ndarray, ndarray]

class hpfracc.solvers.pde_solvers.FractionalReactionDiffusionSolver(method='finite_difference', spatial_order=2, temporal_order=1, derivative_type='caputo')

Bases: FractionalPDESolver

Solver for fractional reaction-diffusion equations.

Solves equations of the form: ∂^α u/∂t^α = D ∂^β u/∂x^β + R(u) + f(x, t, u)

where α and β are fractional orders and R(u) is the reaction term.

Parameters:

• method (str)

• spatial_order (int)

• temporal_order (int)

• derivative_type (str)

__init__(method='finite_difference', spatial_order=2, temporal_order=1, derivative_type='caputo')

Initialize fractional reaction-diffusion solver.

Parameters:

• method (str) – Numerical method

• spatial_order (int) – Order of spatial discretization

• temporal_order (int) – Order of temporal discretization

• derivative_type (str) – Type of fractional derivative

solve(x_span, t_span, initial_condition, boundary_conditions, alpha, beta, diffusion_coeff=1.0, reaction_term=None, source_term=None, nx=100, nt=100, **kwargs)

Solve fractional reaction-diffusion equation.

Parameters:

• x_span (Tuple[float, float]) – Spatial interval (x0, xf)

• t_span (Tuple[float, float]) – Time interval (t0, tf)

• initial_condition (Callable) – Initial condition u(x, 0)

• boundary_conditions (Tuple[Callable, Callable]) – Boundary conditions (left_bc, right_bc)

• alpha (float | FractionalOrder) – Temporal fractional order

• beta (float | FractionalOrder) – Spatial fractional order

• diffusion_coeff (float) – Diffusion coefficient

• reaction_term (Callable | None) – Reaction term R(u)

• source_term (Callable | None) – Source term f(x, t, u)

• nx (int) – Number of spatial grid points

• nt (int) – Number of temporal grid points

• **kwargs – Additional solver parameters

Returns:

Tuple of (t_values, x_values, solution)

Return type:

Tuple[ndarray, ndarray, ndarray]

hpfracc.solvers.pde_solvers.solve_fractional_diffusion(x_span, t_span, initial_condition, boundary_conditions, alpha, beta, diffusion_coeff=1.0, source_term=None, method='finite_difference', nx=100, nt=100, **kwargs)

Solve fractional diffusion equation.

Parameters:

• x_span (Tuple[float, float]) – Spatial interval (x0, xf)

• t_span (Tuple[float, float]) – Time interval (t0, tf)

• initial_condition (Callable) – Initial condition u(x, 0)

• boundary_conditions (Tuple[Callable, Callable]) – Boundary conditions (left_bc, right_bc)

• alpha (float | FractionalOrder) – Temporal fractional order

• beta (float | FractionalOrder) – Spatial fractional order

• diffusion_coeff (float) – Diffusion coefficient

• source_term (Callable | None) – Source term f(x, t, u)

• method (str) – Numerical method

• nx (int) – Number of spatial grid points

• nt (int) – Number of temporal grid points

• **kwargs – Additional solver parameters

Returns:

Tuple of (t_values, x_values, solution)

Return type:

Tuple[ndarray, ndarray, ndarray]

hpfracc.solvers.pde_solvers.solve_fractional_advection(x_span, t_span, initial_condition, boundary_conditions, alpha, beta, velocity=1.0, source_term=None, method='finite_difference', nx=100, nt=100, **kwargs)

Solve fractional advection equation.

Parameters:

• x_span (Tuple[float, float]) – Spatial interval (x0, xf)

• t_span (Tuple[float, float]) – Time interval (t0, tf)

• initial_condition (Callable) – Initial condition u(x, 0)

• boundary_conditions (Tuple[Callable, Callable]) – Boundary conditions (left_bc, right_bc)

• alpha (float | FractionalOrder) – Temporal fractional order

• beta (float | FractionalOrder) – Spatial fractional order

• velocity (float) – Advection velocity

• source_term (Callable | None) – Source term f(x, t, u)

• method (str) – Numerical method

• nx (int) – Number of spatial grid points

• nt (int) – Number of temporal grid points

• **kwargs – Additional solver parameters

Returns:

Tuple of (t_values, x_values, solution)

Return type:

Tuple[ndarray, ndarray, ndarray]

hpfracc.solvers.pde_solvers.solve_fractional_reaction_diffusion(x_span, t_span, initial_condition, boundary_conditions, alpha, beta, diffusion_coeff=1.0, reaction_term=None, source_term=None, method='finite_difference', nx=100, nt=100, **kwargs)

Solve fractional reaction-diffusion equation.

Parameters:

• x_span (Tuple[float, float]) – Spatial interval (x0, xf)

• t_span (Tuple[float, float]) – Time interval (t0, tf)

• initial_condition (Callable) – Initial condition u(x, 0)

• boundary_conditions (Tuple[Callable, Callable]) – Boundary conditions (left_bc, right_bc)

• alpha (float | FractionalOrder) – Temporal fractional order

• beta (float | FractionalOrder) – Spatial fractional order

• diffusion_coeff (float) – Diffusion coefficient

• reaction_term (Callable | None) – Reaction term R(u)

• source_term (Callable | None) – Source term f(x, t, u)

• method (str) – Numerical method

• nx (int) – Number of spatial grid points

• nt (int) – Number of temporal grid points

• **kwargs – Additional solver parameters

Returns:

Tuple of (t_values, x_values, solution)

Return type:

Tuple[ndarray, ndarray, ndarray]

hpfracc.solvers.pde_solvers.solve_fractional_pde(x_span, t_span, initial_condition, boundary_conditions, alpha, beta, equation_type='diffusion', **kwargs)

Generic solver for fractional PDEs.

Parameters:

• x_span (Tuple[float, float]) – Spatial interval (x0, xf)

• t_span (Tuple[float, float]) – Time interval (t0, tf)

• initial_condition (Callable) – Initial condition u(x, 0)

• boundary_conditions (Tuple[Callable, Callable]) – Boundary conditions (left_bc, right_bc)

• alpha (float | FractionalOrder) – Temporal fractional order

• beta (float | FractionalOrder) – Spatial fractional order

• equation_type (str) – Type of PDE (“diffusion”, “advection”, “reaction_diffusion”)

• **kwargs – Additional solver parameters

Returns:

Tuple of (t_values, x_values, solution)

Return type:

Tuple[ndarray, ndarray, ndarray]

SDE Solvers

Fractional Stochastic Differential Equation Solvers

This module provides comprehensive solvers for fractional SDEs including various numerical methods, adaptive step size control, and error estimation.

Performance Note: - Uses FFT-based convolution for O(N log N) history summation instead of O(N²) - Intelligent backend selection for optimal performance - Multi-backend support (PyTorch, JAX, NumPy/Numba)

hpfracc.solvers.sde_solvers._get_gamma_function()

Get gamma function through adapter system.

hpfracc.solvers.sde_solvers._get_intelligent_selector()

Get intelligent backend selector instance.

hpfracc.solvers.sde_solvers._fft_convolution(coeffs, values, axis=0)

Fast convolution using FFT for O(N log N) performance.

Parameters:

• coeffs (ndarray) – Coefficient array (1D)

• values (ndarray) – Value array (can be 1D or 2D)

• axis (int) – Axis along which to perform convolution

Returns:

Convolution result (same shape as values)

Return type:

ndarray

class hpfracc.solvers.sde_solvers.SDESolution(t, y, fractional_order, method, drift_func, diffusion_func, metadata=None)

Bases: object

Solution object for SDE solvers.

Parameters:

• t (ndarray)

• y (ndarray)

• fractional_order (float | FractionalOrder)

• method (str)

• drift_func (Callable)

• diffusion_func (Callable)

• metadata (Dict[str, Any])

t: ndarray

y: ndarray

fractional_order: float | FractionalOrder

method: str

drift_func: Callable

diffusion_func: Callable

metadata: Dict[str, Any] = None

__init__(t, y, fractional_order, method, drift_func, diffusion_func, metadata=None)

Parameters:

• t (ndarray)

• y (ndarray)

• fractional_order (float | FractionalOrder)

• method (str)

• drift_func (Callable)

• diffusion_func (Callable)

• metadata (Dict[str, Any])

Return type:

None

class hpfracc.solvers.sde_solvers.FractionalSDESolver(fractional_order, definition='caputo', adaptive=False, rtol=1e-05, atol=1e-08)

Bases: ABC

Base class for fractional SDE solvers.

A fractional SDE takes the form:

D^α X(t) = f(t, X(t)) dt + g(t, X(t)) dW(t)

where:

• α is the fractional order

• f is the drift function

• g is the diffusion function

• W(t) is a Wiener process

Parameters:

• fractional_order (float | FractionalOrder)

• definition (str)

• adaptive (bool)

• rtol (float)

• atol (float)

__init__(fractional_order, definition='caputo', adaptive=False, rtol=1e-05, atol=1e-08)

Initialize fractional SDE solver.

Parameters:

• fractional_order (float | FractionalOrder) – Fractional order (0 < α < 2)

• definition (str) – Type of fractional derivative (“caputo” or “riemann_liouville”)

• adaptive (bool) – Use adaptive step size

• rtol (float) – Relative tolerance for adaptive stepping

• atol (float) – Absolute tolerance for adaptive stepping

abstractmethod solve(drift, diffusion, x0, t_span, noise_model=None, **kwargs)

Solve fractional SDE.

Parameters:

• drift (Callable[[float, ndarray], ndarray]) – Drift function f(t, x) -> R^d

• diffusion (Callable[[float, ndarray], ndarray]) – Diffusion function g(t, x) -> R^(d x m) where m is dimension of noise

• x0 (ndarray) – Initial condition

• t_span (Tuple[float, float]) – Time interval (t0, tf)

• **kwargs – Additional solver parameters

• noise_model (NoiseModel | None)

Returns:

SDESolution object containing trajectory

Return type:

SDESolution

class hpfracc.solvers.sde_solvers.FastHistoryConvolution(alpha, num_steps, dim)

Bases: object

Helper class for efficient history convolution in fractional SDE solvers.

Features: - Pre-allocated arrays to avoid list overhead (O(N) -> O(1) allocation) - Dynamic switching between direct dot product (small N) and FFT (large N) - JAX/SciPy/NumPy backend support via _fft_convolution

Parameters:

• alpha (float)

• num_steps (int)

• dim (int)

__init__(alpha, num_steps, dim)

Parameters:

• alpha (float)

• num_steps (int)

• dim (int)

update(value)

Add new value to history.

Parameters:

value (ndarray)

convolve()

Compute convolution of weights with history.

Return type:

ndarray

hpfracc.solvers.sde_solvers._milstein_fd_eps(x, base=1e-06)

Scale-aware finite-difference step for ∂g/∂x.

Parameters:

• x (ndarray)

• base (float)

Return type:

float

hpfracc.solvers.sde_solvers._milstein_correction_diagonal(diffusion, t, x, g_vec, dW, dt)

Diagonal / independent-noise Milstein correction:

0.5 * sum_i g_i(t,x) * (∂g_i/∂x_i) * (dW_i^2 - dt)

with ∂g_i/∂x_i approximated by central differences on x_i.

Parameters:

• diffusion (Callable[[float, ndarray], ndarray])

• t (float)

• x (ndarray)

• g_vec (ndarray)

• dW (ndarray)

• dt (float)

Return type:

ndarray

hpfracc.solvers.sde_solvers._milstein_correction_single_wiener(diffusion, t, x, g_vec, dW, dt)

Milstein correction for a single driving Wiener process (m=1, dW scalar):

0.5 * (dW^2 - dt) * sum_j g_j * (∂g_k/∂x_j) for each component k.

Uses central differences in the state directions x_j.

Parameters:

• diffusion (Callable[[float, ndarray], ndarray])

• t (float)

• x (ndarray)

• g_vec (ndarray)

• dW (ndarray)

• dt (float)

Return type:

ndarray

class hpfracc.solvers.sde_solvers.FractionalEulerMaruyama(*args, **kwargs)

Bases: FractionalSDESolver

Fractional Euler-Maruyama method for solving fractional SDEs.

This is a first-order strong convergence method.

__init__(*args, **kwargs)

Initialize fractional SDE solver.

Parameters:

• fractional_order – Fractional order (0 < α < 2)

• definition – Type of fractional derivative (“caputo” or “riemann_liouville”)

• adaptive – Use adaptive step size

• rtol – Relative tolerance for adaptive stepping

• atol – Absolute tolerance for adaptive stepping

solve(drift, diffusion, x0, t_span, num_steps=100, seed=None, noise_model=None, **kwargs)

Solve fractional SDE using Euler-Maruyama method.

Parameters:

• drift (Callable) – Drift function f(t, x)

• diffusion (Callable) – Diffusion function g(t, x)

• x0 (ndarray) – Initial condition

• t_span (Tuple[float, float]) – Time interval (t0, tf)

• num_steps (int) – Number of time steps

• seed (int | None) – Random seed for Wiener process

• noise_model (NoiseModel | None) – Optional noise model for non-Brownian noise

Returns:

SDESolution object

Return type:

SDESolution

class hpfracc.solvers.sde_solvers.FractionalMilstein(*args, **kwargs)

Bases: FractionalSDESolver

Fractional SDE integrator with the same history-convolution path as Euler–Maruyama, plus a Milstein correction where supported.

The correction uses the standard form 0.5 * (dW^2 - dt) multiplied by terms involving spatial derivatives of the diffusion, approximated by central finite differences.

Supported

• Diagonal noise: g(t,x) returns a vector of length d; one Brownian increment per component (same layout as Euler–Maruyama).

• Single Wiener process: g(t,x) returns shape (d, 1); scalar dW.

Not implemented

• General matrix diffusion with m > 1 independent Wiener processes (Levy area / non-commutative terms not included). Use diagonal or scalar noise, or extend the solver.

__init__(*args, **kwargs)

Initialize fractional SDE solver.

Parameters:

• fractional_order – Fractional order (0 < α < 2)

• definition – Type of fractional derivative (“caputo” or “riemann_liouville”)

• adaptive – Use adaptive step size

• rtol – Relative tolerance for adaptive stepping

• atol – Absolute tolerance for adaptive stepping

solve(drift, diffusion, x0, t_span, num_steps=100, seed=None, noise_model=None, **kwargs)

Integrate with fractional history terms and Milstein correction (supported cases).

Parameters:

• drift (Callable)

• diffusion (Callable)

• x0 (ndarray)

• t_span (Tuple[float, float])

• num_steps (int)

• seed (int | None)

• noise_model (NoiseModel | None)

Return type:

SDESolution

hpfracc.solvers.sde_solvers.solve_fractional_sde(drift, diffusion, x0, t_span, fractional_order=0.5, method='euler_maruyama', num_steps=100, seed=None, noise_model=None, **kwargs)

Solve a fractional SDE.

Parameters:

• drift (Callable) – Drift function f(t, x) -> R^d

• diffusion (Callable) – Diffusion function g(t, x) -> R^(d x m)

• x0 (ndarray) – Initial condition

• t_span (Tuple[float, float]) – Time interval (t0, tf)

• fractional_order (float | FractionalOrder) – Fractional order (0 < α < 2)

• method (str) – Solver method (euler_maruyama or milstein). milstein adds a finite-difference Milstein correction for diagonal noise or a single Wiener driver; see FractionalMilstein.

• num_steps (int) – Number of time steps

• seed (int | None) – Random seed for Wiener process

• **kwargs – Additional solver parameters

• noise_model (NoiseModel | None)

Returns:

SDESolution object

Return type:

SDESolution

Example

>>> def drift(t, x):
...     return -0.5 * x
>>> def diffusion(t, x):
...     return 0.2 * np.eye(1)
>>> x0 = np.array([1.0])
>>> sol = solve_fractional_sde(drift, diffusion, x0, (0, 1), 0.5, num_steps=100)
>>> print(sol.y[-1])

hpfracc.solvers.sde_solvers.solve_fractional_sde_system(drift, diffusion, x0, t_span, fractional_order, method='euler_maruyama', noise_type='additive', num_steps=100, seed=None, noise_model=None, **kwargs)

Solve a system of coupled fractional SDEs.

Parameters:

• drift (Callable) – Drift function f(t, x) -> R^d

• diffusion (Callable) – Diffusion function g(t, x) -> R^(d x m)

• x0 (ndarray) – Initial condition

• t_span (Tuple[float, float]) – Time interval (t0, tf)

• fractional_order (float | FractionalOrder | List[float]) – Fractional order(s) for system

• method (str) – Solver method

• noise_type (str) – Type of noise (“additive” or “multiplicative”)

• num_steps (int) – Number of time steps

• seed (int | None) – Random seed

• noise_model (NoiseModel | None) – Optional noise model

• **kwargs – Additional parameters

Returns:

SDESolution object

Return type:

SDESolution

Noise Models

Stochastic Noise Models for Fractional SDEs

This module provides various noise models including Brownian motion, fractional Brownian motion, Lévy processes, and coloured noise.

class hpfracc.solvers.noise_models.NoiseModel

Bases: ABC

Base class for stochastic noise models.

abstractmethod generate_increment(t, dt, size=(), seed=None)

Generate noise increment.

Parameters:

• t (float) – Current time

• dt (float) – Time step

• size (Tuple[int, ...]) – Output shape

• seed (int | None) – Random seed

Returns:

Noise increment array

Return type:

ndarray

prepare(num_steps, dt, size=())

Optional preparation step for pre-computing noise paths. Useful for non-Markovian processes like fBm.

Parameters:

• num_steps (int)

• dt (float)

• size (Tuple[int, ...])

class hpfracc.solvers.noise_models.BrownianMotion(scale=1.0)

Bases: NoiseModel

Standard Brownian motion (Wiener process).

Generates independent Gaussian increments with variance dt.

Parameters:

scale (float)

__init__(scale=1.0)

Initialize Brownian motion.

Parameters:

scale (float) – Scaling factor for noise

generate_increment(t, dt, size=(), seed=None)

Generate Wiener increment.

Parameters:

• t (float)

• dt (float)

• size (Tuple[int, ...])

• seed (int | None)

Return type:

ndarray

variance(dt)

Get variance of increment.

Parameters:

dt (float)

Return type:

float

class hpfracc.solvers.noise_models.FractionalBrownianMotion(hurst=0.5, scale=1.0)

Bases: NoiseModel

Fractional Brownian motion (fBm).

A Gaussian process with long-range dependence characterized by the Hurst exponent H (0 < H < 1).

Parameters:

• hurst (float)

• scale (float)

__init__(hurst=0.5, scale=1.0)

Initialize fractional Brownian motion.

Parameters:

• hurst (float) – Hurst exponent (H=0.5 gives standard Brownian motion)

• scale (float) – Scaling factor for noise

prepare(num_steps, dt, size=())

Pre-compute fBm increments using Davies-Harte method (FFT).

Parameters:

• num_steps (int)

• dt (float)

• size (Tuple[int, ...])

generate_increment(t, dt, size=(), seed=None)

Generate fractional Brownian motion increment.

Uses precomputed Davies-Harte increments if available, otherwise falls back to simplified independent increments (with warning).

Parameters:

• t (float)

• dt (float)

• size (Tuple[int, ...])

• seed (int | None)

Return type:

ndarray

property is_standard_bm: bool

Check if this is standard Brownian motion.

class hpfracc.solvers.noise_models.LevyNoise(alpha=1.5, beta=0.0, scale=1.0, location=0.0)

Bases: NoiseModel

Lévy noise for jump diffusions.

Uses stable distributions to model heavy-tailed noise.

Parameters:

• alpha (float)

• beta (float)

• scale (float)

• location (float)

__init__(alpha=1.5, beta=0.0, scale=1.0, location=0.0)

Initialize Lévy noise.

Parameters:

• alpha (float) – Stability parameter (0 < α ≤ 2)

• beta (float) – Skewness parameter (-1 ≤ β ≤ 1)

• scale (float) – Scale parameter

• location (float) – Location parameter

generate_increment(t, dt, size=(), seed=None)

Generate Lévy noise increment.

Uses stable distribution sampling (simplified implementation).

Parameters:

• t (float)

• dt (float)

• size (Tuple[int, ...])

• seed (int | None)

Return type:

ndarray

class hpfracc.solvers.noise_models.ColouredNoise(correlation_time=1.0, amplitude=1.0, seed=None)

Bases: NoiseModel

Coloured noise (Ornstein-Uhlenbeck process).

Gaussian noise with exponential autocorrelation.

Parameters:

• correlation_time (float)

• amplitude (float)

• seed (int | None)

__init__(correlation_time=1.0, amplitude=1.0, seed=None)

Initialize coloured noise.

Parameters:

• correlation_time (float) – Correlation time constant

• amplitude (float) – Noise amplitude

• seed (int | None) – Random seed for state initialization

generate_increment(t, dt, size=(), seed=None)

Generate coloured noise increment.

Parameters:

• t (float)

• dt (float)

• size (Tuple[int, ...])

• seed (int | None)

Return type:

ndarray

reset()

Reset the noise process state.

hpfracc.solvers.noise_models.numpyro_brownian_noise(t, dt, scale=1.0)

NumPyro model for Brownian noise.

Parameters:

• t (float) – Current time

• dt (float) – Time step

• scale (float) – Noise scale

Returns:

NumPyro sample

hpfracc.solvers.noise_models.numpyro_fractional_brownian_noise(t, dt, hurst=0.5, scale=1.0)

NumPyro model for fractional Brownian noise.

Parameters:

• t (float) – Current time

• dt (float) – Time step

• hurst (float) – Hurst exponent

• scale (float) – Noise scale

Returns:

NumPyro sample

hpfracc.solvers.noise_models.numpyro_levy_noise(t, dt, alpha=1.5, scale=1.0)

NumPyro model for Lévy noise.

Note: NumPyro doesn’t have stable distribution, so uses approximation.

Parameters:

• t (float) – Current time

• dt (float) – Time step

• alpha (float) – Stability parameter

• scale (float) – Noise scale

Returns:

NumPyro sample

class hpfracc.solvers.noise_models.NoiseConfig(noise_type='brownian', hurst=0.5, scale=1.0, alpha=1.5, beta=0.0, correlation_time=1.0, amplitude=1.0)

Bases: object

Configuration for noise models.

Parameters:

• noise_type (str)

• hurst (float)

• scale (float)

• alpha (float)

• beta (float)

• correlation_time (float)

• amplitude (float)

noise_type: str = 'brownian'

hurst: float = 0.5

scale: float = 1.0

alpha: float = 1.5

beta: float = 0.0

correlation_time: float = 1.0

amplitude: float = 1.0

__init__(noise_type='brownian', hurst=0.5, scale=1.0, alpha=1.5, beta=0.0, correlation_time=1.0, amplitude=1.0)

Parameters:

• noise_type (str)

• hurst (float)

• scale (float)

• alpha (float)

• beta (float)

• correlation_time (float)

• amplitude (float)

Return type:

None

hpfracc.solvers.noise_models.create_noise_model(config)

Create a noise model from configuration.

Parameters:

config (NoiseConfig) – Noise configuration

Returns:

NoiseModel instance

Return type:

NoiseModel

hpfracc.solvers.noise_models.generate_noise_trajectory(noise_model, t_span, num_steps, size=(), seed=None)

Generate a complete noise trajectory.

Parameters:

• noise_model (NoiseModel) – Noise model to use

• t_span (Tuple[float, float]) – Time interval (t0, tf)

• num_steps (int) – Number of steps

• size (Tuple[int, ...]) – Shape of noise increments

• seed (int | None) – Random seed

Returns:

Tuple of (time array, noise increments array)

Return type:

Tuple[ndarray, ndarray]

Coupled graph–SDE solvers

Coupled System Solvers for Graph-SDE Dynamics

This module provides numerical solvers for systems of coupled spatial-temporal dynamics, integrating graph-based spatial evolution with fractional SDE temporal evolution.

Modeling scope (read before interpreting results)

• Spatial track: explicit Euler on graph_dynamics(state, adjacency) (one step per substep). No fractional derivative in space unless graph_dynamics implements it.

• Temporal track: history convolutions via FastHistoryConvolution with gamma(alpha+1) scaling and a single dt**alpha factor per step—an engineering approximation to fractional SDE memory, not a full rigorous discretization proof.

• Coupling: both solvers add coupling_strength * (spatial - temporal) into the temporal drift (operator splitting: after the spatial half-step; monolithic: inside the combined drift). The spatial ODE does not receive an symmetric pull from temporal in OperatorSplittingSolver (only MonolithicSolver adds k*(V-U) to dspatial).

• API: coupling_type in solve_coupled_graph_sde() is reserved and not used; only coupling_strength and the solver choice affect dynamics.

• ``multiscale``: not implemented; raises ValueError.

For production experiments, validate against a reference integrator or reduce step size when coupling is strong; splitting error can dominate.

class hpfracc.solvers.coupled_solvers.CoupledSolution(t, spatial, temporal, coupling, metadata=None)

Bases: object

Solution object for coupled graph-SDE systems.

Parameters:

• t (ndarray)

• spatial (ndarray)

• temporal (ndarray)

• coupling (ndarray)

• metadata (Dict[str, Any])

t: ndarray

spatial: ndarray

temporal: ndarray

coupling: ndarray

metadata: Dict[str, Any] = None

__init__(t, spatial, temporal, coupling, metadata=None)

Parameters:

• t (ndarray)

• spatial (ndarray)

• temporal (ndarray)

• coupling (ndarray)

• metadata (Dict[str, Any])

Return type:

None

class hpfracc.solvers.coupled_solvers.CoupledSystemSolver(fractional_orders, coupling_strength=1.0)

Bases: ABC

Base class for coupled system solvers.

Parameters:

• fractional_orders (float | FractionalOrder | Dict[str, float])

• coupling_strength (float)

__init__(fractional_orders, coupling_strength=1.0)

Initialize coupled system solver.

Parameters:

• fractional_orders (float | FractionalOrder | Dict[str, float]) – Fractional order(s) for system

• coupling_strength (float) – Strength of spatial-temporal coupling

abstractmethod solve(graph_dynamics, sde_drift, sde_diffusion, adjacency, node_features, t_span, **kwargs)

Solve coupled system.

Parameters:

• graph_dynamics (Callable)

• sde_drift (Callable)

• sde_diffusion (Callable)

• adjacency (ndarray)

• node_features (ndarray)

• t_span (Tuple[float, float])

Return type:

CoupledSolution

class hpfracc.solvers.coupled_solvers.OperatorSplittingSolver(fractional_orders, coupling_strength=1.0, split_order=2)

Bases: CoupledSystemSolver

Operator splitting solver for graph-SDE dynamics.

Spatial evolution uses graph_dynamics on the spatial track. The temporal (fractional SDE) track adds an explicit coupling term to the drift:

drift_total = sde_drift(t, temporal) + coupling_strength * (spatial - temporal)

using the current spatial state after the preceding spatial substeps in the split. This matches the spirit of MonolithicSolver (coupling in the temporal drift).

Parameters:

• fractional_orders (float | FractionalOrder | Dict[str, float])

• coupling_strength (float)

• split_order (int)

__init__(fractional_orders, coupling_strength=1.0, split_order=2)

Initialize operator splitting solver.

Parameters:

• fractional_orders (float | FractionalOrder | Dict[str, float]) – Fractional order(s)

• coupling_strength (float) – Coupling strength

• split_order (int) – Splitting order (1=Lie-Trotter, 2=Strang)

solve(graph_dynamics, sde_drift, sde_diffusion, adjacency, node_features, t_span, num_steps=100, seed=None, **kwargs)

Solve using operator splitting.

For Strang splitting (order 2): - Half step of spatial dynamics - Full step of temporal dynamics - Half step of spatial dynamics

Parameters:

• graph_dynamics (Callable)

• sde_drift (Callable)

• sde_diffusion (Callable)

• adjacency (ndarray)

• node_features (ndarray)

• t_span (Tuple[float, float])

• num_steps (int)

• seed (int | None)

Return type:

CoupledSolution

_spatial_step(graph_dynamics, adjacency, state, dt)

Single spatial (graph) evolution step.

Parameters:

• graph_dynamics (Callable)

• adjacency (ndarray)

• state (ndarray)

• dt (float)

Return type:

ndarray

_temporal_step(drift, diffusion, state, t, dt, drift_conv, diffusion_conv, gamma_factor, alpha, initial_state, spatial_for_coupling, rng)

Temporal fractional SDE substep with explicit drift coupling to the spatial track.

Parameters:

• drift (Callable)

• diffusion (Callable)

• state (ndarray)

• t (float)

• dt (float)

• drift_conv (FastHistoryConvolution)

• diffusion_conv (FastHistoryConvolution)

• gamma_factor (float)

• alpha (float)

• initial_state (ndarray)

• spatial_for_coupling (ndarray)

• rng (Generator)

Return type:

ndarray

class hpfracc.solvers.coupled_solvers.MonolithicSolver(fractional_orders, coupling_strength=1.0)

Bases: CoupledSystemSolver

Monolithic solver for strongly coupled graph-SDE systems.

Solves the full coupled system simultaneously for better accuracy in strongly coupled regimes, at the cost of higher memory usage.

Parameters:

• fractional_orders (float | FractionalOrder | Dict[str, float])

• coupling_strength (float)

solve(graph_dynamics, sde_drift, sde_diffusion, adjacency, node_features, t_span, num_steps=100, seed=None, **kwargs)

Solve monolithic coupled system.

Parameters:

• graph_dynamics (Callable)

• sde_drift (Callable)

• sde_diffusion (Callable)

• adjacency (ndarray)

• node_features (ndarray)

• t_span (Tuple[float, float])

• num_steps (int)

• seed (int | None)

Return type:

CoupledSolution

hpfracc.solvers.coupled_solvers.solve_coupled_graph_sde(graph_dynamics, sde_drift, sde_diffusion, adjacency, node_features, t_span, fractional_orders=0.5, coupling_type='bidirectional', coupling_strength=1.0, solver='operator_splitting', **kwargs)

Solve coupled graph-SDE system.

See module docstring for limitations (unused coupling_type, approximate temporal fractional path, asymmetric coupling in operator splitting).

Parameters:

• graph_dynamics (Callable) – Spatial dynamics f(spatial, adjacency).

• sde_drift (Callable) – Temporal drift f(t, temporal).

• sde_diffusion (Callable) – Temporal diffusion g(t, temporal).

• adjacency (ndarray) – Graph adjacency matrix.

• node_features (ndarray) – Initial node features (shared as spatial/temporal initial shape).

• t_span (Tuple[float, float]) – Time interval.

• fractional_orders (float | FractionalOrder | Dict[str, float]) – Scalar or dict with 'spatial' / 'temporal' keys; temporal order drives FastHistoryConvolution.

• coupling_type (str) – Reserved for future use; currently ignored.

• coupling_strength (float) – Prefactor on (spatial - temporal) in the temporal drift.

• solver (str) – "operator_splitting" or "monolithic".

• **kwargs – Passed to the solver’s solve (e.g. num_steps, seed).

Returns:

CoupledSolution

Return type:

CoupledSolution

JAX / Diffrax utilities (optional)

Imported explicitly as hpfracc.solvers.solvers_jax_utils (not re-exported from hpfracc.solvers). Autodoc of this submodule is skipped here because optional JAX / Diffrax stacks and version-specific internals often fail on documentation builders; see docs/development/SOLVERS_ARCHITECTURE.md in the repository for exports and extras.