Multi-Stage Multi-Derivative Framework

• The MSMD framework is a method that incorporates multiple stages and higher-order time derivatives to achieve high-order accuracy and extended stability in time discretization.
• It breaks classical Runge–Kutta barriers by enabling robust explicit schemes with enhanced CFL limits and effective shock-capturing through blending and limiting strategies.
• Its application in flux reconstruction and SSP methods demonstrates improved conservation, Fourier stability, and accurate resolution of nonlinear hyperbolic conservation laws.

The Multi-Stage Multi-Derivative (MSMD) framework unifies a broad class of time discretization methods for solving ordinary and partial differential equations, particularly in the context of nonlinear hyperbolic conservation laws. By employing multiple stages and derivatives in the time integration, the MSMD approach enables the construction of high-order numerical schemes with enhanced stability and admissibility properties, breaking classical barriers associated with standard Runge–Kutta and single-derivative methods. This framework underpins recent advances in explicit strong stability preserving (SSP) time-stepping and flux reconstruction (FR) approaches, yielding schemes with optimal accuracy, enlarged CFL stability bounds, and robust preservation of physically admissible states (Babbar et al., 2024, Christieb et al., 2015).

1. General Formulation of the MSMD Framework

The MSMD methodology generalizes explicit time-stepping methods by incorporating higher-order time derivatives in each stage. Consider the system

$u_t = F(u)$

where $u$ is a spatially discretized solution vector. The $s$-stage, $m$-derivative explicit MSMD method advances $u^n \rightarrow u^{n+1}$ as

$$\begin{aligned} Y_i &= u^n + \Delta t \sum_{j=1}^{i-1} a^{(1)}_{ij}\,F(Y_j) + \Delta t^2 \sum_{j=1}^{i-1} a^{(2)}_{ij}\,\dot F(Y_j) + \cdots + \Delta t^m \sum_{j=1}^{i-1} a^{(m)}_{ij}\,F^{(m)}(Y_j), \ u^{n+1} &= u^n + \Delta t \sum_{i=1}^s b^{(1)}_i\,F(Y_i) + \Delta t^2 \sum_{i=1}^s b^{(2)}_i\,\dot F(Y_i) + \cdots + \Delta t^m \sum_{i=1}^s b^{(m)}_i\,F^{(m)}(Y_i). \end{aligned}$$

Here, $\dot F(Y_j)$ and $F^{(m)}(Y_j)$ denote successive time derivatives of $F$ evaluated at different stage values. This structure permits higher-order accuracy with fewer stages compared to conventional Runge–Kutta methods and naturally generalizes to PDE contexts through the method-of-lines approach (Christieb et al., 2015).

2. Order Conditions and Strong Stability Preserving (SSP) Constraints

The accuracy of an MSMD scheme is controlled by a generalized B-series expansion, resulting in rooted-tree order conditions analogous to those governing classical Runge–Kutta methods. For example, up to fourth order, these include: $$\begin{aligned} & \text{Order 1: } b^{(1)}{}^T \mathbf{e} = 1, \ & \text{Order 2: } b^{(1)}{}^T c^{(1)} + b^{(2)}{}^T \mathbf{e} = \frac12, \ & \text{Order 3: } b^{(1)}{}^T (c^{(1)})^2 + 2 b^{(2)}{}^T c^{(1)} = \frac13, \quad b^{(1)}{}^T (A^{(1)} c^{(1)}) + b^{(1)}{}^T c^{(2)} + b^{(2)}{}^T c^{(1)} = \frac16, \ & \text{Order 4: } \ldots \end{aligned}$$ with higher-order relations available up to desired $p$ (Christieb et al., 2015).

For the preservation of strong stability properties (such as total variation diminishing, TVD), componentwise positivity conditions on the method's update matrices must be enforced. If the spatial semi-discretization satisfies forward Euler strong-stability and a second-derivative bound, sufficient conditions involve block positivity: $$R = (I + rS + \frac{r^2}{K^2}\hat{S})^{-1} \geq 0,\quad P = r R S \geq 0,\quad Q = \frac{r^2}{K^2} R \hat{S} \geq 0,$$ where $r$ is proportional to the time-step safety factor. The maximal $r$ defines the SSP coefficient $C$. This forms the basis of the optimization strategies for constructing high-order explicit SSP MSMD schemes (Christieb et al., 2015).

3. Application: MDRK/MSMD in Flux Reconstruction for Conservation Laws

The MSMD approach is concretely realized in the fourth-order, two-stage Multi-Derivative Runge–Kutta (MDRK) scheme for nonlinear conservation laws within the Flux Reconstruction (FR) framework. For the semi-discrete system

$\frac{du}{dt} = L(u) = -\partial_x f(u),$

the two-stage MDRK advances are: $$\begin{aligned} u^* &= u^n + \frac{1}{2}\Delta t\,L(u^n) + \frac{1}{8}\Delta t^2\,L_t(u^n), \ u^{n+1} &= u^n + \Delta t\,L(u^n) + \frac{\Delta t^2}{6}[L_t(u^n) + 2L_t(u^*)]. \end{aligned}$$ Within FR, this is written at solution points $\xi_p$ as

$$\begin{aligned} u_p^* &= u_p^n - \frac{\Delta t}{2\Delta x_e} \partial_\xi \overline{F}(\xi_p), \ u_p^{n+1} &= u_p^n - \frac{\Delta t}{\Delta x_e} \partial_\xi \overline{F}^*(\xi_p), \end{aligned}$$

where $\overline{F}$ and $\overline{F}^*$ are time-averaged flux reconstructions,

$$\overline F := f(u^n)+\frac14 \Delta t\,f_t(u^n),\qquad \overline F^* := f(u^n) + \frac16 \Delta t [f_t(u^n) + 2f_t(u^*)].$$

Time derivatives $f_t$ are approximated with a high-order “approximate Lax–Wendroff” finite difference scheme to yield $O(\Delta t^4)$ accuracy (Babbar et al., 2024).

4. Stability, Admissibility, and Limiting Strategies

Numerical stability is improved in the MSMD-FR approach by defining a D2 (time-averaged) Rusanov numerical flux,

$$\mathcal F_{e+1/2} = \tfrac12[\overline F^- + \overline F^+] - \tfrac12 \lambda_{e+1/2}[\overline u^+ - \overline u^-],$$

with $\lambda_{e+1/2} = \max\{|f'(u_e)|, |f'(u_{e+1})|\}$ and $\overline u^\pm$ appropriately defined stage-averaged states. Fourier analysis reveals a substantial increase in allowable CFL numbers—up to $\sigma_{\max} \approx 0.22$ for $g_2$ correction functions—yielding a 20–40% stability improvement over conventional dissipation mechanisms (Babbar et al., 2024).

For shock-capturing and admissibility preservation, a subcell-based blending limiter is applied. The solution update in each cell is given by

$$u_e = (1 - \theta_e) u_e^H + \theta_e u_e^L,\quad 0 \leq \theta_e \leq 1,$$

where $u_e^H$ and $u_e^L$ are high- and low-order updates respectively, and $\theta_e$ is determined from a smoothness indicator. Admissibility (e.g., positivity of density and pressure for Euler equations) in means is guaranteed provided the low-order scheme is admissible. Local flux limiting is used to enforce positivity adjacent to each cell face, as established in Theorem 3.1 of (Babbar et al., 2024). Dissipation is minimized by employing Gauss–Legendre solution points and MUSCL–Hancock reconstruction on subcells.

5. Concrete Schemes and Optimization

Key MSMD schemes include the fourth-order, two-stage MDRK for FR (Babbar et al., 2024), and the explicit three-stage, two-derivative SSP method achieving fifth-order accuracy, as constructed by Christlieb, Gottlieb, Grant, and Seal (Christieb et al., 2015). The latter has Butcher-like coefficients specifically optimized for maximal SSP coefficient. For the fifth-order, three-stage scheme (for $K=1$),

$$\begin{aligned} Y_1 &= u^n,\ Y_2 &= u^n + 0.7415 \Delta t\,F(Y_1) + 0.2743 \Delta t^2 \dot F(Y_1),\ Y_3 &= u^n + 0.2930 \Delta t\,F(Y_1) + 0.0061 \Delta t^2 \dot F(Y_1) + 0.0368 \Delta t^2 \dot F(Y_2),\ u^{n+1} &= u^n + 1 \cdot \Delta t\,F(Y_1) + 0.0895 \Delta t^2 \dot F(Y_1) + 0.1038 \Delta t^2 \dot F(Y_2) + 0.3067 \Delta t^2 \dot F(Y_3). \end{aligned}$$

The maximal admissible SSP coefficient is $C \approx 0.7851$ (Christieb et al., 2015).

6. Analysis and Verification

Theoretical guarantees include:

• Conservation: Demonstrated by exact quadrature of the update, ensuring

$$\sum w_p u_p^{n+1} = \sum w_p u_p^n - \Delta t [F_{e+1/2} - F_{e-1/2}]$$

for the FR setting.

• Admissibility: The blended scheme is admissible in means under appropriate low-order update conditions, and can be made pointwise admissible by subsequent scaling limiters.
• Fourier Stability: D2 dissipation and $g_2$ correction yield $\sigma_{\max} \approx 0.22$.
• Numerical Results: MSMD-FR achieves fourth-order convergence in 1D advection, Burgers’ equation, and Euler equations for smooth cases, and demonstrates superior shock resolution—particularly with MUSCL–Hancock blending—over classical TVB limiting. In 2D, the scheme successfully resolves complex flows such as double Mach reflection and Rayleigh–Taylor instability, maintaining high-order accuracy and capturing new small-scale structures (Babbar et al., 2024).

7. Significance and Implications

The MSMD approach enables explicit schemes of higher order and stability than attainable with classical Runge–Kutta methods of equal stage number (breaking the RK order barrier). The framework systematically produces SSP methods by explicit construction and optimization; these methods have demonstrably sharper allowable time-steps and broader applicability to nonlinear conservation laws with challenging admissibility constraints. The MSMD-FR schemes combine optimal formal order, provable conservation, enhanced Fourier stability, and robust shock-capturing via blending and local limiting—properties crucial for accurate and reliable simulation of compressible flows and other hyperbolic PDEs (Babbar et al., 2024, Christieb et al., 2015).