Two-Level Strang Splitting for SFAC Equations

• The paper introduces two-level Strang splitting as a second-order accurate method that decomposes Toeplitz matrices into circulant and skew-circulant parts.
• It leverages the Fast Fourier Transform to reduce computational complexity from O(N²) to O(N log N) while maintaining the discrete maximum principle.
• Error analysis confirms optimal convergence with O(τ² + h²) rates, ensuring practical scalability in 2D and 3D simulations.

Two-level Strang splitting is a numerically efficient, second-order accurate approach for solving multi-dimensional spatial fractional Allen–Cahn (SFAC) equations, which are characterized by the presence of nonlocal Riesz fractional derivatives. The method addresses the principal computational hurdle in these problems: the expensive calculation of Toeplitz matrix exponentials arising after spatial semi-discretization. By decomposing the Toeplitz operators into circulant and skew-circulant components and embedding the entire time evolution in a nested Strang splitting framework, the method enables fast implementation via the Fast Fourier Transform (FFT), ensuring both rigorous stability and optimal complexity (Cai et al., 2022).

1. Formulation and Spatial Discretization

Consider the $d$-dimensional SFAC equation:

$$u_t = \varepsilon^2 \mathcal{L}_{x^d}^\alpha u + u - u^3, \quad x \in [a,b]^d, \quad u|_{\partial\Omega} = 0,$$

where $\mathcal{L}_{x^d}^\alpha$ denotes the sum of one-dimensional Riesz derivatives with orders $\alpha_\ell \in (1,2)$. Semi-discretization in space proceeds via uniform grids in each coordinate. Utilizing the second-order weighted-shifted Grünwald approximation for the Riemann–Liouville derivatives,

$_{a}D_x^\alpha u \approx h^{-\alpha} \sum_{k=0}^{i+1} \omega_k^{(\alpha)} u_{i-k+1}, \qquad _{x}D_b^\alpha u \approx h^{-\alpha} \sum_{k=0}^{m-i+1} \omega_k^{(\alpha)} u_{i+k-1},$

yields discrete weights $\omega_k^{(\alpha)}$ defined recursively. The assembled system over all directions gives rise to $N$ coupled ODEs,

$$\frac{d\mathbf{u}}{dt} = A\mathbf{u} + f(\mathbf{u}), \quad f(\mathbf{u}) = \mathbf{u} - \mathbf{u}^3,$$

where $A$ is a sum of Kronecker products with symmetric Toeplitz matrices $B_\alpha$ constructed from the discretization (Cai et al., 2022).

2. Toeplitz Matrix Splitting: Circulant and Skew-Circulant Decomposition

Each Toeplitz matrix $B_\alpha$ is uniquely decomposed as

$B_\alpha = C_\alpha + S_\alpha,$

with $C_\alpha$ circulant and $S_\alpha$ skew-circulant. The explicit formulas for the diagonals ensure that this splitting holds componentwise for all $d$ directions. Consequently, the full spatial operator $A$ is partitioned as

$$A = C + S, \qquad C = \sum_{\ell=1}^d I \otimes \cdots \otimes C_{\alpha_\ell} \otimes \cdots \otimes I, \quad S = \sum_{\ell=1}^d I \otimes \cdots \otimes S_{\alpha_\ell} \otimes \cdots \otimes I.$$

This decomposition is pivotal for the subsequent algorithmic acceleration (Cai et al., 2022).

3. Nested (Two-Level) Strang Splitting Scheme

The two-level Strang splitting combines time integration with an additional split of the linear evolution. The outer (standard) Strang splitting applies to the semi-discrete ODE:

$$u^{n+1} = \mathcal{J}^Q_{\tau/2} \mathcal{J}^P_\tau \mathcal{J}^Q_{\tau/2} u^n,$$

where $\mathcal{J}^Q_\tau$ is the exact nonlinear propagator for the local term $u - u^3$. The core computational obstacle is $\mathcal{J}^P_\tau = \exp(\tau A)$, which involves the full Toeplitz matrix exponential.

To reduce complexity, an inner Strang split is performed on $A = C + S$:

$$e^{\tau A} = e^{\frac{\tau}{2} C} e^{\tau S} e^{\frac{\tau}{2} C} + O(\tau^3).$$

The resulting full two-level Strang splitting update is

$$\hat{u}^{n+1} = \mathcal{J}^Q_{\tau/2}\, \widehat{\mathcal{J}}^P_\tau\, \mathcal{J}^Q_{\tau/2} \hat{u}^n, \qquad \widehat{\mathcal{J}}^P_\tau = e^{\frac{\tau}{2} C} e^{\tau S} e^{\frac{\tau}{2} C}.$$

Both temporal and inner linear splitting individually contribute $O(\tau^3)$ local truncation error, yielding global $O(\tau^2)$ accuracy under bounded commutator conditions (Cai et al., 2022).

4. Fast Fourier Transform-Based Implementation

Since both $C_\alpha$ and $S_\alpha$ are diagonalizable by the discrete Fourier matrix, their exponentials act via diagonal scaling in the frequency domain. The exponential actions are computed as:

$$e^{(\tau/2)C_\alpha} v = F^* \left( e^{(\tau/2)\Lambda_C} \cdot (Fv) \right), \qquad e^{\tau S_\alpha} v = \Psi^* F^* \left( e^{\tau \Lambda_S} \cdot (F(\Psi v)) \right),$$

where $F$ is the discrete Fourier transform and $\Psi$ is a "twist" diagonal matrix for skew-circulant structure. Each FFT or inverse FFT in $d$ dimensions requires $O(N \log N)$ operations, and all diagonalizations are independent and explicit. Overall, this enables each timestep of the full scheme to be executed with $O(d N \log N)$ complexity, a substantial reduction compared to $O(N^2)$ for direct Toeplitz exponentiation (Cai et al., 2022).

5. Stability via the Discrete Maximum Principle

A defining property of the two-level Strang splitting is its unconditional preservation of a discrete maximum principle:

$$\| \hat{u}^n \|_\infty \leq 1 \quad \forall n, \quad \text{provided} \quad \|u^0\|_\infty \leq 1,$$

with no requirement for restrictive time-step (CFL) conditions. This follows from the diagonal dominance and negativity of the diagonal entries in both $C_\alpha$ and $S_\alpha$, combined with the nature of the nonlinear propagator, each individually reducing or preserving the $\infty$-norm at every sub-step (Cai et al., 2022).

6. Error Analysis and Convergence

Assuming $f(u)$ is $C^2$, the ODE solution is in $C^2([0,T])$, and that all relevant commutators remain bounded, the scheme achieves second-order convergence in both time and space. The convergence theorem asserts

$$\| u_{\rm exact}(t_n) - \hat{u}^n \|_\infty \leq C (h^2 + \tau^2), \quad n\tau \leq T,$$

for some $C$ independent of discretization parameters, given $u_{\rm exact} \in C^5(\Omega) \times C^2([0,T])$ (Cai et al., 2022).

7. Computational Efficiency and Numerical Validation

The two-level Strang splitting method provides a dramatic reduction in computational cost for large-scale SFAC simulations. In numerical tests across 2D and 3D domains—with up to several million degrees of freedom—CPU time per step is reduced from tens of seconds (for direct Toeplitz exponentiation) to sub-second scales with the FFT-based method. Crucially, the schemes preserve the discrete maximum principle and monotonically decrease energy, while empirical convergence matches the theoretical $O(\tau^2 + h^2)$ rate.

Representative results from (Cai et al., 2022) include:

τ	$L^\infty$ Error	Order	CPU (s) (2D)
1/100	$1.28 \times 10^{-7}$	—	0.48
1/200	$3.21 \times 10^{-8}$	2.00	0.90
1/400	$8.02 \times 10^{-9}$	2.00	1.83

τ	$L^\infty$ Error (3D)	Order	CPU (s)
1/20	$1.99 \times 10^{-6}$	—	0.30
1/40	$4.96 \times 10^{-7}$	2.00	0.59
1/80	$1.24 \times 10^{-7}$	2.00	1.20

All observations confirm theoretical second-order accuracy and high efficiency (Cai et al., 2022).

---

Two-level Strang splitting establishes an effective framework for discretizing and simulating high-dimensional fractional phase-field models, coupling optimal stability, accuracy, and computational scalability (Cai et al., 2022).