When Lax–Wendroff Shook at a Step — Building a 2nd-Order TVD Advection Scheme with Flux Limiters

In a test case, odd bumps appeared right before and right after a step. I was solving the linear advection equation (pure translation with a constant velocity field) with Lax–Wendroff, and the initial condition was a simple rectangular pulse. As a 2nd-order scheme, I expected it to be sharper than 1st-order upwind — instead the result had an undershoot in front of the step and an overshoot behind it. Logs showed CFL = 0.5, well within the stability limit.

The problem was not that the algorithm was unstable. The algorithm was not monotone.

The limiter buttons start on lax-wendroff. Leave the CFL slider at 0.5 and the orange curve grows bumps around the step. Switch to superbee and those bumps leave.

Orange is the numerical solution. The cyan dashed curve is the exact one. minmod stays dull on both sides, with no bump. superbee can look sharper than the original step. Push CFL near 0.9 and even the quiet limiters spread. Stability does not cancel numerical diffusion.

The Source of the Oscillation — Godunov's Wall#

The 1D linear advection equation is

With constant $u$, the solution is simply the initial profile shifted to the right. Numerically we define a flux $f_{i-1/2}$ at each cell face and write

This is the flux-conserving form, and it preserves conserved quantities (mass, momentum, energy) down to machine precision. Lax–Wendroff is a linear scheme that adds a 2nd-order Taylor correction to this form.

But a theorem Godunov proved in 1959 is unforgiving: no linear scheme for the linear constant-coefficient advection equation can be both monotone and 2nd-order (or higher) accurate. You must give up one of the two. Lax–Wendroff chose 2nd-order accuracy and sacrificed monotonicity — hence the oscillations at the step.

The TVD Condition: Harten's Inequality#

Harten (1983) offered a way out: a scheme whose "total variation decreases".

The left-hand side measures the total wiggle of the numerical solution. If this never grows, no new extrema (overshoot, undershoot) can appear. Schemes that satisfy this are called TVD.

To go around Godunov's wall, the scheme must become nonlinear — it must change shape based on how smooth the solution is locally. 2nd order in smooth regions, 1st-order upwind near sharp gradients. The switch is controlled by a flux limiter.

Flux Limiter in One Line#

The Lax–Wendroff flux is a 1st-order upwind flux plus a correction:

Here $u$ is the advection velocity, $\Delta x$ is the cell width, $\varphi(r)$ is the limiter, and $r$ is the ratio of neighboring slopes.

Smooth solutions give $r \approx 1$; near extrema $r < 0$ or $r \gg 1$. The limiter watches $r$ and dials the correction up or down.

Special choices:

• $\varphi(r) = 0$: donor-cell (1st-order upwind)
• $\varphi(r) = 1$: Lax–Wendroff
• $\varphi(r) = r$: Beam–Warming

In Sweby's diagram, any limiter satisfying $0 \le \varphi(r) \le \min(2r,\,2)$ is TVD. Four representative limiters inside that trapezoid. Setting $\varphi(r) = 0$ is plain donor-cell upwind. What that upwind dulls, in trade for killing the wiggle, is spelled out in the swap when FTCS starts to ring.

Comparing Them in Python#

With roll-based vectorization, this fits in 40 lines. 200 cells, $\mathrm{CFL} = 0.5$, 500 steps (periodic BC; the wave travels around the domain about 2.5 times).

import numpy as np
 
def limiter_phi(name, r):
    r = np.where(np.isfinite(r), r, 0.0)
    if name == 'donor':    return np.zeros_like(r)
    if name == 'lw':       return np.ones_like(r)
    if name == 'minmod':   return np.maximum(0, np.minimum(1, r))
    if name == 'superbee':
        return np.maximum.reduce([np.zeros_like(r),
                                  np.minimum(2*r, 1),
                                  np.minimum(r, 2)])
    if name == 'vanleer':  return (r + np.abs(r)) / (1 + np.abs(r) + 1e-30)
    if name == 'mc':
        return np.maximum(0, np.minimum.reduce([0.5*(1+r),
                                                2*np.ones_like(r),
                                                2*r]))
 
def advect_limited_fvm(q, cfl, name):
    """1D FV advection, periodic BC, u = +1."""
    dq  = q - np.roll(q, 1)                  # face i-1/2: q_i - q_{i-1}
    r   = np.roll(dq, 1) / (dq + 1e-30)      # (q_{i-1}-q_{i-2}) / (q_i-q_{i-1})
    phi = limiter_phi(name, r)
    flx = np.roll(q, 1) + 0.5 * phi * (1 - cfl) * dq
    return q - cfl * (np.roll(flx, -1) - flx)
 
N, cfl, steps = 200, 0.5, 500
x  = np.arange(N) / N
q0 = ((x > 0.2) & (x < 0.4)).astype(float)
for name in ['lw', 'minmod', 'superbee', 'vanleer', 'mc']:
    q = q0.copy()
    for _ in range(steps):
        q = advect_limited_fvm(q, cfl, name)
    l1 = np.abs(q - q0).sum() / N
    print(f'{name:<9s}  L1={l1:.4f}  overshoot={q.max()-1:.3f}  undershoot={q.min():.3f}')

Representative results:

Only lax-wendroff shows a near-20% overshoot; the others are zero. The error ranking is also interesting: superbee has the smallest $L_1$ error, meaning it keeps the step crisp — at the cost of turning smooth sine waves into staircases.

If you repeat the experiment with a smooth wave ($q_0 = \sin(2\pi x)$, 500 steps), the convergence order drops to about 1.5 for minmod, 2.0 for van Leer/MC, and 1.3 for superbee. The "2nd-order accurate" label slips every time the limiter activates near an extremum and locally falls back to 1st order.

At r=1 the four limiters return the same phi#

Plug $r=1$ into the table. minmod, superbee, van Leer, and MC all return 1. $\varphi=1$ is Lax–Wendroff. In the smooth middle they apply the same correction. They split only when $r$ leaves 1.

Try $r=2$. minmod stops at 1. superbee sits on 2, the top edge of Sweby. van Leer is $4/3$. MC is 1.5. That top edge is why superbee keeps the step $L_1$ at 0.042, the smallest in the table. The same choice drops the sine-wave order to about 1.3. van Leer and MC stay near 2.0. minmod lands near 1.5. At an extremum, $r \le 0$, all four fall to 0 and become donor-cell.

The correction coefficient is already fixed by the code. With speed 1 and $\mathrm{CFL}=0.5$, the factor in front of the jump is 0.25. Even with $\varphi=1$, only a quarter of the neighbor jump is added. Five hundred steps travel $0.5\times 500=250$ cells. On a grid of 200 that is 1.25 trips around the periodic domain. The overshoot +0.181 is 18 percent of a unit pulse. The undershoot −0.133 is how far the floor goes negative. minmod's $L_1$ of 0.081 is about 1.9 times superbee's 0.042. Killing the bump costs that much extra smear.

Practical Checklist#

• For shock capturing, start with minmod or van Leer. Zero oscillation is the priority.
• For LES / large-scale turbulence, superbee is risky. It chops smooth turbulent structures into staircases and distorts the energy spectrum. van Leer or MC is the compromise.
• Dropping to 1st-order upwind at boundary cells because the stencil is short drops global accuracy to 1st order. Secure the limiter stencil at boundaries or use ghost cells.
• OpenFOAM's limitedLinear 1 is a one-parameter limiter in the van Leer family. Coefficient 0 = 1st-order upwind; 1.0 = full TVD region.
• Always verify convergence order with a smooth solution. Measuring $L_1$ slopes on a test case that contains a discontinuity pins the order down — because of the discontinuity, not the limiter.

superbee on a step drops the sine wave to order 1.3#

In the table superbee's $L_1$ is 0.042, the smallest. The overshoot is +0.002, almost invisible on the step. Put the same limiter on a sine wave and the order falls to about 1.3. van Leer and MC stay near 2.0. The choice that protects the jump spends the smooth wave's order.

If the scalar is a density, the lax-wendroff undershoot of −0.133 is the door to a negative density. The limiter that closes that door is written up in where a positive value leaks. For a passive step, minmod or van Leer is enough. For a density, open that next door.