Padé-Based ETD2RK-DS Formulation

• The paper introduces a Padé-based ETD2RK-DS formulation that replaces analytic exponentials with low-degree rational approximants to boost computational efficiency.
• It employs tensor slicing and dimension splitting to decouple high-dimensional reaction–diffusion problems, reducing per-step cost significantly.
• Empirical benchmarks and theoretical analysis confirm uniform stability and second-order time accuracy even for large-scale stiff semi-discrete systems.

The Padé-based ETD2RK-DS formulation denotes a class of second-order exponential time-differencing Runge–Kutta schemes for multidimensional reaction–diffusion equations, implemented with dimension (spatial) splitting and rational Padé approximations of operator exponentials and associated $\varphi$-functions. This approach enables stable, accurate, and highly efficient time integration for large-scale stiff semi-discrete systems arising from the spatial discretization of reaction–diffusion problems in two and three spatial dimensions, dramatically reducing computational cost relative to classical unsplit implicit solvers by exploiting both tensor (Kronecker) structure and spectral decompositions (Sarumi, 11 Jan 2026, &&&1&&&).

1. Padé Approximation for Matrix Exponentials and $\varphi$-Functions

The foundation of the method is the replacement of analytic exponential and $\varphi$-functions by low-degree Padé rational approximants. For an argument $z$, the $\varphi$-functions required by exponential integrators are: $$\varphi_0(z)=e^z, \quad \varphi_1(z)=\frac{e^z-1}{z}, \quad \varphi_2(z) = \frac{e^z-z-1}{z^2}, \;\dots$$ These are approximated, for example, by the $(0,2)$ Padé approximant for $\varphi_0(z)$: $$P_{0,2}(z)=\frac{2}{z^2+2z+2}\quad \text{with}\quad e^z - P_{0,2}(z) = O(z^3)\quad \text{as}\ z\to 0.$$ For practical recursion and linear solvers, three rational functions are defined: $$Q_1(z) = P_{0,2}(z),\quad Q_2(z) = \frac{z+2}{z^2+2z+2},\quad Q_3(z) = \frac{z+1}{z^2+2z+2},$$ with partial fraction decompositions at poles $s = -1 + i$, $\bar{s} = -1 - i$ and residues $r_1=-i$, $r_2=(1-i)/2$, $r_3=1/2$. The rational actions $Q_\ell(z)g$ correspond to solutions of shifted linear systems of the form $(z-sI)v = r_\ell g$, with final results given as $2\Re(v)$ (Sarumi, 11 Jan 2026).

2. Dimension-Split ETD2RK Scheme

After discretizing space, the semi-discrete initial value problem for $u\in\mathbb{R}^{m^d}$ is

$$u'(t) + A_hu(t) = f(t,u(t)), \qquad A_h = A_1 + A_2 \quad (\text{in 2D}),$$

where $A_1$ and $A_2$ are Kronecker-structured matrices representing directional Laplacians. The second-order, dimension-split ETD2RK-DS time-stepping advances the solution via

$$\begin{aligned} W^{n+1} &= e^{-\tau A_1}e^{-\tau A_2}U^n + A_1^{-1}(I-e^{-\tau A_1})e^{-\tau A_2}f(t_n,U^n), \ U^{n+1} &= W^{n+1} + (A_1^{-1}\tau - A_1^{-2} + A_1^{-2}e^{-\tau A_1})\frac{f(t_{n+1},W^{n+1}) - e^{-\tau A_2}f(t_n,U^n)}{\tau}, \end{aligned}$$

with each matrix exponential and $\varphi$-function product replaced by an equivalent Padé rational operator (Sarumi, 11 Jan 2026). The fully discrete, Padé-based forms read: $$\begin{aligned} \widehat W^{n+1} &= Q_1(\tau A_1) Q_1(\tau A_2) \widehat U^n + \tau Q_2(\tau A_1) Q_1(\tau A_2)f(t_n,\widehat U^n), \ \widehat U^{n+1} &= \widehat W^{n+1} + \tau Q_3(\tau A_1) \left[ f(t_{n+1},\widehat W^{n+1}) - Q_1(\tau A_2)f(t_n,\widehat U^n) \right]. \end{aligned}$$ Each use of $Q_\ell(\tau A_\nu)$ acts via solutions of complex-shifted linear systems as above, with tensor structure enabling reduction to independent 1D solves.

3. Stability, Convergence, and Theoretical Guarantees

For locally Lipschitz $f$ in $u$ and $A_i$ satisfying $\|e^{-\tau A_i}\| \leq 1$, uniform stability is established: $\|U^n\| \leq M_*,\quad 0 \leq n \leq N,$ where $M_*$ is independent of $\tau$. Second-order accuracy in time is proven for sufficiently regular exact solutions: $\|u(t_n) - U^n\| = O(\tau^2),\quad n\tau \leq T,$ by an analysis leveraging the mild solution representation, Lipschitz bounds, and careful control of the Padé and $\varphi$-function remainder terms (Sarumi, 11 Jan 2026).

The use of real or complex Padé rational approximations preserves $L$-stability, contingent on suitable pole-residue parametrizations. For schemes using real distinct poles (RDP), the scalar stability function

$$R(z) = \frac{9}{1+\frac{1}{3}z} - \frac{8}{1+\frac{1}{4}z} + z\left(\frac{4}{1+\frac{1}{3}z} - \frac{3}{1+\frac{1}{4}z}\right) + z^2\left(\frac{2}{1+\frac{1}{3}z} - \frac{1}{1+\frac{1}{4}z}\right)$$

satisfies $\lim_{\Re z \to -\infty}R(z)=0$, confirming the strong damping of stiff modes (Asante-Asamani et al., 2020).

4. Implementation: Tensor Slicing and Sylvester Equation Solvers

The Kronecker sum structure of $A_1,\ldots,A_d$ allows dimension splitting and explicit tensor rearrangement: solutions of $(\tau A_\nu - sI)v = r_\ell g$ are broken into $O(m^{d-1})$ independent tridiagonal problems of size $m$ each, for a total application cost of $O(m^d)$ per direction.

For instance, in 2D:

• For $A_1=I_y\otimes A_x$, slice $g$ into $g^{[j]}\in\mathbb{R}^{m_x-1}$ and solve

$$(\tau A_x - sI)v^{[j]} = r_\ell g^{[j]},\quad j=1,\dots,m_y-1.$$

• For $A_2=A_y\otimes I_x$, index swapping produces $\tilde g^{[i]}$, similarly solved.

Pseudocode for a 2D $Q_\ell(\tau A_2)g$ operation:

for i = 1 to m_x-1: extract tilde g^{[i]} solve (τ A_y – s I) · tilde v^{[i]} = r_ℓ · tilde g^{[i]} scatter tilde v^{[i]} back into v return 2 · Re(v)

This slicing reduces the 2D per-step cost from $O(m^3)$ (for full unsplit LU) to $O(m^2)$; similarly $O(m^5)$ to $O(m^3)$ in 3D (Sarumi, 11 Jan 2026).

Each 1D shifted system may be recast as a Sylvester equation or diagonalized explicitly if $A_x = P_x D_x P_x^T$. Then, with stored eigendecomposition,

$$(\tau D_x - sI)\bar v = r_\ell P_x^T g,\quad v = P_x \bar v,$$

where all right-hand sides are processed via two matrix–vector multiplies plus a Hadamard division. These operations are $O(m^2)$ in 2D, $O(m^3)$ in 3D (Sarumi, 11 Jan 2026).

5. Empirical Performance and Benchmark Problems

Numerical experiments validate the method's accuracy and cost reductions. For the Allen–Cahn equation and FitzHugh–Nagumo systems in 2D and 3D, empirical error rates confirm $O(\tau^2)$ time accuracy—e.g., $E(N)\approx C N^{-2}$.

Performance benchmarks highlight substantial efficiency gains:

Solver	Cost per Step (2D)	CPU Time ($m=512$)	Scaling
Sylvester (Padé)	$O(m^2)$	$0.8$–$13.8$ s	Mild increase
LU-slice (Padé)	$O(m^2)$	$1.1$–$18$ s	Mild increase
Full LU (unsplit)	$O(m^3)$	Tens of seconds	Steep increase

The Sylvester-based solver incurs only moderate additional cost when switching to higher-order Padé, in contrast to LU-based approaches which approximately double in cost due to more linear solves (Sarumi, 11 Jan 2026).

6. Comparison with Alternative ETD2RK-DS Formulations

A related dimensionally-split ETD2RK-DS with real-pole rational approximants (RDP) has been developed independently (Asante-Asamani et al., 2020). This uses second-order rational approximants of the form: $$e^{-z} \approx r(z)=\frac{1-\frac{5}{12}z}{(1+\frac{1}{3}z)(1+\frac{1}{4}z)} = \frac{9}{1+\frac{1}{3}z} - \frac{8}{1+\frac{1}{4}z},$$ yielding sparse real-shifted linear solves and favorable $L$-stability properties. The implementation benefits from Thomas-type (tridiagonal) solvers ($O(m)$), and is naturally parallelizable over directions.

Empirical studies report that such dimensionally split Padé/RDP-based ETD2RK schemes outperform classical IMEX and ETD–Crank–Nicolson solvers in both accuracy and computational efficiency, with observed efficiency orders of convergence $\mathrm{EOC}\approx2$. Large-scale 3D problems (e.g., $200^3$ nodes) are solved in practical wall-clock time on standard desktops (Asante-Asamani et al., 2020).

7. Scope, Generalizations, and Practical Considerations

The Padé-based ETD2RK-DS approach is applicable to semilinear parabolic PDEs on tensor-product domains, with Dirichlet, Neumann, or periodic boundary conditions. The method assumes Kronecker-structured (separable) discretizations, with spectral (eigendecomposition-based), or tridiagonal solvers, applicable to both coupled and single-species systems. Advantages include:

• Uniform stability and second-order convergence for mildly restricted nonlinearities (locally Lipschitz $f$; practical stability observed in locally bounded, non-Lipschitz $f$).
• Scalability: per-step cost reduction from $O(m^3)$ to $O(m^2)$ in 2D (and $O(m^5)$ to $O(m^3)$ in 3D) via explicit tensor-slicing and spectral solvers.
• Capability to use higher-order rational approximants with only mild cost increase.
• Compatibility with parallel and spectral methods (e.g., FFT for periodic cases).

A plausible implication is that for high-dimensional, large-scale stiff diffusive–reactive problems where operator splitting and effective exploitation of structure is possible, Padé-based ETD2RK-DS formulations provide a near-optimal practical strategy among exponential integrator schemes (Sarumi, 11 Jan 2026, Asante-Asamani et al., 2020).