Mittag-Leffler Euler Integrator for Fractional SPDEs

• The Mittag-Leffler Euler Integrator is a numerical method for fractional evolution equations, using the Mittag-Leffler function to capture memory effects and anomalous diffusion.
• It generalizes the exponential Euler scheme to fractional-order settings, employing a two-parameter Mittag-Leffler function for improved convergence in stochastic Volterra equations.
• The method achieves nearly first-order strong convergence and supports rigorous error and large-deviation analyses, making it effective for complex fractional diffusion problems.

The Mittag-Leffler Euler integrator (MLEI) is a time-stepping numerical method tailored for evolution equations with fractional time dynamics, particularly those involving stochasticity and spatial operators with nonlocal memory kernels. MLEI generalizes the exponential Euler scheme to fractional-order settings, exploiting the smoothing and decay properties of the Mittag-Leffler function to achieve enhanced convergence rates in the discretization of semilinear stochastic Volterra equations and related space-time fractional diffusion problems (Kovács et al., 2018, Dai et al., 2022).

1. Mathematical Formulation and Definition

The MLEI is principally designed for SPDEs of the form

$$u(t) + \int_{0}^{t} b(t-s)\,Au(s)\,ds = F(u(t))\,dt + dW(t),$$

where $u(t)$ is the unknown process in a Hilbert space, $A$ is a (typically unbounded) linear operator, $b(t) = t^{\alpha-1}/\Gamma(\alpha)$ with $0<\alpha<1$ defines the memory kernel, $F$ is a nonlinear mapping, and $W(t)$ is a cylindrical Wiener process. When expressed with a Riemann--Liouville or Caputo fractional integral, the model becomes

$u(t) + J_0^\alpha(Au)(t) = F(u(t))\,dt + dW(t).$

For space-time fractional diffusion equations (with Caputo derivative and potentially nontrivial spatial operator powers), the archetype is

$$\partial_t^\alpha X(t) + A^\beta X(t) = F(X(t)) + {}_t^\gamma\dot{W}(t),\quad X(0)=X_0,$$

with $A=-\Delta$, $\beta\in(0,1]$, and the noise term given as a Riemann--Liouville fractional integral of Wiener process increments (Dai et al., 2022).

Discrete MLEI Scheme

Let $0=t_0<t_1<...<t_M=T$ be a uniform temporal grid with $\Delta t = T/M$. The MLEI algorithm computes approximations $U^N_m$ (after spectral Galerkin projection onto $N$ modes) as

$$U_{m+1}^N = S(\Delta t)\, U_m^N + \int_{t_m}^{t_{m+1}} S(t_{m+1}-s) P_N F(U_m^N) ds + \int_{t_m}^{t_{m+1}} S(t_{m+1}-s) P_N dW(s),$$

with the resolvent family $S(t) = E_{\alpha+1}(-A t^{\alpha+1})$ defined in terms of the Mittag-Leffler function (see Section 2). For implementation, the deterministic integral can be compressed via the two-parameter Mittag-Leffler function: $$U_{m+1}^N = E_{\alpha+1}(-A \Delta t^{\alpha+1}) U_m^N + \Delta t\, E_{\alpha+1,\alpha+2}(-A \Delta t^{\alpha+1}) P_N F(U_m^N) + \int_{t_m}^{t_{m+1}} E_{\alpha+1}(-A (t_{m+1}-s)^{\alpha+1}) P_N dW(s).$$ For time-fractional equations not involving explicit memory (e.g., space-time fractional diffusion), the mild MLEI variant reads

$$Y_0 = X_0,\qquad Y_m = {}_{1-\alpha}(t_m) Y_0 + \sum_{j=0}^{m-1} \int_{t_j}^{t_{j+1}} {}_0(t_m-s) F(Y_j)\,ds + \Lambda(t_m),$$

where the operators ${}_{\eta}(t)$ act via spectral expansions involving the two-parameter Mittag-Leffler function, and $\Lambda(t_m)$ is a (precomputable) Gaussian convolution (Dai et al., 2022).

2. Mittag-Leffler Functions and Fractional Resolvent Families

The Mittag-Leffler function $E_{a,b}(z)$ is defined by

$$E_{a,b}(z) = \sum_{n=0}^\infty \frac{z^n}{\Gamma(an+b)},$$

where $a>0$, $b>0$. This function generalizes the exponential ($a=1$), naturally encoding anomalous diffusion and fractional-order relaxation processes.

The fractional resolvent family associated with the operator $A$ is

$$S(t) = \sum_{k=1}^\infty E_{\alpha+1}(-\lambda_k t^{\alpha+1}) (\cdot,e_k) e_k,$$

with $(e_k)$ an orthonormal basis of eigenvectors and $\lambda_k$ the corresponding eigenvalues.

Essential smoothing estimates for the analysis are, for $\rho = \alpha + 1 \in (1,2)$, and $s\in[0,1/\rho]$,

$$\|A^s S(t)\| \le C t^{-s\rho},\quad \|A^s \dot{S}(t)\| \le C t^{-s\rho-1},\quad \|A^{-s} \dot{S}(t)\| \le C t^{s\rho-1}.$$

These are vital in establishing temporal and spatial error bounds in the strong convergence analysis (Kovács et al., 2018).

3. Error Analysis and Super-Convergence

The MLEI achieves a strong convergence order that is (almost) twice that of backward Euler-based convolution quadrature (BE-CQ) schemes for stochastic fractional problems. Theoretical results assert that, under additive noise of regularity $\beta\in(0,1/\rho]$,

$$\sup_{0\le m\le M} \|u(t_m)-U_m^N\|_{L^2(\Omega;H)} \le C\left(\lambda_{N+1}^{-\gamma/2} + \Delta t^{\gamma \rho}\right),$$

for initial data $u_0 \in L^4(\Omega; \dot{H}^\gamma)$, $\gamma<1/\rho$, and suitable conditions on $F$ (Kovács et al., 2018).

In the trace-class noise scenario ($\beta=1/\rho$), the time-convergence rate is $O(\Delta t^{1-\varepsilon})$, nearly first order. The BE-CQ method only provides $O(\Delta t^{1/2})$ for comparable noise regularity.

For space-time fractional diffusion equations, a refined decomposition of the singular time-fractional integral yields "super-convergence"—that is, the strong error order strictly exceeds the temporal Hölder regularity of the solution. The key is splitting the singular convolution into regular and singular parts, carefully exploiting the cancellation and the continuity properties of the nonlinearity. Analysis via telescopic splitting of local errors and singular Grönwall-type inequalities establishes enhanced order for both linear and nonlinear equations (Dai et al., 2022).

4. Comparative Assessment with Other Schemes

When compared with classical BE-CQ time discretizations, the MLEI leverages the explicit structure of fractional evolution via the resolvent family. The table below summarizes key contrasts:

Method	Temporal Rate (Trace-Class Noise)	Kernel Handling	Nonlinearity Freezing
MLEI	$O(\Delta t^{1-\varepsilon})$	Exact	Once per step
BE-CQ	$O(\Delta t^{1/2})$	Convolution	Implicit via linearization

MLEI thereby attains nearly first-order strong convergence in time, enabled by the full smoothing of the Mittag-Leffler kernel, but introduces computational requirements such as the evaluation of matrix Mittag-Leffler functions and efficient sampling of the resulting correlated stochastic increments (Kovács et al., 2018).

5. Numerical Implementation and Experimental Findings

Numerical validation involves single-mode and multi-dimensional tests for fractional orders $\rho = 1.2, 1.5, 1.75$, nonlinearities (e.g., $f(u) = \sin u$), and additive noise. Results consistently show, on log–log plots, a slope close to 1 for MLEI as $\Delta t \to 0$, while BE-CQ slopes approximate 0.5. Simulations use up to 100 Monte Carlo paths and high-resolution temporal reference solutions, confirming the theoretically predicted rates (Kovács et al., 2018).

6. Large Deviations and $\Gamma$-Convergence of Rate Functions

Extending beyond strong error analysis, the large-deviation behavior of MLEI for discretized small-noise SPDEs is established using the weak convergence approach of Budhiraja–Dupuis. Both continuous and discrete Freidlin-Wentzell type large-deviation principles (LDPs) are derived, characterizing the probability of rare events under vanishing noise:

• The continuous LDP rate function corresponds to the skeleton equation of the original SPDE.
• The discrete rate function for the MLEI converges to the continuous one in the sense of $\Gamma$-convergence as the timestep refines.

This provides a rigorous probabilistic foundation for the asymptotic behavior of trajectories computed via MLEI, ensuring that large-deviation probabilities are faithfully approximated in the discretization (Dai et al., 2022).

7. Context, Applications, and Significance

MLEI is naturally suited for the temporal discretization of models with fractional-order memory—such as those arising in viscoelasticity, anomalous diffusion, and heat propagation in heterogeneous media. The approach generalizes exponential Euler methods beyond integer-order systems, offering considerable improvements in temporal accuracy and permitting advanced analysis, such as sharp large-deviation estimates and the study of rare events. The framework supports both additive and fractionally integrated noise, accommodating a broad range of SPDEs of current research interest (Kovács et al., 2018, Dai et al., 2022).