Software 2.5.1 (The breaking time). For Burgers’ equation with \(u_0 = e^{-x^2}\), the breaking time from Theorem 2.1.2, and a direct check that the characteristic map \(x_0 \mapsto x_0 + u_0(x_0)t\) stops being monotone:

>>> import numpy as np
>>> u0 = lambda x: np.exp(-x**2)           # u_t + u u_x = 0
>>> du0 = lambda x: -2 * x * np.exp(-x**2)
>>> xs = np.linspace(-3, 3, 600001)
>>> tb = -1 / du0(xs).min(); tb, 1 / (np.sqrt(2) * np.exp(-0.5))
(1.1658219908227196, 1.165821990798562)
>>> X = lambda x0, t: x0 + u0(x0) * t      # straight characteristics
>>> for t in (1.0, tb, 1.5):
...     print(round(t, 4), bool(np.all(np.diff(X(xs, t)) > 0)))   # ordered?
...
1.0 True
1.1658 False
1.5 False

Reading the output. The minimum of \(u_0'\) on a fine grid gives \(t_b = 1.16582\), agreeing with \(e^{1/2}/\sqrt 2\) to ten digits. At \(t = 1\) the characteristic positions are still in order, so the solution is single-valued; at \(t_b\) and after, they are not, and the classical solution has ceased to exist.

MATLAB/Octave: -1/min(gradient(exp(-x.^2), x)).

Software 2.5.2 (Conservative and non-conservative schemes). Godunov’s method updates cell averages using the exact flux of the Riemann problem at each cell interface. It is compared with a naive upwind discretisation of \(u_t + uu_x = 0\), which is consistent with the equation for smooth solutions but is not in conservation form. Riemann data \(u_l = 2\), \(u_r = 0\) until \(T = 2\); then the reverse data:

>>> def godunov(ul, ur):                   # Riemann flux, f(u) = u^2/2
...     shock = np.where(ul + ur > 0, ul**2, ur**2) / 2
...     fan = np.where(ul > 0, ul**2 / 2, np.where(ur < 0, ur**2 / 2, 0.0))
...     return np.where(ul > ur, shock, fan)
...
>>> def solve(uL, uR, conservative, N=800, L=8.0, T=2.0):
...     dx = L / N; x = -L / 2 + dx * (np.arange(N) + 0.5)
...     u = np.where(x < 0, uL, uR); dt = 0.4 * dx / max(abs(uL), abs(uR))
...     for _ in range(int(round(T / dt))):
...         if conservative:
...             F = godunov(u[:-1], u[1:])
...             u[1:-1] -= dt / dx * (F[1:] - F[:-1])
...         else:                          # upwind on u u_x: not conservative
...             u[1:] -= dt / dx * u[1:] * (u[1:] - u[:-1])
...     return x, u
...
>>> for cons in (True, False):             # uL = 2, uR = 0, until T = 2
...     x, u = solve(2.0, 0.0, cons)
...     print(cons, round(x[np.argmax(u < 1.0)], 3))   # where u drops
...
True 2.005
False 0.005
>>> x, u = solve(0.0, 2.0, True)           # the reverse: a rarefaction fan
>>> [round(float(np.interp(q, x, u)), 3) for q in (-1.0, 1.0, 2.0, 3.0, 5.0)]
[0.0, 0.512, 1.007, 1.498, 2.0]

Reading the output. The exact shock speed is \(\frac 12(2 + 0) = 1\), so at \(T = 2\) the shock is at \(x = 2\). Godunov’s method places it at \(2.005\), within one cell. The non-conservative scheme leaves it at the origin, \(0.005\): it converges, but to the wrong weak solution, because it does not respect the integral form of the law on which the Rankine–Hugoniot condition rests. This is the Lax–Wendroff theorem in negative form: conservative schemes that converge, converge to weak solutions; non-conservative ones need not. For the reverse data, Godunov’s method produces the fan \(u = x/t\): at \(x = 1, 2, 3\) the values \(0.51\), \(1.01\), \(1.50\) against the exact \(0.5\), \(1\), \(1.5\).

MATLAB/Octave: the same vectorised update in a loop.

Software 2.5.3 (A red light). For a road with \(v_{\max } = 100\) km/h and jam density \(150\) cars/km, traffic at \(50\) cars/km meets a queue at a red light:

>>> vmax, rmax = 100.0, 150.0              # km/h, cars per km (jammed)
>>> flux = lambda r: r * vmax * (1 - r / rmax)
>>> speed = lambda r: vmax * (1 - 2 * r / rmax)       # characteristic speed
>>> rL, rR = 50.0, 150.0                   # queue at a red light ahead
>>> (flux(rR) - flux(rL)) / (rR - rL), speed(rL), speed(rR)
(-33.333333333333336, 33.333333333333336, -100.0)
>>> rc = rmax / 2; rc, flux(rc)            # density of maximum flow, capacity
(75.0, 3750.0)

Reading the output. The back of the queue moves backwards at \(33.3\) km/h, equal to \(v_{\max }\rho _l/\rho _{\max }\). The entropy condition holds: the characteristic speeds on the two sides, \(+33.3\) and \(-100\) km/h, both run into the shock. The capacity of the road is \(3750\) cars an hour, at the critical density \(75\) cars/km, where cars travel at \(50\) km/h.

MATLAB/Octave: q = @(r) r*100.*(1 - r/150); (q(150) - q(50))/(150 - 50).