Implicit-Explicit Time Discretization

• Implicit-explicit (IMEX) schemes are a time discretization method that splits stiff and nonstiff components, applying implicit methods to the stiff parts and explicit methods to the nonstiff parts.
• They enhance computational efficiency by permitting larger, stable time steps for stiff dynamics while maintaining accuracy in the nonstiff regime.
• IMEX methods are rigorously analyzed for stability and error, with successful applications in fluid dynamics, porous media flow, and broader multiscale problems.

An implicit-explicit (IMEX) time discretization scheme is a temporal integration methodology used to solve systems of differential equations containing both stiff and nonstiff terms by treating each component differently: stiff contributions are handled implicitly (allowing large stable time steps), while nonstiff components are advanced explicitly (maintaining computational efficiency). IMEX schemes are prevalent across a broad range of applications, including incompressible and compressible fluid dynamics, kinetic models, porous media, optimal control problems, and multiscale systems where disparate temporal scales coexist.

1. Core Principles and Mathematical Formulation

An IMEX scheme splits the right-hand side of a semidiscrete ODE or PDE system,

$\frac{d}{dt} u(t) = F(u(t), t) + G(u(t), t),$

into the nonstiff portion $F$ (treated explicitly) and the stiff portion $G$ (treated implicitly). For time discretization, a common approach is to use combinations such as explicit Runge–Kutta, Adams–Bashforth, or extrapolation methods for $F$, and diagonally implicit Runge–Kutta (DIRK), backward differentiation formulas (BDF), or Crank–Nicolson type schemes for $G$.

A prototypical two-level IMEX BDF scheme for time step $n \to n+1$ is: $$\frac{1}{\tau} (u^{n+1} - u^n) = F(\hat{u}^{n+1}) + G(u^{n+1}),$$ where $\hat{u}^{n+1}$ is an explicit extrapolation, e.g., $\hat{u}^{n+1} = 2u^n - u^{n-1}$. More generally, IMEX Runge–Kutta methods define each stage through

$$U^{(i)} = u^n + \Delta t \sum_j a_{ij} G(U^{(j)}) + \Delta t \sum_{j<i} \tilde{a}_{ij} F(U^{(j)}),$$

with the implicit and explicit tableau coefficients $(a_{ij}), (\tilde{a}_{ij})$ tailored for stability, accuracy, or SSP (strong-stability-preserving) properties. The final update is a convex combination of stage values.

This splitting philosophy is leveraged in a wide array of couplings—spatially or functionally—across stiff and nonstiff model contributions.

2. Scheme Construction, Stability, and Order Conditions

To achieve high-order accuracy and robust performance in the presence of stiffness, IMEX schemes are meticulously constructed to balance local truncation error, absolute monotonicity, and region of stability. For many systems, e.g., stiff kinetic or reaction terms, an all-stages-implicit (ASI) construction is necessary, particularly when singular limits ($\varepsilon \to 0$) arise. Here, each stage is solved implicitly and the scheme becomes unified across the singularly perturbed regime (Duan, 2016).

The stability constraints depend on the explicit portion: for nonstiff parts, CFL-type restrictions may apply, but implicit handling of stiff terms relaxes this constraint. Unconditional stability is attainable if the implicit integrator, such as BDF2 or backward Euler, is $A$-stable and the explicit–implicit coupling satisfies order conditions. For example, in the stabilized IMEX Crank–Nicolson method for the Oseen equation,

$$\tau^{-1}(u_h^{n+1} - u_h^n) + C_h^{n+\frac{1}{2}}\hat{u}_h^{n+1} + A_h\bar{u}_h^{n+1} + G_hp_h^{n+1} = L_h^{n+1},$$

the pressure–velocity coupling and viscous terms are discretized implicitly to ensure robustness at high Reynolds numbers, while convection is extrapolated explicitly (Burman et al., 21 May 2024).

Stiff accuracy—where the final stage weights are drawn from the last row of the implicit tableau—ensures that stage-order reduction does not destroy global accuracy, even in the presence of relaxation or singular perturbations (Albi et al., 2013).

3. Spatial Discretization and Stabilization

IMEX schemes are typically coupled with advanced spatial discretizations—finite element, finite volume, or discontinuous Galerkin methods—which interact intimately with the temporal splitting. Notably:

• Symmetric gradient-jump penalization (or interior penalty) is used to enforce $\inf\text{-}\sup$ stability of equal-order velocity–pressure pairs in incompressible flows, or to control oscillations in convection-dominated regimes (Burman et al., 21 May 2024).
• For multiscale or fractured media, coarsening/projection by nonlocal multicontinuum (NLMC) methods or meshfree generalized finite element methods reduces the system to a lower-dimensional subspace amenable to efficient IMEX time stepping (Nikiforov, 5 Feb 2024, Vasilyeva, 25 Apr 2024).

Spatial stabilization terms participate in error estimates and may induce suboptimal rates (e.g., $O(h^{k+1/2})$ instead of $O(h^{k+1})$) due to consistency–stability trade-offs.

4. Error Analysis, Superconvergence, and Unconditional Stability

Rigorous a priori error analysis is standard. For the IMEX Crank–Nicolson method applied to Oseen’s equation with symmetric jump stabilization, the error in the velocity satisfies

$$\|u(t^n) - u_h^n\| \le C( h^{\,k+\frac{1}{2}} + \tau^2 ),$$

where $h$ is the mesh size and $\tau$ the time step (Burman et al., 21 May 2024). For second-order BDF–IMEX schemes, spatial and temporal supercloseness and superconvergence properties can be established by decomposing the error into temporal and spatial components and employing special interpolants or postprocessing (Altmann et al., 23 Aug 2025). Fully discrete error bounds hold without step size restrictions, i.e., unconditional stability, by exploiting energy dissipation arguments or contractive mapping properties.

5. Applications: Fluid Dynamics, Multiscale Media, and Wave Propagation

IMEX discretizations see systematic use in computational fluid dynamics, particularly for high Reynolds number flows, evolutionary Hamilton–Jacobi–BeLLMan equations, multiscale diffusion, and wave propagation in media with high contrast or geometric stiffness:

• In optimal control governed by the Goldstein–Taylor or Hamilton–Jacobi equations, IMEX–Runge–Kutta or Eulerian IMEX schemes attain asymptotic preserving properties, permitting time steps independent of the parabolic or advective CFL constraint (Vladimirsky et al., 2013, Albi et al., 2013).
• For nonlinear thermistor or semilinear wave equations with strong damping, BDF–Galerkin or midpoint–explicit midpoint pairings maintain efficiency and high-order accuracy, with splitting of nonlinear source terms via extrapolation (Eckhardt et al., 28 Jun 2024, Altmann et al., 23 Aug 2025).
• In multicontinuum or fractured porous media, IMEX decoupling schemes partition the operator so that dominant block-diagonal (self-coupling) parts are handled implicitly and the coupling (off-diagonal) parts explicitly, leading to efficient, physically meaningful upscaled models (Vasilyeva, 25 Apr 2024).
• For wave equations and Maxwell’s equations, IMEX and implicit leapfrog methods enable unconditionally stable and high-order time discretizations compatible with domain decomposition, overset grids, or phase-field models (Carson et al., 22 Apr 2024, Arya et al., 7 Apr 2024).

6. Fractional Step and Partitioned Algorithms

A fundamental structural feature of many IMEX schemes is their amenability to operator splitting and fractional step (projection) methods:

• For the incompressible Oseen or Navier–Stokes equations, the IMEX Crank–Nicolson discretization naturally leads to a splitting where an explicit velocity update is followed by a pressure Poisson problem and then a velocity correction. This allows reusing efficient solvers for subcomponents, e.g., Stokes or Poisson solvers (Burman et al., 21 May 2024).
• For multiscale flow and diffusion, schemes may operate on split subspaces—rapid, high-contrast components are advanced implicitly, while smoother or slowly evolving modes are updated explicitly, enabling large time steps decoupled from the most severe local stiffness (Wang et al., 21 Mar 2024, Nikiforov, 5 Feb 2024).

Such approaches are essential for balancing computational efficiency and stability in settings where system coupling or localized stiffness would render monolithic treatment (fully implicit or explicit) prohibitively expensive.

7. Numerical Validation and Computational Implications

Extensive numerical testing accompanies IMEX scheme development:

• Benchmarking on canonical problems (Taylor–Green vortex, Kelvin–Helmholtz instability, flow around a cylinder) verifies optimal or near-optimal convergence rates and quantifies dissipation, robustness, or conservation (Burman et al., 21 May 2024).
• Adaptivity in time and/or space is frequently supported, with runtime gains observed—IMEX methods often yield speedup factors of 4–15 over explicit solvers in industrial-scale turbulent flow simulations when combined with tailored spatial hybridization techniques (Pereira et al., 30 Jan 2024).
• For problems with parameter-dependent stiffness, e.g., problems tending toward singular limits or high-contrast coefficients, IMEX schemes maintain uniform accuracy and stability across the parameter regime (Albi et al., 2013, Duan, 2016).
• Superconvergence and postprocessing techniques enable extraction of higher-order accuracy from IMEX–Galerkin schemes, especially in finite element settings with suitable mesh and projection operators (Altmann et al., 23 Aug 2025).

A notable limitation arises in the loss of accuracy for certain vanilla IMEX variants when coefficients (e.g., damping) are time-dependent; revised schemes using midpoint quadrature for such terms can restore the expected order of convergence (Jiao et al., 28 Oct 2024).

---

In summary, implicit–explicit time discretization schemes provide a unifying framework for the robust and efficient integration of evolutionary PDEs exhibiting both stiff and nonstiff mechanisms. Their design incorporates a nuanced balance of stability, accuracy, and computational tractability, enabling wide applicability to high-dimensional, multiscale, and multiphysics systems while supplying rigorous a priori guarantees on convergence and qualitative behavior. Their ongoing development continues to address emerging challenges in multiscale modeling, domain decomposition, and parallel-in-time integration for stiff and geometrically complex physical systems.