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

In 1948 the Dutch physicist J. M. Burgers tossed out a one-line equation as a model for turbulence. A smooth curve breaks into a vertical cliff in finite time, with no external force — the shortest drama in nonlinear hyperbolic PDEs. This post reads off the shock-formation time by hand from the geometry of characteristics (the spacetime curves along which information travels), shows how a weak solution (defined integrally where differentiability fails) cuts the multivalued region with a single line, and reproduces the same picture in 60 lines of Python. At the end you can grab the (x,t) plane and watch lines collide for yourself.

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.

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.

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#

Use the simulation below to play with amplitude and wavenumber.

Push amplitude $A$ to 0.45 and $t^*$ shrinks to about 0.35; the red dashed line moves up. Set wavenumber $k = 2$ and breaking happens at two locations at once — once formed, the two shocks merge into a larger one. The cyan curve and the faint dashed initial curve enclose equal areas on each side of the shock. That is the equal-area rule, made visible.

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.

Next Time Traffic Slows Down#

Cars merge into one lane and slow. The car ahead is slow, the car behind is fast. The instant the fast car touches the slow car's bumper, the head of the jam is born. That is the meeting point of two trajectories in the (x,t) plane. A traffic jam is essentially a Burgers shock, and it propagates upstream at the Rankine–Hugoniot speed until something resolves.

One-line takeaway. Nonlinear hyperbolic PDEs create discontinuities in finite time even from smooth initial data. Accepting that fact and recovering meaning through the detour of a weak solution is where the numerical theory of shocks begins.