Limiters, the TVD Condition, and Shock Capture

A first-order upwind scheme smears a contact discontinuity over many cells and provides no way to sharpen it. A second-order scheme preserves the sharpness and adds oscillations around the jump that do not decay as the grid refines: the computed solution has undershoots below the initial minimum and overshoots above the maximum. Both behaviors follow from the same analysis, and the result that reconciles them imposes a precise constraint on what a scheme may do locally.

The oscillation and where it comes from

Interpolate cell averages with a linear profile and the reconstruction is second order. Its truncation error has a leading term proportional to the second derivative, and at a discontinuity the second derivative is a delta function: the reconstruction overshoots on one side and undershoots on the other. The amplitude of the overshoot does not shrink as the grid refines, because the discontinuity remains a discontinuity. This is the Godunov barrier in its empirical form, and it has a theorem behind it.

A linear scheme; constant coefficients, the same stencil everywhere; that is monotone, in the sense of preserving the ordering of initial data, can be at most first order accurate in the maximum norm. The statement is usually quoted as a consequence of Godunov’s theorem: there is no linear scheme of second or higher order that is monotone. Monotonicity is desirable because it forbids new extrema, and forbidding new extrema is exactly what eliminates the oscillation at a discontinuity. The barrier says that linearity and higher order and monotonicity cannot coexist, so one must go.

Total variation and what the TVD condition buys

Relax monotonicity to a weaker but still useful property. Define the total variation of a discrete solution as

$$ TV(u) = \sum_j \left| u_{j+1} - u_j \right| $$

For a scalar conservation law with convex flux, the exact solution is total-variation diminishing: $TV(u(t))$ does not increase with time. A numerical scheme is TVD if $TV(u^{n+1}) \le TV(u^n)$ for all admissible data. The property forbids new extrema, since a new local maximum would require an additional rise and fall, which increases the total variation. It also provides a compactness argument for convergence of the scheme’s sequence of grid solutions.

The condition is checked through a local form. Harten’s sufficient condition for a three-point explicit scheme $u_j^{n+1} = u_j^n - C_{j-1}\left(u_j^n - u_{j-1}^n\right) + D_j\left(u_{j+1}^n - u_j^n\right)$ is that the coefficients satisfy $C \ge 0$, $D \ge 0$ and $C + D \le 1$ for every cell. The condition is sufficient, not necessary, and it is the workhorse of limiter design because it reduces to inequalities on the coefficients that can be enforced by construction.

For a scheme built from a reconstruction, the condition becomes a condition on the slope. Write the reconstruction as $u_j + \sigma_j(x - x_j)$ and the update as the first-order upwind flux plus a correction involving the slopes. The correction is TVD if the slope is bounded by twice the smaller of the one-sided slopes, which is the core of every slope limiter.

The limiter function and its admissible region

The standard formulation introduces a limiter $\phi(r)$ applied to the ratio of consecutive slopes. With $r$ the ratio of the upwind to the downwind difference,

$$ r_j = \frac{u_j - u_{j-1}}{u_{j+1} - u_j} $$

the reconstructed slope is $\sigma_j = \phi(r_j),(u_{j+1} - u_j)$, and the TVD condition is that $\phi$ satisfy, for all $r > 0$,

$$ 0 \le \phi(r) \le 2, \qquad 0 \le \phi(r) \le 2r $$

The region is the intersection of the strip $0 \le \phi \le 2$ with the wedge $\phi \le 2r$. Within it, each named limiter makes a different choice. Minmod takes $\phi(r) = \max(0, \min(1, r))$, the most restrictive option, which reverts to first order wherever the slopes have opposite signs and produces a smooth but diffusive result. Superbee takes the upper envelope, $\phi(r) = \max(0, \min(2r, 1), \min(r, 2))$, which is the least diffusive TVD limiter and compresses discontinuities toward a single cell, sometimes at the cost of turning a smooth profile into steps. Van Leer’s harmonic limiter, $\phi(r) = (r + |r|)/(1 + |r|)$, sits between them and is smooth in $r$, which makes it convenient for implicit schemes where a nondifferentiable limiter upsets the Newton iteration discussed in Newton iterations for implicit operators.

The second-order accuracy of the limited scheme is not uniform. In smooth regions $r \approx 1$ and every admissible limiter returns $\phi(1) = 1$, so the scheme is second order there. At extrema $r \le 0$ and the limiters return zero, degrading the scheme to first order in a neighborhood of every smooth extremum. That degradation clips the peaks of smooth profiles. A sine wave advected by a TVD scheme loses amplitude at its maxima, and the loss is a genuine defect of the TVD class, not an implementation artifact.

Implementing a limited reconstruction

The flux form is the natural implementation. For a scalar advection with speed $c > 0$, each face between cells $j-1$ and $j$ receives the upwind cell’s reconstructed value, and the face flux is the speed times that value:

/* MUSCL reconstruction with minmod limiter, upwind flux, c > 0. */
double slopes[N];
for (size_t j = 1; j + 1 < N; j++) {
    double d_up   = u[j] - u[j - 1];        /* upwind difference */
    double d_down = u[j + 1] - u[j];        /* downwind difference */
    slopes[j] = (d_up * d_down > 0.0)       /* same sign: limit */
              ? copysign(fmin(fabs(d_up), fabs(d_down)), d_up)
              : 0.0;                        /* extremum: fall back to first order */
}
for (size_t j = 1; j < N; j++) {
    /* value on the right edge of the upwind cell j-1 */
    double u_face = u[j - 1] + 0.5 * slopes[j - 1];
    flux[j] = c * u_face;
}

Two details govern the result. The slope must be computed per cell and stored, because each face uses the reconstruction of the cell on its upwind side; recomputing it at the face from the two adjacent values would lose the limiter’s information and reduce the scheme to the centered stencil. And the limiter must be applied to the upwind cell’s slope, not to an average of the two cells, or the scheme stops being upwind and the TVD argument does not apply to it.

The Courant condition from the linear analysis, $\nu = c\Delta t/\Delta x \le 1$, still bounds the step, and for the limited scheme the bound is unchanged for the linear case. The limiter constrains the reconstruction, not the time integration. A scheme that runs a limited reconstruction with a multistage time integrator must satisfy the bound for the integrator’s own stability region rather than for the first-order update, since a stage evaluates the flux at an intermediate state; the stability region analysis for those methods appears in Butcher tableaux and order conditions.

Where TVD stops being the right answer

Four limitations of the class are worth naming, because each motivates a different extension.

The peak clipping described above follows directly from the second-order degeneracy at extrema, and it is what motivates the total variation bounded (TVB) modification, which relaxes the constraint where the solution is smooth enough that the local variation is $O(\Delta x^2)$. The relaxation restores second-order accuracy at smooth extrema at the cost of a solution-dependent threshold, and the choice of threshold is a tuning parameter with consequences for both accuracy and resilience.

The second is that the TVD property is a statement about scalar problems. For systems, applying a limiter component by component to the conserved variables does not preserve the crucial property that the reconstruction stays in the physical state space, and it does not respect the characteristic structure. The fix is to limit in characteristic variables, transforming to the local eigenbasis, limiting each wave, and transforming back. The cost is a small matrix operation per face and per step; the benefit is correct behavior at shocks where the conserved variables change by different amounts.

The third is multidimensionality. The TVD condition does not extend to finite-volume schemes on unstructured meshes in a way that preserves accuracy near discontinuities, and genuinely multidimensional limiters or the alternative class of genuinely multidimensional schemes must be used. A common compromise is a limiter applied along each face’s normal direction, which is dimensionally split and inherits the one-dimensional properties only approximately.

The fourth is accuracy in smooth regions. TVD is a strong constraint, and a scheme that merely satisfies it will be diffuse. The class of weighted essentially non-oscillatory (WENO) schemes relaxes the constraint to a weaker one: instead of forbidding new extrema, it detects non-smooth stencils and weights them down, obtaining higher order in smooth regions while keeping the oscillation bounded at discontinuities. The design trade shifts from a hard bound to a weighting function that must be chosen with care, and the price is more stencil evaluations per face.

The connection to the hyperbolic step limit

Everything above presupposes that the time integration does not itself introduce oscillations. The stability condition for the explicit time step is the CFL bound derived in global error accumulation and the CFL condition, and it is independent of the limiter. When a scheme is run at a Courant number near 1 with a limiter that is too aggressive, the limiter’s first-order behavior spreads over the whole profile and the result is not a sharp discontinuity but a smeared one; when the limiter is too permissive, the oscillations return. The scheme’s behavior is the joint result of the reconstruction, the limiter, and the step, and each must be validated against the others.

The final observation is about steady states. A limiter enforces the TVD condition along the marching direction, and a steady-state computation reaches the same condition after many time steps. The steady solution’s accuracy near a shock is then governed by the limiter, and the choice between minmod, superbee, and the harmonic limiter is a choice about how many cells the shock occupies in the converged answer. That choice is the practical content of the trade the TVD condition leaves open.