inviscid Burgers equation in one dimension. The function represents the velocity of a compressible fluid. The equation can model a shock wave: for discontinuous initial data such as
the solution (understood in appropriate weak sense) has the jump from to 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 as a sequence of numbers in the top row of a spreadsheet, I use the difference of adjacent entries as approximation to , multiply it by , and use the result to calculate the solution at the next moment of time. The time step is set at . On the first try, I use the forward difference scheme for , 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
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
Going back to the PDE, rewrite it as
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 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
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
Here is the plot of the data:
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 . Recalling that every row corresponds to time interval of , we can estimate the propagation speed to be about , which is what the analysis of PDE would predict.
