Crank-Nicolson Approach in Fractional PDEs

• The Crank-Nicolson approach is a family of implicit, time-centered numerical schemes that achieve stability and accuracy in solving differential and integro-differential equations.
• It employs a trapezoidal time discretization combined with fourth-order spatial finite differences to efficiently handle fractional derivatives and stiff dynamics.
• The method ensures unconditional stability and optimal error bounds through discrete fractional-energy arguments and predictor-corrector schemes.

The Crank-Nicolson Approach is a family of implicit, time-centered numerical schemes renowned for their stability and accuracy in solving diverse classes of differential and integro-differential equations, especially those involving fractional derivatives, dispersive and advective terms, and stiff or degenerate parabolic dynamics. Originating as a second-order method for parabolic PDEs, the Crank-Nicolson methodology now encapsulates a wide spectrum of extensions—including multi-step, high-order, and fractional-derivative variants—engineered for demanding applications in transport, anomalous diffusion, and memory-dominated physical models. Its appeal arises from unconditional stability, controlled accuracy, and flexibility in spatial and temporal discretization, particularly in regimes where explicit schemes are hampered by stringent CFL-type restrictions or inefficient low-order convergence.

1. Theoretical Foundation and Discrete Formulation

At its core, the Crank-Nicolson approach centers on a trapezoidal, time-centered discretization of differential equations of the form

$u_t = \mathcal{L}u + f,$

producing the canonical time-stepping stencil

$$\frac{U^{n+1} - U^n}{k} = \frac12 \left( \mathcal{L}U^{n+1} + \mathcal{L}U^n \right) + \frac12\left(f^{n+1} + f^n\right),$$

where $U^n$ approximates $u(t_n)$ and $k$ is the time step. In the classical setting, this yields second-order accuracy in time and, with proper spatial discretization, high-order accuracy in space.

For equations involving memory or non-local fractional time derivatives, such as the time-variable fractional mobile-immobile advection-dispersion equation

$$u_t(x,t) = -c\,D_t^\beta u + u_{xx} - u_x + f(x,t), \quad 0 < \beta(x,t) < 1,$$

the Crank-Nicolson framework is nontrivial. Here, repeated convolution-type sums and shifted quadrature rules are invoked to approximate nonlocal terms, with precise weights and remainder bounds required for stability and consistency. In the approach of (Ngondiep, 2022), the temporal grid is enriched by half-steps and predictor-corrector sublevels to align the discrete scheme with the internal structure of the Caputo derivative and to facilitate optimal error propagation.

2. Fourth-Order Two-Step Explicit Euler/Crank-Nicolson Hybrid Scheme

The scheme analyzed in (Ngondiep, 2022) achieves higher spatial accuracy via fourth-order finite-difference stencils: $$u_x(x_j) \approx \delta_x^4 u_j = \frac{-u_{j+2} + 8u_{j+1} - 8u_{j-1} + u_{j-2}}{12h},$$

$$u_{xx}(x_j) \approx \delta_{2x}^4 u_j = \frac{-u_{j+2} + 16u_{j+1} - 30u_j + 16u_{j-1} - u_{j-2}}{12h^2},$$

paired with a two-step, time-shifted Euler/Crank-Nicolson hybridization:

• Predictor: an explicit-Euler step employs L1-type convolution quadrature for the Caputo derivative and spatial fourth-order discretization, solving a pentadiagonal linear system at the fractional time level $t_{n+\frac12+\alpha}$.
• Corrector: a Crank-Nicolson-type step at $t_{n+1+\alpha}$ further balances temporal errors and refines the solution by including additional weighted fractional-derivative terms and quadrature remainders.

Notationally, primary unknowns $U_j^n \approx u(x_j,t_n)$ propagate on grids $x_j = L_0 + jh$, $t_n = nk$, and the Caputo derivative is discretized at fractional time levels using precomputed, decaying, strictly positive weights $a_{n+\frac12,\ell+\frac12}^{\alpha,\beta}$.

3. Stability and Error Analysis

The main theoretical result is unconditional stability and optimal error bounds in the strong $L^\infty(0,T;L^2)$ norm, even for non-smooth fractional solutions. The crucial analytic components are:

• Quadrature Control: The remainder in the fractional-derivative quadrature is provably $O(k^{2-\beta})$ (see Lemma 1 of (Ngondiep, 2022)), with the weights $a_{n+\frac12,\ell+\frac12}^{\alpha,\beta}$ shown to be positive and strongly decaying (Lemma 2).
• Discrete Fractional-Energy Argument: By employing telescoping in time and "fractional-energy" inequalities,

$$(u^{n+s}, \Delta_t^{\beta} u^{n+s+\alpha}) \geq \frac12 \Delta_t^\beta ( (u^{n+s+\alpha})^2 ),$$

all history terms due to the Caputo derivative are controlled.

• Operator Dissipativity: The fourth-order spatial operator

$L_h = \delta_{2x}^4 - \delta_x^4$

is dissipative, with

$$(L_h u, u) \leq 0, \qquad (-L_h u, u)\geq c_1 \|\delta_x u\|_2^2,$$

enforcing energy decay at the discrete level.

• Convergence: For $0 < \beta < 2/3$, $0 < \alpha < 1/2$, one obtains

$$\max_{0\leq n\leq N-1} \|U^{n+\frac12}\|_2, \ \max_{0\leq n\leq N} \|U^n\|_2 \leq \|u\|_{L^\infty(0,T;L^2)} + C(k + k^{2-\beta} + k^2 + h^4)$$

and—using discrete Grönwall arguments—

$$\max_n \|e^n\|_2 + \max_n \|e^{n+\frac12}\|_2 \leq C\,(k + h^4)$$

where $e^n = U^n - u(x_j, t_n)$. Thus, the method achieves first-order accuracy in time and fourth-order in space, with unconditional stability ensured by the dissipativity and energy telescoping.

4. Computational Implementation and Practical Considerations

Algorithmic features include:

• Linear Algebra: All implicit solves are for pentadiagonal, nonsymmetric systems arising from the high-order spatial discretization. These are efficiently handled by direct band solvers or Krylov subspace methods (e.g., GMRES) with band-circulant preconditioners.
• Fractional Weights: The a priori computation of all required weights $a_{n+\frac12,\ell+\frac12}^{\alpha,\beta}$ is recommended; these can be reused across time steps due to the time-invariance for constant $\beta$.
• No CFL Constraint: The unconditional stability of the method obviates any explicit time-step restriction; for optimal accuracy, it is natural to choose $k \sim h^4$.
• Boundary Conditions: Dirichlet conditions are enforced directly; ghost values may be set to prescribed boundary data without affecting the solvability or stability of the pentadiagonal systems.

The step-by-step procedure involves initialization, computation of a shifted start value $U^\alpha$, predictor-corrector sweeps at each time level, and direct enforcement of boundary values for all temporal sub-levels.

5. Comparison with Classical Crank-Nicolson and Other Schemes

Relative to standard integer-order Crank-Nicolson (CN) methods:

• Classical CN: Second-order in time and usually second-order in space unless augmented with compact higher-order stencils. Unconditional stability is standard but proven only for linear, constant-coefficient problems.
• Modified Two-Step CN: The present approach features a hybrid explicit Euler (L1-type) predictor and a CN-type corrector, shifted by an adjustable parameter $\alpha$. This adaptation is necessary for fractional time derivatives, which induce nonlocal temporal memory, and achieves first-order temporal accuracy (optimal for fractional dynamics with limited time regularity) and fourth-order spatial accuracy, with unconditional stability and pentadiagonal algebraic structure. It avoids step-size constraints that often accompany explicit/singly-implicit handling of the fractional term.

A comparative summary:

Feature	Classical CN (PDE)	Modified Euler/CN (Fractional)
Time Accuracy	Second-order	First-order (optimal for fractional)
Space Accuracy	Second-order	Fourth-order
Stability	Unconditional (linear)	Unconditional (fractional, any k)
Solve Structure	Tridiagonal	Pentadiagonal
Fractional Support	Not directly	Yes (Caputo, variable β)

6. Extensions and Generalizations

The explicit–implicit predictor-corrector and the energy-based unconditional-stability framework are extensible to broader classes of nonlocal, subdiffusive, or multi-term fractional PDEs. The technical machinery—discrete fractional-convolution quadrature, telescopic energy inequalities, high-order spatial discretization, and two-step start-up—carries over to systems with spatially and temporally variable fractional order, nonlinear memory, and stiff multi-physics couplings. For further details regarding the development and application of these high-order fractional schemes, and precise stability arguments, refer to the full theoretical and numerical discussion in (Ngondiep, 2022).