When Characteristics Collide, a Shock Is Born — Burgers and the (x,t) Plane

In 1948 Burgers tossed out a one-line model of turbulence. A smooth curve becomes a cliff in finite time, with no external force. The birthday of that cliff is the instant characteristics (spacetime curves that carry information) cross. A weak solution — the integral form that survives after differentiability dies — cuts the overlap with one line.

Raise amplitude $A$ or wavenumber $k$ and the red dashed breaking time moves up.

That dashed line is the first time a crest line steps on the valley line ahead. Larger $A$ pulls the meeting forward. At $k = 2$ the crossings come in pairs.

Two Lines Meet at One Point#

Hyperbolic equations are about signal propagation. As long as information flows along characteristics at fixed speeds, life is peaceful. For linear advection $\partial_t q + u\,\partial_x q = 0$ with constant $u$, the characteristics are parallel straight lines and the initial profile simply slides sideways. There is no room for a shock in this picture.

Nonlinearity changes everything. When $u$ depends on $q$ itself, large values move fast and small values move slow. The instant a fast piece overtakes a slow piece, two values of $q$ arrive at the same $(x,t)$ point. The function becomes multi-valued — physically nonsense — and that nonsense is the birth of a shock.

A Self-Advecting Equation#

The conservative form of Burgers' equation reads:

$q$ is the conserved scalar and $f(q) = q^2/2$ is the flux. The chain rule gives the non-conservative form:

The advection speed is $u(x,t) = q(x,t)$. The quantity transports itself. That single nonlinearity is the source of every drama to follow. The Jacobian $\partial f/\partial q = q$ is positive where signals run rightward and negative where they run leftward. An initial profile whose sign varies in space — a sine wave is the canonical example — forces fast points to chase slow ones.

The (x,t) Plane Drawn by Characteristics#

A characteristic is a curve obeying $dx/dt = q$. Since the comoving derivative $Dq/Dt$ vanishes on it, $q$ is frozen at its starting value along that curve. The characteristic is a straight line in the (x,t) plane with slope $1/q_0(x_0)$; the starting value $q_0(x_0)$ is the line's velocity.

For the initial condition $q_0(x_0) = 0.55 + A\sin(2\pi k x_0)$:

• Lines launched from peaks ($q_0 > 0.55$) tilt sharply rightward (fast).
• Lines launched from valleys ($q_0 < 0.55$) tilt gently (slow).

A peak line will eventually step on the toe of the valley line right ahead of it. The point where the two lines meet is the first shock.

Breaking Time t* — The First Collision#

Two neighboring characteristics from $x_0$ and $x_0 + \Delta x$ meet at time $t$ when:

Taking the limit $\Delta x \to 0$:

Only negative slopes produce a meeting. The earliest collision happens where $q_0'(x_0)$ is most negative:

Larger amplitude or higher wavenumber breaks faster. For the sine above, $q_0'(x_0) = 2\pi k A\cos(2\pi k x_0)$, with minimum $-2\pi k A$:

With $A = 0.35$, $k = 1$, $t^* \approx 0.455$. Before that moment the curve smoothly skews. The instant $t^*$ is crossed, a vertical wall stands up.

The most negative slope is $-2\pi k A$. Since $2\pi \approx 6.28$, $2\pi \times 0.35 \approx 2.20$, and the reciprocal is $0.455$. The crest is $0.55 + 0.35 = 0.90$ and the valley is $0.20$, so the initial jump is $0.70$. That matches the later shock height of about $0.7$. Out to $t = 1.2$ the jump does not grow. Only the corners round off.

Push $A$ to $0.45$ and the denominator rises to about $2.83$, so $t^*$ drops near $0.35$. Set the wavenumber to $2$ as well and the denominator doubles, $t^* \approx 0.227$. Same amplitude, shorter wavelength, earlier crossing.

Weak Solutions and Rankine–Hugoniot — The Knife That Cuts#

Past $t^*$, the straight reading hands you three values of $q$ at one $x$ — impossible. The fix is a weak solution: weaken the definition by dropping the differential form for the integral conservation law.

This makes sense even when $q$ is discontinuous. If two flat states $q_L$ and $q_R$ are joined by a shock moving at speed $s$:

This is the Rankine–Hugoniot condition. The shock speed is the arithmetic mean of the two states. The condition slices the multivalued region with a single line and arranges it so the two cut areas are equal — the equal-area rule.

A separate entropy condition is needed too. A shock is physical only when $q_L > q_R$. If $q_L < q_R$, the answer is a fan-shaped rarefaction. Violating this leaves you with a mathematically valid but unphysical expansion shock.

Plug $0.90$ and $0.20$ into $s = (q_L + q_R)/2$ and $s = 0.55$, the background of the sine. The shock slides at that mean speed and leaves only the two states. Reverse the order and the speed formula still returns $0.55$. Only the answer changes, from shock to rarefaction. One arithmetic mean does not tell the two waves apart.

The mesh spacing is $1/256 \approx 0.00391$. With CFL $0.4$ and the initial $\max|q| = 0.90$, $\Delta t \approx 0.00174$. Reaching $t = 1.2$ takes about $690$ steps. After the shock, a local spike in $\max|q|$ makes the same formula chop the step smaller. That knife at each face is the same one as in a shock tube in a single cell.

Python: From Tracing Characteristics to Godunov Shocks#

Trace characteristics and capture the weak shock with a first-order Godunov scheme.

import numpy as np
 
def burgers_breaking_time(q0_prime_min):
    # q0_prime_min: minimum of q0'(x) over x (must be negative for a shock)
    if q0_prime_min >= 0:
        return float('inf')
    return -1.0 / q0_prime_min
 
def trace_characteristics(q0_func, x0_arr, t_arr):
    # x[i, j] = x0_arr[i] + q0(x0_arr[i]) * t_arr[j]
    return x0_arr[:, None] + q0_func(x0_arr)[:, None] * t_arr[None, :]
 
def rankine_hugoniot_speed(qL, qR):
    # shock speed s = (qL + qR) / 2
    return 0.5 * (qL + qR)
 
def godunov_burgers_step(q, dx, dt):
    # conservation form: q^{n+1} = q^n - (dt/dx) * (F_{i+1/2} - F_{i-1/2})
    qL = q                          # left state at right interface
    qR = np.roll(q, -1)             # right state at right interface
    s = 0.5 * (qL + qR)
    f_shock = np.where(s >= 0, 0.5 * qL**2, 0.5 * qR**2)
    f_rare  = np.where(qL >= 0, 0.5 * qL**2,
              np.where(qR <= 0, 0.5 * qR**2, 0.0))
    F_iph = np.where(qL > qR, f_shock, f_rare)   # interface flux on the right
    F_imh = np.roll(F_iph, 1)                    # left interface = shift
    return q - (dt / dx) * (F_iph - F_imh)
 
# Simulation parameters --------------------------------------
N, L = 256, 1.0
A, k = 0.35, 1
x = (np.arange(N) + 0.5) * (L / N)
q0 = lambda xx: 0.55 + A * np.sin(2 * np.pi * k * xx)
q  = q0(x)
 
t_star = burgers_breaking_time(-2 * np.pi * k * A)
print(f"breaking time t* = {t_star:.3f}")
 
t_end, t = 1.2, 0.0
while t < t_end:
    dt = 0.4 * (L / N) / np.max(np.abs(q))    # CFL = 0.4
    q  = godunov_burgers_step(q, L / N, dt)
    t += dt
print(f"shock height (max - min) at t = {t:.2f}: {q.max() - q.min():.3f}")

You should see breaking time t* = 0.455 and a shock height of about 0.7. The Godunov scheme automatically slices the multivalued region into a single line — the knife is the line np.where(qL > qR, f_shock, f_rare). That one expression solves the exact Riemann problem at every cell interface and guarantees the weak solution.

Watch the Lines Collide#

In the figure above, the cyan curve and the faint initial dashes enclose equal areas on each side of the shock. That is the equal-area rule. At $k = 2$ the two crossings soon merge, and the jump grows.

Equal areas mean that, after the cut, $\int q\,dx$ still matches the flux difference at the two ends. If the shock walks at the mean of $0.90$ and $0.20$, which is $0.55$, the two clipped triangles cancel. The height stays $0.70$. Only the location moves at $s = 0.55$.

Two Reasons Numerical Shock Capture Fails#

(1) CFL violation. Break $\Delta t \le \Delta x / \max|q|$ and the numerical domain outruns the physical one. If amplitude appears to grow with time, suspect $dt$ first. The code above recomputes $dt$ every step from np.max(np.abs(q)), which matters even more once a shock forms because $q$ can spike locally.

(2) Wrong entropy. Rankine–Hugoniot alone admits unphysical expansion shocks. Simple linearized solvers like Roe make an "entropy glitch" in transonic rarefactions (where the fan crosses zero). Godunov, HLLE, and exact Riemann solvers enforce the entropy condition automatically, but Roe needs an explicit Harten-style entropy fix.

When a compressible LES code shows shock positions oscillating or negative pressures, almost always it is one of these two.

Whether a jam counts as a shock#

Cars pack into one lane and slow down. The one ahead is slow, the one behind is fast. The bumper contact is a crossing in the (x,t) plane. That head is a Burgers shock, and until it clears it drifts upstream at the Rankine–Hugoniot speed.

Hold the jam as a continuous curve and the multi-valued fold stays. Whether to cut it with one line or open it as a rarefaction is decided by which state is larger. If a step starts ringing, the first place to look is a flux limiter.