Second-Order Fast Finite Difference Schemes

• A second-order fast finite difference scheme is a discretization framework that achieves O(h²+τ²) accuracy using implicit time discretizations and second-order spatial stencils.
• It employs optimized FFT-accelerated preconditioned Krylov methods to solve structured matrices, reducing computational complexity to nearly O(N log N).
• These schemes are robust for fractional, distributed-order, and variable-coefficient PDEs, ensuring unconditional stability and reliable convergence even with non-smooth initial data.

A second-order fast finite difference scheme refers to a discretization and linear algebraic framework for time-dependent PDEs in which the local truncation error is $O(h^2 + \tau^2)$ (grid spacing and timestep), and where the resulting linear systems at each timestep are solved by nearly optimal algorithms—typically leveraging structured matrices (e.g., Toeplitz, circulant, Sylvester forms) and FFT-accelerated preconditioned Krylov subspace methods to ensure $O(N \log N)$ or comparably efficient complexity. This paradigm is particularly impactful for non-integer-order (fractional) diffusion, distributed-order, and variable-coefficient PDEs arising in anomalous transport, wave, and subdiffusion models, where explicit schemes or naive direct solvers are computationally prohibitive. Robust analysis covers accuracy, unconditional stability, and convergence, even in challenging settings such as distributed-order derivatives, Riesz space operators, and non-smooth initial data.

1. Core Principles and Problem Scope

Second-order fast finite difference schemes combine high-order accuracy with optimal or quasi-optimal solution complexity. The distinguishing attributes are:

• Second-order approximation in both time and space (e.g., $O(h^2 + \tau^2)$ global error in maximum or energy norm), often extended to distributed or variable-order terms to yield $O(\Delta\alpha^2)$ accuracy when discretizing integrals over order parameters (Jian et al., 2020, Jian et al., 2018, Zhao et al., 2017, Hao et al., 15 Jun 2024).
• Fully implicit time discretizations—typically Crank–Nicolson or specially constructed weighted and shifted Grünwald/Alikhanov interpolatory schemes for fractional/integral/time distributed-order operators.
• Fast linear algebra based on the exploitation of structured coefficient matrices (Toeplitz, block-Toeplitz-toeplitz-blocks, Sylvester, or banded) and the use of FFT-based preconditioned Krylov subspace solvers.
• Wide applicability: includes classical, Caputo, Riesz, distributed-order, and variable-order fractional equations; anisotropic and multidimensional settings; variable coefficients and weak singularities at initial time (Huang et al., 17 Dec 2025, Zhang et al., 2021, Jian et al., 2020).

These methods address high-cost challenges intrinsic to nonlocal operators in fractional PDEs, which otherwise yield dense matrices and require $O(N^2)$–$O(N^3)$ computation per timestep if solved naïvely.

2. Discretization Techniques

Time Discretization

• For fractional Caputo or distributed-order derivatives, weighted and shifted Grünwald formulas or Alikhanov ($L2$–$1_\sigma$) interpolatory quadratures deliver second-order accuracy uniformly over the fractional order. For the distributed-order integral $\int_{a}^{b} \omega(\alpha)\,^C D_t^\alpha u(t)\,d\alpha$, a composite trapezoidal quadrature reduces to a sum of uniform-fractional subproblems with $O(\Delta\alpha^2)$ error (Jian et al., 2020, Jian et al., 2018).
• For nonsmooth initial data (e.g., weak singularities), graded time meshes $t_n = T (n/N)^r$ restore second-order temporal accuracy under regularity restrictions, in conjunction with second-order Lagrange interpolation for local and nonlocal contributions (Huang et al., 17 Dec 2025).
• In variable/fractional order cases, exponential sum approximations (ESA) for the convolution kernel enable storage and computationally efficient second-order schemes. The fast $L2$–$1_\sigma$ (FL2–$1_\sigma$) approach reduces the computation to $O(\log N)$ per time step/component for variable-order kernels while preserving accuracy and stability (Zhang et al., 2021).

Space Discretization

• Riesz/Caputo/variable-order fractional derivatives are discretized using second-order (or fourth-order, if needed) Grünwald–Letnikov, weighted-shifted Grünwald, or centered difference stencils. For the Riesz derivative of order $\beta$, discretizations are second-order consistent in $h$ via symmetric Toeplitz stencils, for example,

$$\frac{\partial^\beta u}{\partial|x|^\beta}(x_i) = -h^{-\beta} \sum_{k} \hat{g}_k^{(\beta)} u(x_{i-k}) + O(h^2)$$

with explicit coefficient formulas (Jian et al., 2020, Lin et al., 2019).

• For variable-order or multi-dimensional/anisotropic fractional Laplacians, semi-discrete Fourier analysis leads to stencils whose weights are obtained by numerical quadrature and Chebyshev interpolation in the order parameter. The resulting schemes retain $O(h^2)$ local error (Hao et al., 15 Jun 2024).

3. Algebraic Structure and Fast Solvers

Matrix Structure and Preconditioning

• The system matrices are typically structured:
  • Toeplitz or symmetric Toeplitz in 1D (fractional derivatives).
  • Block-Toeplitz-with-Toeplitz-blocks (BTTB) in 2D/3D.
  • Sylvester operator form for coupled spatial directions (Jian et al., 2020).
• Preconditioners:
  • Classical circulant or block circulant preconditioners (R. Chan type) for Toeplitz/BTTB matrices, diagonalizable via the FFT, producing spectra clustered around unity and enabling rapid convergence (Jian et al., 2018, Lin et al., 2019).
  • Truncated banded Toeplitz approximations for variable/coupled coefficient problems.
• For distributed-order and variable-order schemes, the Chebyshev interpolation or sum-of-exponentials approximations allow rewriting the operator as a sum of constant-order convolution stencils with local coefficient multipliers, retaining the fast FFT-based application (Hao et al., 15 Jun 2024, Zhang et al., 2021).

Iterative Solution and Complexity

• Solving the resulting linear system at each step employs preconditioned CG or BiCGSTAB (1D) or global CG in Frobenius/Hilbert-Schmidt inner product (2D Sylvester) (Jian et al., 2020).
• Each Krylov iteration is $O(M\log M)$ for $M$ grid points.
• Empirically, the number of iterations is mesh-independent ($O(1)$) due to eigenvalue clustering.
• Total time per step is $O(M\log M)$ (1D) or $O(M_1M_2\log(M_1M_2))$ (2D).
• In explicit schemes (e.g., Lax–Wendroff), operator-matching and explicit moment-matching yield classical $O(N)$ fast update per step when the CFL condition is satisfied (Chin, 2013).

4. Stability, Convergence, and Robustness

• All referenced schemes are rigorously proven unconditionally stable in appropriate discrete energy, $L_2$, or maximum norms, even for nonlocal/fractional operators (Jian et al., 2020, Jian et al., 2018, Lin et al., 2019, Zhao et al., 2017, Liu et al., 2017). The fundamental strategy is multiplying the discrete scheme by a suitable test function (often the solution increment), leveraging symmetric/negative definiteness of space operators, and applying discrete Grönwall's inequalities.
• The convergence rate is $O(h^2+\tau^2)$, or $O(h^2 + \tau^2 + \Delta\alpha^2)$ for distributed order, under solution regularity matching the discretization (Jian et al., 2020, Jian et al., 2018, Hao et al., 15 Jun 2024).
• Weak regularity (e.g., non-smooth initial data) is typically handled by mesh grading; provided the grading parameter is sufficiently large ($r \ge \max\{3,\,2/(\delta-1)\}$ where $\delta$ is the solution regularity exponent), second-order convergence persists (Huang et al., 17 Dec 2025).

5. Applications and Numerical Performance

• These methods are applied to distributed-order and Riesz/fractional PDEs, advection-diffusion equations with variable coefficients, tempered time-fractional equations with initial singularity, and multi-dimensional variable-order fractional Laplacians (Jian et al., 2020, Huang et al., 17 Dec 2025, Hao et al., 15 Jun 2024).
• Numerical experiments (e.g., for $u(x,t)=t^4 x^3(1-x)^3$ or multidimensional benchmarks) verify second-order rates and show that FFT-based preconditioned Krylov solvers offer per-step CPU times up to $5$–$8\times$ faster than direct LU/Cholesky on large-scale problems, with per-iteration counts nearly independent of grid resolution (Jian et al., 2020, Jian et al., 2018, Lin et al., 2019).
• For variable-order/tempered models, the fast ESA schemes maintain full accuracy and drastically reduce both computational work and memory requirements, scaling as $O(N \log N)$ or better (Zhang et al., 2021, Huang et al., 17 Dec 2025).
• In multidimensional settings with non-rectangular geometry, domain embedding and appropriate enforcement of Dirichlet boundary conditions ensure maximal regularity is exploited for second-order convergence (Hao et al., 15 Jun 2024).

Reference	Problem Type	Fast Solver	Accuracy
(Jian et al., 2020)	Time distributed-order + Riesz	GSF+PCG, Global PCG (2D)	$O(h^2 + \tau^2 + \Delta\alpha^2)$
(Jian et al., 2018)	Distributed-order fractional diff.	Circulant PCG (1D/BTTB)	$O(h^2 + \tau^2 + \Delta\alpha^2)$
(Lin et al., 2019)	OSFDE, variable coefficient	Toeplitz PCG	$O(h^2 + \tau^2)$
(Chen et al., 2013)	2D space-fractional convection-diff	ADI splitting	$O(h^2 + \tau^2)$
(Huang et al., 17 Dec 2025)	Tempered time-fractional	SOE approach, tridiagonal	$O(h^2 + \tau^2)$ (graded)

6. Extensions and Limitations

• Operator-matching frameworks permit the systematic design of higher-order explicit and implicit schemes for more general PDEs (linear and nonlinear, constant or variable coefficients) (Chin, 2013).
• The structure-exploiting fast solvers remain effective in multi-dimensional and variable-order settings provided the underlying grid and operator are compatible with FFT or matrix Kronecker-product structure (Hao et al., 15 Jun 2024).
• Main limitations arise when the coefficient matrices lose Toeplitz/BTTB structure (e.g., for non-uniform grids or nonperiodic boundary conditions), but in such cases, local embedding and truncated preconditioners restore much of the acceleration (Jian et al., 2020).
• For subdiffusion with non-smooth initial data, special graded grids and operator quadrature/recursion are essential to preserving accuracy (Huang et al., 17 Dec 2025).

7. Summary and Perspective

Second-order fast finite difference schemes represent a robust and theoretically justified class of methods for time-dependent PDEs with nonlocal or fractional structure, offering second-order accuracy in time and space, unconditional stability, and computational cost scaling as $O(N \log N)$ or better. Their development relies on rigorous discrete analysis, advanced quadratures for noninteger derivatives, and fast iterative solvers for structured matrices. The combination of high-order accuracy and near-linear complexity positions these schemes as foundational for modern simulation of anomalous transport, viscoelasticity, fractional diffusion, and related fields. Future work may further extend these paradigms to adaptive, irregular meshes and fully nonlinear systems while preserving these computational and theoretical advantages.