Chapter 10.3: LWR Theory & Spinodal Decomposition

From Lattice Gas to Continuum: The LWR Equation

The TASEP mean-field equation, when taken to the continuum limit, yields the Lighthill-Whitham-Richards (LWR) equation—the foundational PDE of macroscopic traffic flow theory. This is a scalar conservation law whose characteristics are the kinematic waves of traffic. When the fundamental diagram is concave, shocks form spontaneously: free-flowing traffic suddenly transitions to a jam. The analogy with spinodal decomposition in phase-separating mixtures provides deep insight into why traffic jams have a preferred spacing.

10.3.1 Deriving LWR from TASEP Mean-Field

The TASEP mean-field equation for site $i$ is:

$$\frac{d\rho_i}{dt} = \rho_{i-1}(1-\rho_i) - \rho_i(1-\rho_{i+1})$$

This is a discrete conservation law: $d\rho_i/dt = J_{i-1,i} - J_{i,i+1}$ with current $J_{i,i+1} = \rho_i(1-\rho_{i+1})$. Taking the continuum limit with lattice spacing$a$ and $x = ia$:

$$\frac{\partial \rho}{\partial t} + \frac{\partial q(\rho)}{\partial x} = 0$$

This is the LWR equation (Lighthill & Whitham, 1955; Richards, 1956), with the flux function $q(\rho) = \rho v(\rho)$. For the Greenshields model:

$$q(\rho) = v_{\max} \rho\!\left(1 - \frac{\rho}{\rho_{\max}}\right)$$

Characteristics and Wave Speed

The characteristic speed is:

$$c(\rho) = q'(\rho) = v_{\max}\!\left(1 - \frac{2\rho}{\rho_{\max}}\right)$$

Information propagates forward in free flow ($\rho < \rho_c$) and backward in congestion ($\rho > \rho_c$). Since $q''(\rho) = -2v_{\max}/\rho_{\max} < 0$, the flux is concave: characteristics from higher density travel slower, meaning they converge and shocks form inevitably.

10.3.2 Shock Waves: Rankine-Hugoniot Condition

When characteristics converge, a shock (discontinuity in density) forms. The shock speed is given by the Rankine-Hugoniot jump condition, derived from conservation of vehicles across the shock:

$$s = \frac{q(\rho_R) - q(\rho_L)}{\rho_R - \rho_L} = \frac{\Delta q}{\Delta \rho}$$

where $\rho_L, \rho_R$ are the densities to the left and right of the shock. For the Greenshields model:

$$s = v_{\max}\!\left(1 - \frac{\rho_L + \rho_R}{\rho_{\max}}\right)$$

Lax Entropy Condition

Not every discontinuity is a physically admissible shock. The Lax entropy conditionrequires that characteristics must converge onto the shock from both sides:

$$c(\rho_L) > s > c(\rho_R)$$

For a concave flux ($q'' < 0$), this means$\rho_L < \rho_R$: density increases across a shock. A transition from high density to low density would instead be a rarefaction wave(a smooth expansion fan). This is why traffic jams have a sharp upstream front (shock) and a gradual downstream tail (rarefaction).

Physical Interpretation

A driver approaching a jam encounters a sudden wall of stopped traffic (the shock). After passing through the jam and reaching the front, traffic gradually accelerates back to free-flow speed (the rarefaction). This asymmetry is a direct consequence of the Lax entropy condition.

10.3.3 Payne-Whitham: Second-Order Traffic Model

The LWR equation is first-order: velocity is instantaneously determined by density. The Payne-Whitham model adds a velocity relaxation equation, creating a second-order system:

$$\frac{\partial \rho}{\partial t} + \frac{\partial (\rho v)}{\partial x} = 0$$

$$\frac{\partial v}{\partial t} + v \frac{\partial v}{\partial x} = \frac{V_e(\rho) - v}{\tau_r} - \frac{c_0^2}{\rho}\frac{\partial \rho}{\partial x}$$

The terms on the right represent: (1) relaxation toward the equilibrium velocity$V_e(\rho)$ with time scale$\tau_r$, and (2) an anticipation pressure $c_0^2$ representing drivers reacting to density gradients ahead of them.

The key insight: when the effective diffusion coefficient becomes negative, small perturbations grow rather than decay. Linearising around uniform flow$(\rho_0, v_0)$:

$$D_{\text{eff}} = c_0^2 - \rho_0 |V_e'(\rho_0)| \cdot v_0 \, \tau_r$$

When $D_{\text{eff}} < 0$, the uniform state is unstable—this is the traffic analog of spinodal decomposition.

10.3.4 Cahn-Hilliard Analogy and Spinodal Instability

Adding a fourth-order regularisation to the LWR equation (analogous to surface tension in phase separation) yields the Cahn-Hilliard form:

$$\frac{\partial \rho}{\partial t} + \frac{\partial q}{\partial x} = \frac{\partial}{\partial x}\!\left(D \frac{\partial \rho}{\partial x} - \kappa \frac{\partial^3 \rho}{\partial x^3}\right)$$

where $D = q''(\rho)/2$ is the effective diffusion coefficient and $\kappa > 0$ is the regularisation (viscosity/anticipation). When$q''(\rho) < 0$ (which holds for the Greenshields model everywhere), $D < 0$ and the system is in the spinodal region.

Linear Stability Analysis

Perturbing around uniform density $\rho_0$with $\rho = \rho_0 + \delta\rho \, e^{ikx + \sigma t}$:

$$\sigma(k) = -D k^2 - \kappa k^4 = |D| k^2 - \kappa k^4$$

The growth rate $\sigma(k)$ is positive for wavenumbers $0 < k < k_c = \sqrt{|D|/\kappa}$. The most unstable (fastest growing) wavenumber is:

$$k^* = \frac{k_c}{\sqrt{2}} = \sqrt{\frac{|D|}{2\kappa}} \quad \Longrightarrow \quad \lambda^* = 2\pi\sqrt{\frac{2\kappa}{|D|}}$$

Predicted Jam Spacing

For typical traffic parameters ($v_{\max} \approx 120$ km/h,$\rho_{\max} \approx 150$ veh/km,$\kappa \approx 0.01$ km³/h), the predicted most unstable wavelength is:

$$\lambda^* \approx 1\text{--}3 \text{ km}$$

This matches the empirically observed spacing between stop-and-go traffic jams on highways remarkably well. The Cahn-Hilliard analogy explains why jams have a characteristic spacing: just as phase-separating alloys form domains of a preferred size, traffic instabilities select a preferred wavelength through the competition between the destabilising negative diffusion and the stabilising higher-order viscosity.

10.3.5 Simulation: Godunov Scheme for LWR

We solve the LWR equation numerically using the Godunov scheme, which respects the conservation law structure and correctly captures shocks. The Godunov flux for a concave flux function is:

$$F_{i+1/2} = \begin{cases} \min_{\rho_L \leq \rho \leq \rho_R} q(\rho) & \text{if } \rho_L \leq \rho_R \\ \max_{\rho_R \leq \rho \leq \rho_L} q(\rho) & \text{if } \rho_L > \rho_R \end{cases}$$

LWR Godunov Scheme: Shock Formation & Spinodal Instability

Python

script.py160 lines

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('Agg')

# ── LWR with Greenshields using Godunov scheme ──
rho_max = 150.0   # veh/km
v_max = 120.0     # km/h

def q_flux(rho):
    """Greenshields flux."""
    return v_max * rho * (1 - rho / rho_max)

def dq_drho(rho):
    """Characteristic speed."""
    return v_max * (1 - 2 * rho / rho_max)

def godunov_flux(rho_L, rho_R):
    """Godunov numerical flux for concave q."""
    rho_c = rho_max / 2.0
    if rho_L <= rho_R:
        # Rarefaction or no wave
        if rho_L >= rho_c:
            return q_flux(rho_L)
        elif rho_R <= rho_c:
            return q_flux(rho_R)
        else:
            return q_flux(rho_c)  # sonic point
    else:
        # Shock
        if rho_L <= rho_c:
            return q_flux(rho_L)
        elif rho_R >= rho_c:
            return q_flux(rho_R)
        else:
            return max(q_flux(rho_L), q_flux(rho_R))

# ── Setup ──
L_road = 20.0     # km
Nx = 500
dx = L_road / Nx
x = np.linspace(dx/2, L_road - dx/2, Nx)
CFL = 0.8
T_final = 0.15    # hours (~9 minutes)

# Initial condition: traffic signal creates a jam
rho = np.where((x > 8) & (x < 12), 0.8 * rho_max, 0.2 * rho_max)
# Add small sinusoidal perturbation
rho += 5.0 * np.sin(2 * np.pi * x / 2.5)
rho = np.clip(rho, 0.01, rho_max - 0.01)

# Store snapshots
snapshots = [(0.0, rho.copy())]
t = 0.0
dt_max = CFL * dx / v_max

while t < T_final:
    dt = min(dt_max, T_final - t)

# Compute Godunov fluxes
    F = np.zeros(Nx + 1)
    for i in range(Nx + 1):
        rL = rho[max(i-1, 0)]
        rR = rho[min(i, Nx-1)]
        F[i] = godunov_flux(rL, rR)

# Update
    rho_new = rho - (dt / dx) * (F[1:] - F[:-1])
    rho = np.clip(rho_new, 0.0, rho_max)
    t += dt

if t >= 0.03 and len(snapshots) < 2:
        snapshots.append((t, rho.copy()))
    elif t >= 0.07 and len(snapshots) < 3:
        snapshots.append((t, rho.copy()))
    elif t >= 0.11 and len(snapshots) < 4:
        snapshots.append((t, rho.copy()))

snapshots.append((T_final, rho.copy()))

# ── Plotting ──
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
fig.patch.set_facecolor('#0a0a0a')

colors_t = ['#34d399', '#2dd4bf', '#fbbf24', '#f87171', '#a78bfa']

# Panel 1: Density evolution
ax1 = axes[0, 0]
ax1.set_facecolor('#0a0a0a')
ax1.tick_params(colors='white')
for spine in ax1.spines.values():
    spine.set_color('#334155')
for i, (ti, ri) in enumerate(snapshots):
    ax1.plot(x, ri, color=colors_t[i], linewidth=1.8, label=f't = {ti*60:.1f} min')
ax1.set_xlabel('Position (km)', color='white')
ax1.set_ylabel('Density (veh/km)', color='white')
ax1.set_title('LWR: Shock Formation', color='#6ee7b7', fontsize=13)
ax1.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white', fontsize=9)

# Panel 2: Flow vs density (fundamental diagram with trajectory)
ax2 = axes[0, 1]
ax2.set_facecolor('#0a0a0a')
ax2.tick_params(colors='white')
for spine in ax2.spines.values():
    spine.set_color('#334155')
rho_plot = np.linspace(0, rho_max, 300)
ax2.plot(rho_plot, q_flux(rho_plot), color='#34d399', linewidth=2.5, label='q(ρ)')
ax2.scatter(snapshots[-1][1][::10], q_flux(snapshots[-1][1][::10]),
            color='#f87171', s=10, alpha=0.5, label='Final state')
ax2.set_xlabel('Density (veh/km)', color='white')
ax2.set_ylabel('Flow (veh/h)', color='white')
ax2.set_title('Fundamental Diagram', color='#6ee7b7', fontsize=13)
ax2.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')

# Panel 3: Characteristic speed
ax3 = axes[1, 0]
ax3.set_facecolor('#0a0a0a')
ax3.tick_params(colors='white')
for spine in ax3.spines.values():
    spine.set_color('#334155')
ax3.plot(x, dq_drho(snapshots[0][1]), color='#34d399', linewidth=1.5, alpha=0.6, label='t=0')
ax3.plot(x, dq_drho(snapshots[-1][1]), color='#f87171', linewidth=1.5, label='t=final')
ax3.axhline(y=0, color='#fbbf24', linestyle='--', alpha=0.5)
ax3.set_xlabel('Position (km)', color='white')
ax3.set_ylabel('Wave speed dq/dρ (km/h)', color='white')
ax3.set_title('Characteristic Speed', color='#6ee7b7', fontsize=13)
ax3.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')

# Panel 4: Spinodal instability (growth rate)
ax4 = axes[1, 1]
ax4.set_facecolor('#0a0a0a')
ax4.tick_params(colors='white')
for spine in ax4.spines.values():
    spine.set_color('#334155')

k_arr = np.linspace(0.01, 10, 500)
D_eff = v_max / rho_max  # |D| = |q''|/2
kappa_values = [0.005, 0.01, 0.02, 0.05]
colors_k = ['#34d399', '#2dd4bf', '#fbbf24', '#f87171']

for i, kappa in enumerate(kappa_values):
    sigma = D_eff * k_arr**2 - kappa * k_arr**4
    ax4.plot(k_arr, sigma, color=colors_k[i], linewidth=2, label=f'κ = {kappa}')
    # Mark most unstable
    k_star = np.sqrt(D_eff / (2 * kappa))
    sigma_star = D_eff**2 / (4 * kappa)
    if k_star < 10:
        ax4.plot(k_star, sigma_star, 'o', color=colors_k[i], markersize=6)

ax4.axhline(y=0, color='white', linestyle='--', alpha=0.3)
ax4.set_xlabel('Wavenumber k (1/km)', color='white')
ax4.set_ylabel('Growth rate σ (1/h)', color='white')
ax4.set_title('Spinodal Growth Rate σ(k)', color='#6ee7b7', fontsize=13)
ax4.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white', fontsize=9)

plt.tight_layout(pad=2.0)
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("LWR Godunov simulation with shock formation and spinodal analysis complete.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

10.3.6 Computing Jam Spacing

The most unstable wavelength $\lambda^*$predicts the spacing between stop-and-go waves. We can compute this for different traffic conditions:

Jam Spacing from Spinodal Instability

Python

script.py101 lines

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('Agg')

# ── Jam spacing prediction from Cahn-Hilliard ──
v_max = 120.0    # km/h
rho_max = 150.0  # veh/km

# D_eff = |q''(rho)|/2 depends on rho for general flux
# For Greenshields: q''(rho) = -2*v_max/rho_max = const
D_abs = v_max / rho_max  # |q''|/2

kappa_range = np.linspace(0.001, 0.1, 200)
lambda_star = 2 * np.pi * np.sqrt(2 * kappa_range / D_abs)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
fig.patch.set_facecolor('#0a0a0a')
for ax in [ax1, ax2]:
    ax.set_facecolor('#0a0a0a')
    ax.tick_params(colors='white')
    for spine in ax.spines.values():
        spine.set_color('#334155')

ax1.plot(kappa_range, lambda_star, color='#34d399', linewidth=2.5)
ax1.axhspan(1, 3, alpha=0.15, color='#fbbf24', label='Observed: 1-3 km')
ax1.set_xlabel('Anticipation parameter κ (km³/h)', color='white')
ax1.set_ylabel('Jam spacing λ* (km)', color='white')
ax1.set_title('Predicted Jam Spacing', color='#6ee7b7', fontsize=13)
ax1.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
ax1.set_ylim(0, 8)

# ── Solve Cahn-Hilliard-type traffic equation ──
# Simulate density evolution with negative diffusion + 4th order term
L_road = 30.0
Nx = 600
dx = L_road / Nx
x = np.linspace(dx/2, L_road - dx/2, Nx)

rho_0 = 0.5 * rho_max  # In the unstable regime
kappa = 0.015

# Small random perturbation around uniform density
np.random.seed(42)
rho = rho_0 + 3.0 * np.random.randn(Nx)
rho = np.clip(rho, 5, rho_max - 5)

# Explicit time-stepping with 4th-order term
dt = 0.5 * dx**2 / (4 * D_abs + 4 * kappa / dx**2)
n_steps = 3000

for step in range(n_steps):
    # Second derivative (diffusion)
    d2rho = np.zeros(Nx)
    d2rho[1:-1] = (rho[2:] - 2*rho[1:-1] + rho[:-2]) / dx**2
    d2rho[0] = (rho[1] - 2*rho[0] + rho[-1]) / dx**2  # periodic BC
    d2rho[-1] = (rho[0] - 2*rho[-1] + rho[-2]) / dx**2

# Fourth derivative
    d4rho = np.zeros(Nx)
    d4rho[2:-2] = (rho[4:] - 4*rho[3:-1] + 6*rho[2:-2] - 4*rho[1:-3] + rho[:-4]) / dx**4
    d4rho[0] = d4rho[2]
    d4rho[1] = d4rho[2]
    d4rho[-1] = d4rho[-3]
    d4rho[-2] = d4rho[-3]

# Advective flux derivative (conservative)
    dq = np.zeros(Nx)
    q_val = v_max * rho * (1 - rho / rho_max)
    dq[1:-1] = (q_val[2:] - q_val[:-2]) / (2 * dx)
    dq[0] = (q_val[1] - q_val[-1]) / (2 * dx)
    dq[-1] = (q_val[0] - q_val[-2]) / (2 * dx)

rho = rho - dt * dq + dt * D_abs * 0.5 * d2rho - dt * kappa * d4rho
    rho = np.clip(rho, 1, rho_max - 1)

# Plot final density profile
ax2.plot(x, rho, color='#34d399', linewidth=1.5)
ax2.axhline(y=rho_0, color='#fbbf24', linestyle='--', alpha=0.5, label=f'ρ₀ = {rho_0:.0f}')

# Find peaks to measure spacing
from scipy.signal import find_peaks
peaks, _ = find_peaks(rho, height=rho_0 + 10, distance=20)
if len(peaks) > 1:
    spacings = np.diff(x[peaks])
    mean_spacing = np.mean(spacings)
    ax2.set_title(f'Jam Pattern (mean spacing: {mean_spacing:.2f} km)', color='#6ee7b7', fontsize=13)
    for p in peaks:
        ax2.axvline(x=x[p], color='#f87171', alpha=0.3, linewidth=0.8)
else:
    ax2.set_title('Jam Pattern Formation', color='#6ee7b7', fontsize=13)

ax2.set_xlabel('Position (km)', color='white')
ax2.set_ylabel('Density (veh/km)', color='white')
ax2.legend(facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')

plt.tight_layout(pad=2.0)
plt.savefig('output.png', dpi=150, bbox_inches='tight', facecolor='#0a0a0a')
plt.close()
print("Jam spacing analysis complete.")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Key Takeaways

The LWR equation $\partial_t \rho + \partial_x q = 0$ is the continuum limit of the TASEP mean-field.
Shocks (traffic jams) form when characteristics converge; the Rankine-Hugoniot condition$s = \Delta q / \Delta \rho$ gives the jam propagation speed.
The Lax entropy condition explains why jams have a sharp upstream front and gradual downstream tail.
The Cahn-Hilliard/spinodal decomposition analogy predicts a most unstable wavelength$\lambda^* = 2\pi\sqrt{2\kappa/|D|} \approx 1\text{--}3$ km, matching observed jam spacings.
The Payne-Whitham model adds velocity relaxation, producing negative effective diffusion in the spinodal region.

← TASEP Phase Diagram General ASEP →