Exponential-Trapezoidal Discretization

• Exponential-trapezoidal discretization is a numerical method that blends exponential integrator concepts with the trapezoidal rule to achieve high precision.
• It yields A-stable, second-order schemes with rigorous error control and rapid convergence for stiff ODEs, PDEs, and integro-differential equations.
• The approach extends to Fourier and operator integrals using Padé approximations and exponential transforms for efficient, high-order integration.

Exponential-trapezoidal discretization encompasses a family of numerical techniques that combine the classical trapezoidal rule with exponential weighting or change-of-variable transformations to achieve high efficiency and accuracy in the integration of differential equations, operator functions, and oscillatory integrals. These methods are especially relevant for stiff ODEs/PDEs, semilinear integro-differential equations, and quadrature involving Fourier or operator-theoretic integrals. Central to these approaches is the synergy of exponential integrator ideas (exact or rationally approximated action of matrix exponentials or associated semigroups) with the trapezoidal rule’s time-centering, yielding schemes that are A-stable, of at least second order, and that often admit rigorous error control and rapid convergence.

1. Exponential-Trapezoidal Method for Evolution Equations

The exponential-trapezoidal method discretizes time evolution of (semi)linear problems of the form $u'(t) = Au(t) + g(t,u(t))$ by making use of the variation-of-constants (Duhamel) formula. This formula expresses the solution at $t_{n+1}=t_n+h$ as:

$$u(t_{n+1}) = e^{hA}u(t_n) + \int_0^h e^{(h-\tau)A}g(t_n+\tau,u(t_n+\tau))\,d\tau.$$

The integral term is replaced by a trapezoidal approximation, but with the exponential weighting exactly preserved:

$$\int_0^h e^{(h-\tau)A}g(t_n+\tau,u(t_n+\tau))\,d\tau \approx \frac{h}{2}\big[ e^{hA}g(t_n,u(t_n)) + g(t_{n+1},u(t_{n+1})) \big].$$

This yields the implicit exponential-trapezoidal scheme:

$$U_{n+1} = e^{hA}U_n + \frac{h}{2} \left[ e^{hA}g(t_n, U_n) + g(t_{n+1}, U_{n+1}) \right].$$

A Taylor expansion reveals local truncation error $O(h^3)$, so the method is second order in time. Each step is implicit, but fixed-point iteration (Picard) rapidly converges due to the stiffness being addressed by the exponential factors (Ostermann et al., 2024).

2. Exponential-Trapezoidal and Padé-Based Time Integration

For linear systems, especially those arising in structural dynamics and wave propagation, exponential–trapezoidal schemes can be systematically generalized using Padé approximations to the matrix exponential. Consider $\dot{x}(t)=Ax(t)+f(t)$; the homogeneous solution part involves $e^{hA}x(t)$. The diagonal Padé approximant $[M/M]$ of $e^z$,

$e^z \approx \frac{P_M(z)}{Q_M(z)},$

can be factorized to admit a partial fraction decomposition, reducing the action of $R(hA)$ to a sequence of $M$ sparse linear solves. For $M=1$, the method coincides with the classical trapezoidal/Newmark scheme (second-order, unconditionally A-stable). Higher $M$ yields $2M$-th order methods, requiring $M$ linear solves per step (Song et al., 2021).

The generic update for time step $n\to n+1$ has the structure:

$x_{n+1} = \sum_{j=1}^M \beta_j y^{(j)}_{n+1},$

where each $y^{(j)}_{n+1}$ satisfies

$$(1-\alpha_j h A) y^{(j)}_{n+1} = x_n + \sum_{k=0}^M C_k^{(j)} f_n^{(k)},$$

with $\alpha_j$ and $\beta_j$ defined by Padé theory. When applied to second-order systems (e.g., with mass, damping, and stiffness matrices), this leads to $M$ decoupled sparse linear systems for displacement updates per timestep.

Key attributes:

• For $M=1$: equivalence with trapezoidal/Newmark-$\gamma=1/2$, $\beta=1/4$.
• Global order: $2M$ for Padé($M$, $M$).
• A-stability for all $M$; not L-stable (no artificial high-frequency numerical damping).

3. Exponential-Trapezoidal Quadrature for Oscillatory and Operator Integrals

For integrals on the half-line or unbounded domains, such as Fourier type integrals and operator fractional powers, exponential–trapezoidal discretizations combine variable transformations inducing (double-)exponential decay with the classical trapezoidal rule. For the cosine transform,

$$F^{(c)}(\omega) = \int_0^\infty f(x)\cos(\omega x)\,dx,$$

Ooura–Mori transformations of the form $x = \tau/\omega\,\phi_i(\xi)$, with

$$\phi_1(\xi) = \frac{\xi}{1-\exp(-2\pi\sinh\xi)}, \quad \phi_2(\xi) = \frac{\xi}{1-\exp(-2\xi-\alpha(1-e^{-\xi})-\beta(e^{\xi}-1))}, \; 0<\alpha<\beta<1,$$

yield rapid decay and mesh compression, resulting in quadrature sums with doubly exponential convergence. This technique enables a priori mesh and truncation parameter selection, and error decomposition (discretization, left/right truncation) with rigorous strip-analytic or residue-based bounds (Denich et al., 2023).

For operator functions such as fractional powers,

$$A^{-\alpha} = \frac{2\sin(\alpha\pi)}{\pi} \int_0^\infty t^{2\alpha-1} (I+t^2A)^{-1} dt,$$

single-exponential (SE: $t = e^y$) and double-exponential (DE: $t = \exp(\pi/2 \sinh y)$) transformations, paired with the trapezoidal rule, yield approximations with nearly $\exp(-c\sqrt{n})$ and $\exp(-c\sqrt{n/\log n})$ convergence rates, respectively, where $n$ is the number of quadrature nodes (Aceto et al., 2021).

4. Error Analysis and Convergence Properties

For the ODE/PDE time-stepping case, the fully discrete exponential–trapezoidal scheme admits the error bound (in an abstract Hilbert space $H$, with spectral Galerkin approximation using $N$ eigenmodes and time step $h$):

$$\|u(t_n) - U_n^N \|_V \leq C\left[t_n^{-\alpha\rho}\lambda_{N+1}^{-\alpha-\beta} + h\,t_n^{-\nu\rho}\lambda_{N+1}^{-\beta} + \lambda_{N+1}^{\nu-\alpha-\gamma} + h^{2-(\nu-\eta)\rho}\sup_{t\leq T} \|A^\eta g'(t)\| + h^2 \sup_{t\leq T}\|A^{-\delta}g''(t)\| \right]$$

with regularity exponents $\beta, \gamma, \eta, \delta$ as in (Ostermann et al., 2024). For fixed spatial discretization ($N$ large, spectral errors negligible), temporal error is $O(h^2)$.

In the quadrature context, both SE and DE exponential–trapezoidal rules achieve exponential convergence depending on analyticity strip width. For the SE scheme,

$$\left\|A^{-\alpha} - Q_{M,N,h}^{SE}(A)\right\| \leq C \exp\left( -\pi\sqrt{2\alpha(1-\alpha)}\sqrt{n} \right),$$

with $n = M+N+1$, and for DE,

$$\|A^{-\alpha} - Q_{n,h}^{DE}(A)\| \lesssim \overline K_\alpha\exp\left( -c\sqrt{n/\log n}\right).$$

For Fourier-type integrals, error decomposes into discretization and left/right truncation parts, with explicit parameter-choice algorithms available for guaranteed accuracy given tolerances (Denich et al., 2023, Aceto et al., 2021).

5. Implementation Strategies

Efficient evaluation of the exponential operator $e^{hA}$ for time-stepping may be realized by:

• Diagonalizing $A$ (spectral methods) when the eigenbasis is available.
• Krylov subspace approximations (Arnoldi/Lanczos), with per-vector cost $O(m\,\text{cost}(A \cdot \text{vector}))$.
• Rational and polynomial interpolation methods (Leja points, contour integrals).

For Padé-based schemes, all requisite factorizations can be precomputed for fixed $M$ to accelerate repeated solves. For operator quadrature, computation reduces to a set of resolvent solves at positive scalar or operator shifts.

The per-step computational cost for Padé-based time integrators scales with $M$: for $M=1$ it matches the Newmark method, while higher-order versions require proportionally more solves but can allow larger timesteps for the same accuracy (Song et al., 2021).

For Fourier or operator quadrature, the complexity scales with the number of quadrature nodes ($O(N)$ function or resolvent evaluations), with automatic parameter selection algorithms enabling black-box use (Denich et al., 2023).

6. Numerical Experiments and Practical Results

Benchmark tests of exponential–trapezoidal integrators include:

• Heat equations with memory (using Riesz and exponential kernels) and nonlinearity, showing second-order time convergence unaffected by step size restriction (Ostermann et al., 2024).
• Structural dynamic systems up to nearly $10^6$ DOF; higher-order Padé-based exponential–trapezoidal methods observed to outperform conventional schemes for long-time integration and high accuracy (Song et al., 2021).
• Fractional operator powers on high-contrast spectra, with exponential or super-exponential decay of errors versus number of resolvent solves, confirming theoretical rates (Aceto et al., 2021).
• Oscillatory quadrature with analytically known and singular integrands, demonstrating accuracy in line with predicted bounds with automatic parameter selection (Denich et al., 2023).

7. Scope of Applicability and Future Prospects

Exponential–trapezoidal discretization schemes generalize classical time integration and quadrature, providing A-stable and arbitrarily high-order accurate solvers for stiff ODE/PDE evolution, fractional operator computation, and highly oscillatory or improper integrals. The methods are particularly impactful where robust error control, uncondi­tional stability, or stiff source handling is critical (e.g., semilinear evolutionary PDEs, large-scale structural dynamics, and fractional diffusion).

A plausible implication is continued extension to nonlinear and multiphysics PDEs, improved operator-function quadrature, and rigorous adaptive parameter selection frameworks. The non-L-stability of standard Padé-based variants indicates continuing need for hybrid or low-damping alternatives in contexts requiring suppression of nonphysical high-frequency modes (Song et al., 2021, Ostermann et al., 2024).

---

References:

• (Ostermann et al., 2024): The exponential trapezoidal method for semilinear integro-differential equations
• (Song et al., 2021): High-order implicit time integration scheme based on Padé expansions
• (Denich et al., 2023): Some notes on the trapezoidal rule for Fourier type integrals
• (Aceto et al., 2021): Exponentially convergent trapezoidal rules to approximate fractional powers of operators