Backward Euler Convolution Quadrature Method

• Backward Euler CQ is a numerical framework that approximates Riemann–Liouville and Caputo fractional derivatives by discretizing the Laplace operator using the backward Euler rule.
• It employs efficient algorithms such as FFT-based convolution and sum-of-exponentials compression to accelerate time-stepping while preserving first-order accuracy.
• The method is widely applied in simulating anomalous diffusion, wave propagation, dispersive media, and fractional variational integrators, ensuring robust stability and convergence.

The backward Euler convolution quadrature (CQ) method is a numerical framework for approximating fractional derivatives—specifically Riemann–Liouville and Caputo derivatives of order $0<\alpha\le1$—using time-stepping based on the backward Euler linear multistep scheme. It delivers a stable, first-order accurate discretization by systematically replacing the Laplace variable in the continuous operator with a symbol inherited from the underlying time integrator. CQ is extensively used in simulating anomalous diffusion, subdiffusion, wave propagation, electromagnetic phenomena in dispersive media, and fractional variational integrators (Hariz et al., 2024, &&&1&&&, Hassell et al., 2014, Dölz et al., 2020, López-Fernández et al., 2023, Melenk et al., 2024).

1. Algebraic and Operator Foundations of Backward Euler CQ

The core of the backward Euler CQ method is the mapping of continuous fractional operators to discrete analogues via the symbol $\delta(\zeta)$ of the backward Euler rule:

$\delta(\zeta) = 1 - \zeta,$

where $\zeta$ is the generating variable in the complex plane. The fractional power is raised as

$$(\delta(\zeta)/h)^\alpha = h^{-\alpha}(1-\zeta)^\alpha.$$

The approximation of the left-sided Riemann–Liouville derivative $D_{-}^\alpha y(t_n)$ or the Caputo derivative is then

$$\mathcal{D}_h^\alpha y_n := h^{-\alpha} \sum_{k=0}^{n} \omega_k^{(\alpha)} y_{n-k},$$

where the weights $\omega_k^{(\alpha)}$ are the binomial coefficients

$\omega_k^{(\alpha)} = (-1)^k \binom{\alpha}{k},$

with

$$\binom{\alpha}{k} = \frac{\alpha (\alpha-1) \cdots (\alpha-k+1)}{k!}.$$

For fractional integrals, the explicit formula becomes

$$\omega_n = h^{-\alpha} \frac{\Gamma(n+\alpha)}{\Gamma(\alpha)\Gamma(n+1)},$$

and the full discrete operator reads

$$(I_h^\alpha f)(t_n) = \sum_{j=0}^n \omega_{n-j} f(t_j),$$

where $h$ is the uniform time step and $t_n = n h$ (Liu et al., 2019).

2. Algorithmic Implementation and Computational Acceleration

Implementation requires precomputing the weights $\omega_k^{(\alpha)}$ via a three-term recurrence:

$$\omega_0^{(\alpha)} = 1,\quad \omega_k^{(\alpha)} = -\frac{\alpha-(k-1)}{k} \omega_{k-1}^{(\alpha)},\quad k\ge1.$$

Time-stepping at $n$ involves evaluating the convolution sum, costing $O(n)$ per step and $O(N^2)$ for $N$ total steps.

Several acceleration techniques are available:

• FFT-based convolution: Zero-padding the time series and weights, then performing a discrete Fourier transform accelerates convolution to $O(N\log N)$ (Hassell et al., 2014).
• Sum-of-exponentials compression: Using parameter integral representations for weights and Gauss–Jacobi quadrature, history variables can be updated via geometric recurrences at $O(\log N)$ cost per step, dramatically reducing work and storage (Sun et al., 2019, López-Fernández et al., 2023).
• Parsimonious CQ: Piecewise polynomial (Chebyshev) interpolation of the Laplace domain symbol enables the method to require only $O(\sqrt{N}\log N)$ evaluations of the Laplace-domain operator in some settings, notably in acoustic/electromagnetic wave boundary integral equations (Melenk et al., 2024).
• Oblivious convolution quadrature: Blocks history sums and uses contour quadrature in the Laplace domain for $O(\log N)$ work and $O(\log N)$ storage, retaining first-order accuracy.

3. Integration with Fractional PDEs and Variational Principles

Backward Euler CQ is commonly applied to fractional PDEs such as subdiffusion and fractional Fokker-Planck equations. In such cases, the discrete operator is substituted directly for the fractional derivative in the semi-discrete or fully discrete numerical scheme. For example, in a CQ-FEM discretization of a quasilinear subdiffusion equation, time-stepping updates combine CQ evaluation for the Caputo derivative and finite element assembly for the spatial part:

$$(\partial_h^\alpha U^n, \chi) + \int_\Omega D(x, t_n, U^{n-1}) \nabla U^n\cdot \nabla\chi\,dx = (f(x, t_n, U^{n-1}), \chi),$$

where $\partial_h^\alpha$ is the CQ operator (López-Fernández et al., 2023).

For fractional variational integrators, the CQ approximation replaces all appearances of fractional derivatives and fractional integrals in the discrete action, and the discrete Hamilton principle yields difference equations that couple standard variational updates with nonlocal fractional damping contributions. The global accuracy in time of such variational integrators is limited by the order of the underlying CQ (first-order for backward Euler) (Hariz et al., 2024).

4. Stability, Convergence, and Saturation Phenomena

Backward Euler CQ inherits the $A$-stability of the backward Euler integrator. Specifically, stability regions are characterized via the mapping of the spectrum of the sectorial operator to the symbol $\delta(\zeta)/h$, guaranteeing boundedness and damping of solutions with Re$\,\lambda<0$ (Liu et al., 2019). The method yields

$$\|D^\alpha y - \mathcal{D}_h^\alpha y\|_{\infty,[h,T]} = O(h),$$

for sufficiently smooth data, uniformly away from $t=0$.

Saturation effects arise when the data possesses singularities at $t=0$; if $y(t) \sim t^{\beta-1}$ near the origin, the local error is $O(h^{\min(1,\beta)})$ owing to the loss of regularity. Start-up corrections—modifying initial weights to match polynomial data—restore optimal global convergence even for nonsmooth initial data (Liu et al., 2019).

5. Applications: Wave Propagation, Diffusion, and Dispersive Media

Typical applications include

• Wave simulations via boundary integral equations: The CQ weights are constructed using operator-valued Laplace transforms (e.g., of single-layer or boundary potential operators), and discrete convolution solves are performed at each time step (Hassell et al., 2014, Betcke et al., 2016).
• Fractional Fokker–Planck and subdiffusion equations: The method allows efficient discretization even for solutions with low regularity (Sun et al., 2019, López-Fernández et al., 2023).
• Maxwell’s equations in dispersive media: The polarization term is discretized via CQ, preserving the energy-dissipation structure and yielding stability and first-order accuracy in the time step (Dölz et al., 2020).
• Fractional variational mechanics: CQ enables geometric integrators with variational structure for systems with fractional damping (Hariz et al., 2024).

6. Fast, Parsimonious, and Oblivious CQ Variants

Recent developments include fast and parsimonious CQ algorithms:

• Fast CQ compresses the historical memory term using sum-of-exponentials and history recursions (Sun et al., 2019).
• Parsimonious CQ applies piecewise Chebyshev interpolation in the Laplace-phase domain, reducing Laplace operator evaluations to $O(\sqrt{N}\log N)$ while maintaining convergence order, applicable to scattering and other wave problems (Melenk et al., 2024).
• Oblivious CQ blocks time-history into logarithmically sized intervals, using contour quadrature and running sums for minimal memory footprint (López-Fernández et al., 2023, Dölz et al., 2020).

7. Practical Considerations and Extensions

Implementation demands careful computation of convolution weights, consideration of start-up corrections for nonsmooth problems, and efficient data structures to manage historical sums. Quadrature error and time discretization error must be balanced, especially in overresolving strategies—where the Laplace-domain contour integral is discretized with more points than time steps to achieve exponential convergence in the aliasing error (Betcke et al., 2016).

A plausible implication is that, for memory-intensive simulations on large time intervals or with stiff kernels, employing fast, parsimonious, or oblivious CQ significantly reduces computational resource requirements without sacrificing stability or convergence order.

• "Fractional variational integrators based on convolution quadrature" (Hariz et al., 2024)
• "Fast algorithms for convolution quadrature of Riemann-Liouville fractional derivative" (Sun et al., 2019)
• "Convolution Quadrature for Wave Simulations" (Hassell et al., 2014)
• "The unified theory of shifted convolution quadrature for fractional calculus" (Liu et al., 2019)
• "Overresolving in the Laplace domain for convolution quadrature methods" (Betcke et al., 2016)
• "Convolution Quadrature for the quasilinear subdiffusion equation" (López-Fernández et al., 2023)
• "Parsimonious convolution quadrature" (Melenk et al., 2024)
• "A convolution quadrature method for Maxwell's equations in dispersive media" (Dölz et al., 2020)