Rarefaction colliding with two shocks

The Burgers equation {v_t+vv_x =0} is a great illustration of shock formation and propagation, with {v(x,t)} representing velocity at time {t} at position {x}. Let’s consider it with piecewise constant initial condition

\displaystyle u(x,0) = \begin{cases} 1,\quad & x<0 \\ 0,\quad & 0<x<1 \\ 2,\quad & 1<x<2 \\ 0,\quad & x>2 \end{cases}

The equation, rewritten as {(v_t,v_x)\cdot (1,v)=0}, states that the function {v} stays constant on each line of the form {x= x_0 + v(x_0,0)t}. Since we know {v(x_0,0)}, we can go ahead and plot these lines (characteristics). The vertical axis is {t}, horizontal {x}.


Clearly, this picture isn’t complete. The gap near {1} is due to the jump of velocity from {0} to {2}: the particles in front move faster than those in back. This creates a rarefaction wave: the gap gets filled by characteristics emanating from {(1,0)}. They have different slopes, and so the velocity {v} also varies within the rarefaction: it is {v(x,t)=(x-1)/t}, namely the velocity with which one has to travel from {(1,0)} to {(x,t)}.


The intersecting families of characteristics indicate a shock wave. Its propagation speed is the average of two values of {v} to the left and to the right. Let’s draw the shock waves in red.


Characteristics terminate when they run into shock. Truncating the constant-speed characteristics clears up the picture:


This is already accurate up to the time {t=1}. But after that we encounter a complication: the shock wave to the right separates velocity fields of varying intensity, due to rarefaction to the left of the shock. Its propagation is now described by the ODE

\displaystyle \frac{dx}{dt} = \frac12 \left( \frac{x-1}{t} + 0 \right) = \frac{x-1}{2t}

which can be solved easily: {x(t) = 2\sqrt{t} +1}.

Similarly, the shock on the left catches up with rarefaction at {t=2}, and its position after that is found from the ODE

\displaystyle \frac{dx}{dt} = \frac12 \left( 1 + \frac{x-1}{t} \right) = \frac{x-1+t}{2t}

Namely, {x(t) = t-\sqrt{2t}+1} for {t>2}. Let’s plot:


The space-time trajectories of both shocks became curved. Although initially, the shock on the right moved twice as fast as the shock on the left, the rarefaction wave pulls them toward each other. What is this leading to? Yet another collision.

The final collision occurs at the time {t=6+4\sqrt{2} \approx 11.5} when the two shock waves meet. At this moment, rarefaction ceases to exist. The single shock wave forms, separating two constant velocity fields (velocities {1} and {0}). It propagates to the right at constant speed {1/2}. Here is the complete picture of the process:


I don’t know where the holes in the spacetime came from; they aren’t predicted by the Burgers equation.

What difference a difference scheme makes

One of the simplest nonlinear PDE is

{\displaystyle u_t = -uu_x}

inviscid Burgers equation in one dimension. The function {u} represents the velocity of a compressible fluid. The equation can model a shock wave: for discontinuous initial data such as

{\displaystyle u(x,0) = \begin{cases} 1, \quad x\le 0 \\ 0, \quad x>0 \end{cases}}

the solution (understood in appropriate weak sense) has the jump from {1} to {0} at all times, but the position of the jump moves to the right. The equation can also demonstrate the creation of shocks from continuous initial data (nonlinear steepening of waves), but I am not going into that here.

Let’s try to solve this PDE numerically. Given the initial {u(x,0)} as a sequence of numbers {\dots 1 1 1 1 1 0 0 0 0 0 \dots} in the top row of a spreadsheet, I use the difference of adjacent entries as approximation to {u_x}, multiply it by {u}, and use the result to calculate the solution at the next moment of time. The time step is set at {\Delta t=0.1}. On the first try, I use the forward difference scheme for {u_x}, which in spreadsheet terms amounts to formulas such as

F2 = F1 - 0.1*F1*(G1-F1)

The result of running this scheme is nothing short of disastrous:

Forward (downwind) difference
Forward (downwind) difference

Let’s try the backward (upwind) difference (“wind direction” is positive, since we work with positive values of velocity). This means using formulas such as

F2 = F1 - 0.1*F1*(F1-E1)

The numerical solution is less ridiculous, but is still quite wrong: it is constant in time.

Backward (downwind) difference
Backward (downwind) difference

Going back to the PDE, rewrite it as

{\displaystyle u_t = -(u^2/2)_x}

This is the conservation form of the PDE, which is more suitable for interpretation of discontinuous solutions. Indeed, it says exactly that the vector field {V(x,t) = (u^2/2,u)} has zero divergence. The property of having zero divergence can be understood without taking derivatives: it means that the flux integral over every closed curve (e.g., over every rectangle) is zero.

Discretizing the conservation form of the PDE leads to finite differences of squares, such as

F2 = F1 - 0.1*(G1^2-F1^2)/2

(conservative downwind method):

Conservative downwind method
Conservative downwind method

Still not good. But the conservative upwind method, with formulas such as

F2 = F1 - 0.1*(F1^2-E1^2)/2

turns out to work:

Conservative upwind method
Conservative upwind method

Here is the plot of the data:

Numerical shock propagation with the conservative upwind method.
Numerical shock propagation with the conservative upwind method.

The shock front looks a bit diffuse, but it sustains itself: note how the value at “4” converges to {1}. Recalling that every row corresponds to time interval of {0.1}, we can estimate the propagation speed to be about {0.5}, which is what the analysis of PDE would predict.