---
title: Runge-Kutta Discontinuous Galerkin Schemes
url: https://www.emergentmind.com/topics/runge-kutta-discontinuous-galerkin-schemes
type: topic
---

# Runge-Kutta Discontinuous Galerkin Schemes

A Runge-Kutta Discontinuous Galerkin (RKDG) scheme is a high-order finite element method for numerically solving time-dependent partial differential equations, notably hyperbolic conservation laws. The RKDG framework applies a spatial discretization by discontinuous Galerkin (DG) finite elements, combined with strong-stability-preserving (SSP) explicit Runge-Kutta time integration. When the flux function contains spatially discontinuous coefficients—i.e., $u_t + f\big(u, \theta(x)\big)_x = 0$ where $\theta(x)$ is a piecewise constant parameter—the DG scheme requires special treatment at interfaces of discontinuity to retain steady-state preservation, high accuracy, and non-oscillatory shock resolution. The hybrid RKDG scheme augmented with $\delta$-mapping achieves these goals by unifying left/right states at discontinuous interfaces and adjusting limiters to maintain entropy and supply-demand consistency [1511.01346].

## 1. Problem Formulation and Flux Discontinuity

Consider the initial-value problem for a hyperbolic conservation law with spatially variable flux parameters,
\[
u_t + \partial_x f(u, \theta(x)) = 0, \qquad x \in [0, L],~t > 0,
\]
where $u \in \mathbb{R}^r$ and $\theta(x) \in \mathbb{R}^q$ is piecewise constant, possibly discontinuous at mesh interfaces $x_{j+1/2}$. At such interfaces, denote the left and right limits by $\theta^-(x_{j+1/2}) = \theta_j$ and $\theta^+(x_{j+1/2}) = \theta_{j+1}$, respectively.

Classical steady or stationary-shock solutions across jumps satisfy the local Rankine–Hugoniot condition,
\[
[f(u, \theta_j)] = 0 \quad\text{across each jump of }\theta,
\]
ensuring conservation is retained in the presence of spatially discontinuous flux.

## 2. Spatial Discretization via Discontinuous Galerkin

Partition $[0, L]$ into cells $I_j = (x_{j-1/2}, x_{j+1/2})$. On each cell, approximate $u$ by polynomials of degree $\leq k$,
\[
u_h(x, t)\big|_{I_j} = \sum_{l=0}^k u_j^l(t)\, \varphi_j^l(x),
\]
where $\{\varphi_j^l(x)\}$ is a local basis (e.g., scaled Legendre polynomials).

The weak DG form, for each basis index $l$, reads
\[
\int_{I_j} (u_h)_t\, \varphi_j^l\,dx
- \int_{I_j} f(u_h, \theta)\, \partial_x \varphi_j^l\,dx
+ \widehat{f}_{j+1/2}\, \varphi_j^l(x_{j+1/2}^-)
- \widehat{f}_{j-1/2}\, \varphi_j^l(x_{j-1/2}^+)
= 0.
\]

After expansion, the method yields a finite-dimensional ODE system for each coefficient,
\[
\frac{d}{dt} u_j^l(t) = \frac{2l + 1}{\Delta_j} \left[
\int_{I_j} f(u_h, \theta)\, \partial_x \varphi_j^l\,dx
- \left\{ (-1)^l \widehat{f}_{j-1/2} + \widehat{f}_{j+1/2} \right\}
\right].
\]

## 3. Numerical Flux Definition via δ-Mapping

Direct use of standard numerical fluxes at interfaces where $\theta$ jumps fails to preserve steady states. Instead, the $\delta$-mapping algorithm maps left/right DG traces $(u^-_{j+1/2}, u^+_{j+1/2})$ to unified states $(\delta u^-_{j+1/2}, \delta u^+_{j+1/2})$ evaluated at an intermediate parameter $\bar\theta$, defined to satisfy:

- **Supply-Demand Equalization**: maximize $\gamma \leq 1$ subject to
  \[
  f(\delta u^\mp, \theta_{j+1/2}) = \gamma f(u^\mp, \theta_{j+1/2\mp 1/2}),
  \]
- **Entropy and Eigenvalue Compatibility**:
  \[
  \lambda_l(\delta u^\mp, \theta_{j+1/2})\, \lambda_l(u^\mp, \theta_{j+1/2\mp 1/2}) \geq 0,
  \]
  ensuring no sign change for characteristic speeds.

The mapped states are used in a classical Riemann solver,
\[
\widehat{f}_{j+1/2} = \widetilde{f}(\delta_{j+1/2} u^-_{j+1/2}, \delta_{j+1/2} u^+_{j+1/2}, \theta_{j+1/2}).
\]

## 4. Slope Limiting with δ-Mapping

To control nonphysical oscillations near shocks and discontinuities, a TVB-type minmod slope limiter is applied before each RK stage. For $k=1$,
\[
u_h|_{I_j}(x) = \bar{u}_j + \frac{2u_j^1}{\Delta_j}(x - x_j),
\]
with slope replaced by a δ-mapped minmod,
\[
\Lambda\Pi_h^1 u_h = \bar{u}_j + \frac{2}{\Delta_j} m(u_j^1, \delta_j \bar{u}_{j+1} - \bar{u}_j, \bar{u}_j - \delta_j \bar{u}_{j-1})(x - x_j),
\]
where $\delta_j$ maps neighboring cell averages to $\theta_j$, guaranteeing steady-state preservation.

## 5. SSP Runge–Kutta Time Integration

Let the semi-discrete DG operator be $L_h(u_h)$. Time integration is performed by an SSP RK scheme of order $k+1$,
\[
u_h^{(0)} = u_h^n, \quad
u_h^{(i)} = \sum_{m=0}^{i-1} \alpha_{im}\, u_h^{(m)}
+ \Delta t \sum_{m=0}^{i-1} \beta_{im}\, L_h(u_h^{(m)}),
\]
setting $u_h^{n+1} = u_h^{(k)}$ at the completion of the step.

The coefficients $(\alpha_{im}, \beta_{im})$ are taken from the TVD RK (Cockburn–Shu) Butcher tableau.

CFL stability constraint is
\[
\Delta t \leq C\, \min_j \frac{\Delta_j}{\max_i |\lambda_i|}, \quad C = \frac{1}{2k+1},
\]
with characteristic speeds evaluated on δ-mapped states. The formal accuracy is $\mathcal{O}(h^{k+1})$.

## 6. Algorithmic Workflow

The full algorithm proceeds as:

1. Decompose the domain into cells $I_j$ and project $u$ onto $P^k(I_j)$.
2. At each discontinuous interface, compute δ-mapped states by supply-demand and entropy constraints.
3. Construct the DG spatial operator $L_h$ using the unified flux.
4. Apply the δ-mapped slope limiter to suppress oscillations near discontinuities.
5. Advance in time via SSP RK$(k+1)$ under the CFL constraint.

This process achieves high-order accuracy, preserves nontrivial steady states exactly, and suppresses non-physical oscillations arising from flux discontinuities.

## 7. Representative Applications

The RKDG + δ-mapping methodology has been successfully applied to:

- **Nonlinear Elasticity in Layered Media**: For systems
  \[
  \varepsilon_t - v_x = 0, \quad (\rho v)_t - \sigma(\varepsilon, K)_x = 0, \quad \sigma(\varepsilon, K) = K \varepsilon + \beta K^2 \varepsilon^2,
  \]
  with periodically varying $(\rho, K)$, incident pulses split into amplitude-dependent solitary waves; the scheme robustly captures wave speed and structure in agreement with analytic theory.

- **Multi-Class Inhomogeneous Traffic Flow**: For multiple vehicle classes $u_l = a(x) \rho_l$, with flux $f_l = b_l(x)\, u_l\, v(\rho)$, the δ-mapping equalizes supply/demand at interface jumps (lane number/speed limit), cleanly resolving up to $m+2$ nonlinear waves (shocks, rarefactions, contact discontinuities).

## 8. Generalization and Extensions

The δ-mapping concept can be integrated into any classical high-order shock-capturing scheme (WENO, TVD finite difference, central schemes), by replacing cell values $u_i$ by mapped values at discontinuous interfaces. This approach extends naturally to complex network problems, higher-order model couplings, and multidimensional domains with junction or interface conditions, providing a systematic route to high-order, entropy-consistent numerical simulation in the presence of spatially discontinuous coefficients.

---

For further technical details, including explicit construction of δ-mapped states, entropy constraints, and specific numerical benchmarks, see [1511.01346].

Source: https://www.emergentmind.com/topics/runge-kutta-discontinuous-galerkin-schemes