Stabilized IMEX Crank–Nicolson Method

• The paper introduces a stabilized IMEX Crank–Nicolson scheme that splits stiff linear terms and non-stiff nonlinear components to achieve second-order temporal accuracy.
• It employs targeted stabilization techniques such as continuous interior penalty, biharmonic modifications, and doubly-stabilized splittings to ensure robust energy dissipation under severe conditions.
• Numerical experiments confirm its effectiveness in high-Reynolds, convection–diffusion, and phase-field models, offering uniform error estimates and computational efficiency.

The stabilized implicit–explicit (IMEX) Crank–Nicolson method is a class of time-stepping schemes for parabolic, convection–diffusion, incompressible flow, and phase-field problems, designed to achieve second-order temporal accuracy, unconditional or robust stability, and efficient linear algebra at each step. The essential idea is to split the governing equation into linear (stiff) and nonlinear or advective (non-stiff) parts, treating the former implicitly via Crank–Nicolson averaging, and the latter explicitly, often with extrapolation and stabilization. To maintain stability and structure-preservation under severe time step or spatial mesh constraints, targeted stabilization mechanisms such as continuous interior penalty (CIP) terms, biharmonic modification, or doubly-stabilized splittings are incorporated. These methods have seen extensive analysis and application, notably for high-Péclet or high-Reynolds convection–diffusion, Oseen/Navier–Stokes flows, Cahn–Hilliard, and Allen–Cahn equations (Burman et al., 2020, Wang et al., 2017, Orizaga et al., 2024, Wang et al., 2020, Burman et al., 2024, Hou et al., 2023).

1. Abstract Formulation and Model Problems

The stabilized IMEX Crank–Nicolson method applies to evolution equations of the abstract form

$$\partial_t u = F(u) = \mathcal{L}(u) + \mathcal{N}(u),$$

where $\mathcal{L}$ is a linear (possibly stiff, such as diffusive or biharmonic) component and $\mathcal{N}$ is a nonlinear, advective, or otherwise explicitly treated part.

Typical model equations include:

• Convection–Diffusion: $$\partial_t u + \beta\cdot\nabla u - \mu \Delta u = f$$ (Burman et al., 2020)
• Oseen/Navier–Stokes: $$\partial_t u + \beta\cdot\nabla u + \nabla p - \nu \Delta u = f,\quad\nabla\cdot u=0$$ (Burman et al., 2024)
• Cahn–Hilliard and Thin-Film: $$\partial_t u = \nabla\cdot[M(u)\nabla(W'(u) - \Delta u)]$$ (Orizaga et al., 2024)
• Allen–Cahn with Mobility: $$\partial_t \phi = -M(\phi) (-\varepsilon^2 \Delta \phi + F'(\phi))$$ (Hou et al., 2023)

In each case, stabilization is introduced to control explicit terms and improve robustness at large time steps or high Reynolds/Péclet numbers.

2. Temporal Discretization Scheme

Stabilized IMEX Crank–Nicolson schemes combine implicit Crank–Nicolson averaging for stiff linear terms with explicit, often extrapolated or Adams–Bashforth-type, evaluation for nonlinear/convective components. Stabilizing corrections of $O(\tau^2)$ (where $\tau$ is the time step) are added to preserve energy dissipation, mass conservation, or maximum principles as required.

Generic two-level IMEX Crank–Nicolson scheme: $$(u^{n+1}-u^n, v) + \tau\,\mathcal{C}_{\text{stab}}(\hat{u}^{n+1}, v) + \tau\,\mathcal{A}\bigg(\frac{u^{n+1}+u^n}{2}, v\bigg) = \tau(f^{n+1}, v)$$ where:

• $\mathcal{L}$0 is the symmetric, implicit (e.g., diffusive) part,
• $\mathcal{L}$1 is the explicitly stabilized advective/nonlinear term, typically evaluated at an extrapolation $\mathcal{L}$2 (for example, $\mathcal{L}$3).

For problems with additional nonlinear bulk forces or degeneracies (as in Cahn–Hilliard or Allen–Cahn systems), further linear stabilization terms are introduced, often involving higher-order time differences and Laplacians (Wang et al., 2017, Hou et al., 2023, Wang et al., 2020).

3. Stabilization Techniques and Spatial Discretization

Continuous Interior Penalty (CIP) and Jump Stabilization

For finite element methods, especially in convection–dominated regimes, symmetric stabilization terms such as continuous interior penalty are required for robustness. For example, the stabilized bilinear form for the convection term is: $\mathcal{L}$4 Here, $\mathcal{L}$5 is the stabilization parameter, and $\mathcal{L}$6 is the jump of the gradient across interior faces (Burman et al., 2020).

Biharmonic and Doubly-Stabilized Splittings

For fourth-order problems (e.g., Cahn–Hilliard, thin-film), the splitting $\mathcal{L}$7 is adopted and further stabilization is enforced via biharmonic terms (with $\mathcal{L}$8 chosen relative to the maximum mobility) (Orizaga et al., 2024).

In phase-field models, doubly-stabilized IMEX CN schemes insert both first-difference Laplacian and second-difference (leapfrog-type) terms, denoted as $\mathcal{L}$9, with $\mathcal{N}$0, $\mathcal{N}$1 chosen for unconditional stability and maximum bound principle (MBP) preservation (Wang et al., 2017, Hou et al., 2023).

4. CFL Conditions, Stability, and Error Analysis

Stability of the scheme is governed by coupled spatial (mesh) and temporal (time step) constraints, generally expressed through CFL-type numbers related to mesh Reynolds/Péclet: $\mathcal{N}$2 For $\mathcal{N}$3th-degree finite elements:

• Piecewise affine $\mathcal{N}$4: standard hyperbolic CFL suffices.
• Higher degree $\mathcal{N}$5: a stricter $\mathcal{N}$6-CFL, $\mathcal{N}$7, is required for stability (Burman et al., 2020, Burman et al., 2024).

Main stability and error results:

• Discrete energy dissipation: For appropriate stabilization parameters, modified discrete energies decrease monotonically for any $\mathcal{N}$8, i.e., unconditional stability (Wang et al., 2017, Wang et al., 2020, Hou et al., 2023, Orizaga et al., 2024).
• Uniform bounds: Solutions are bounded in natural norms (e.g., $\mathcal{N}$9, $$\partial_t u + \beta\cdot\nabla u - \mu \Delta u = f$$0, $$\partial_t u + \beta\cdot\nabla u - \mu \Delta u = f$$1 for MBP).
• Error order: $$\partial_t u + \beta\cdot\nabla u - \mu \Delta u = f$$2- and $$\partial_t u + \beta\cdot\nabla u - \mu \Delta u = f$$3-error estimates of $$\partial_t u + \beta\cdot\nabla u - \mu \Delta u = f$$4 for convection–diffusion/Oseen; $$\partial_t u + \beta\cdot\nabla u - \mu \Delta u = f$$5 (with $$\partial_t u + \beta\cdot\nabla u - \mu \Delta u = f$$6 depending on spatial discretization) for phase-field models (Burman et al., 2020, Burman et al., 2024, Wang et al., 2017, Hou et al., 2023).

5. Algorithmic Details and Practical Implementation

At each time step, the stabilized IMEX Crank–Nicolson method necessitates solving systems with constant-coefficient, self-adjoint (possibly block) elliptic operators, as nonlinearity is handled explicitly. This enables the use of efficient direct or preconditioned iterative solvers, becoming particularly advantageous for spectral or uniform grid discretizations.

• In mixed formulations (Cahn–Hilliard), the coupled $$\partial_t u + \beta\cdot\nabla u - \mu \Delta u = f$$7 system reduces to a block $$\partial_t u + \beta\cdot\nabla u - \mu \Delta u = f$$8 system with time-invariant structure, admitting efficient block preconditioners (Wang et al., 2020).
• Adaptive time-stepping, based either on interface velocity or energy variation, further enhances efficiency during slow coarsening or equilibrium (Wang et al., 2020, Hou et al., 2023).

Stabilization parameters may be taken as minimal as possible in practice once empirical stability is observed, even though theoretical bounds may prescribe large values (e.g., $$\partial_t u + \beta\cdot\nabla u - \mu \Delta u = f$$9, $$\partial_t u + \beta\cdot\nabla u + \nabla p - \nu \Delta u = f,\quad\nabla\cdot u=0$$0 for phase-field models, with $$\partial_t u + \beta\cdot\nabla u + \nabla p - \nu \Delta u = f,\quad\nabla\cdot u=0$$1 a global Lipschitz constant) (Wang et al., 2017, Wang et al., 2020).

6. Applications, Numerical Experiments, and Extensions

The stabilized IMEX Crank–Nicolson method supports a broad spectrum of computational PDEs:

• High Reynolds/Péclet flows: Demonstrated with Oseen and Navier–Stokes equations for vortex dynamics, mixing layers, and wake flows at high Reynolds number, yielding robust error orders (e.g., $$\partial_t u + \beta\cdot\nabla u + \nabla p - \nu \Delta u = f,\quad\nabla\cdot u=0$$2 for $$\partial_t u + \beta\cdot\nabla u + \nabla p - \nu \Delta u = f,\quad\nabla\cdot u=0$$3 elements) and stable energy dissipation (Burman et al., 2024).
• Phase-field dynamics: For Cahn–Hilliard and Allen–Cahn systems with degenerate mobility or sharp interface limits, the method ensures unconditional MBP, energy decay, and spectral/higher-order spatial accuracy, supporting adaptive time stepping for long-term evolution (Wang et al., 2017, Wang et al., 2020, Orizaga et al., 2024, Hou et al., 2023).
• Thin-film equations: Adaptation of the biharmonic stabilization yields energy-stable and mass-conserving solvers for dewetting, coarsening, and other fourth-order gradient flows (Orizaga et al., 2024).

Extensions include fractional step (“projection”) methods in incompressible flow, where the IMEX CN scheme underlies operator-splittings for pressure-velocity coupling without loss in order or structure preservation (Burman et al., 2024).

7. Theoretical Results and Parameter Selection

Central theorems and guiding results across model problems include:

• Energy stability: For stabilization parameters $$\partial_t u + \beta\cdot\nabla u + \nabla p - \nu \Delta u = f,\quad\nabla\cdot u=0$$4 or $$\partial_t u + \beta\cdot\nabla u + \nabla p - \nu \Delta u = f,\quad\nabla\cdot u=0$$5 satisfying explicit algebraic bounds (e.g., $$\partial_t u + \beta\cdot\nabla u + \nabla p - \nu \Delta u = f,\quad\nabla\cdot u=0$$6 with $$\partial_t u + \beta\cdot\nabla u + \nabla p - \nu \Delta u = f,\quad\nabla\cdot u=0$$7 for IMEX2), unconditional discrete energy dissipation is guaranteed for any time step (Wang et al., 2017, Orizaga et al., 2024, Hou et al., 2023).
• Error estimates: For smooth solutions and appropriate time/space steps, global error bounds are established as $$\partial_t u + \beta\cdot\nabla u + \nabla p - \nu \Delta u = f,\quad\nabla\cdot u=0$$8 (convection–diffusion), $$\partial_t u + \beta\cdot\nabla u + \nabla p - \nu \Delta u = f,\quad\nabla\cdot u=0$$9 (phase-field), with prefactors depending algebraically (not exponentially) on small parameters such as $$\partial_t u = \nabla\cdot[M(u)\nabla(W'(u) - \Delta u)]$$0 (Burman et al., 2020, Wang et al., 2017, Wang et al., 2020).
• Maximum principle preservation: For Allen–Cahn dynamics with variable mobility, the scheme preserves $$\partial_t u = \nabla\cdot[M(u)\nabla(W'(u) - \Delta u)]$$1 unconditionally (for sufficiently large stabilizers), with associated $$\partial_t u = \nabla\cdot[M(u)\nabla(W'(u) - \Delta u)]$$2-error estimates (Hou et al., 2023).

Tables summarizing recommended stabilization parameters and their empirical effects are reported in (Orizaga et al., 2024, Hou et al., 2023).

---

References:

• (Burman et al., 2020) E. Burman, J. Guzmán. "Implicit-explicit multistep formulations for finite element discretisations using continuous interior penalty"
• (Burman et al., 2024) E. Burman, D. Garg, J. Guzmán. "Implicit-explicit Crank-Nicolson scheme for Oseen's equation at high Reynolds number"
• (Wang et al., 2017) X. Wang, J. Yu. "Convergence Analysis of an Unconditionally Energy Stable Linear Crank-Nicolson Scheme for the Cahn-Hilliard Equation"
• (Wang et al., 2020) P. Cheng, F. Yang, H. Wu. "An Energy Stable Linear Diffusive Crank-Nicolson Scheme for the Cahn-Hilliard Gradient Flow"
• (Orizaga et al., 2024) M. Orizaga, T. Witelski. "IMEX methods for thin-film equations and Cahn-Hilliard equations with variable mobility"
• (Hou et al., 2023) X. Xu, P. Gao, Z. Qiao, T. Tang. "A linear doubly stabilized Crank-Nicolson scheme for the Allen-Cahn equation with a general mobility"