MATHEMATICS 456/556

Computer Project: The viscous Burgers’ equation and diffuse shocks

The “inviscid” Burgers equation

u_{t} + uu_{x} = 0

can admit shocks, places where the solution develops a discontinuity. One way of dealing with this is to prescribe how the shock moves. Writing the equation as the conservation law

u_{t} + J_{x} = 0, J(u) = \frac{1}{2}u^{2},  (1)

we can insist that, in a frame of reference moving with the shock, the flux J must be the same on either side (otherwise, mass u would pile up at the shock). This leads to Rankine-Hugoniot condition for the shock velocity

V = \frac{J(u_{+}) − J(u_{−})}{u_{+} − u_{−}}

where u_{−}, u_{+} are values of the solution on either side of the shock.

Of course, a discontinuous solution will not solve (1) at the point where the shock occurs; indeed, the derivatives are not even defined there. It has proved advantageous to modify this equation slightly so that shocks do not form, but that solutions are approximations to shock solutions of (1). One way of doing this is to add a diffusive term

u_{t} + uu_{x} = ϵ^{2}u_{xx}.

Without the uu_{x} term, we have the heat equation, which we know smooths things out. The result is a competition between shock formation and diffusion, which gives a “diffusive shock” - not an actual discontinuity, but a zone where the solution transitions rapidly from one value to another.

(1) Download the file “burgers.m”. This is a lot like the explicit algorithm we used for the heat equation. The spatial derivatives are all replaced by their second-order accurate finite-difference formulas. Thus the approximation u_{i,j} ≈ u(iΔx, jΔt) satisfies

\frac{u_{i,j+1} − u_{i,j}}{Δt} + u_{i,j} \frac{u_{i+1,j} − u_{i,j−1}}{2Δx} = ϵ^{2} \frac{−2u_{i,j} + u_{i+1,j} + u_{i,j−1}}{Δx^{2}} .

Note that we still need to pick Δt with care, for the same stability reasons this was necessary for the explicit heat-equation algorithm. The code assumes boundary conditions u(0) = 2 and u(50) = 0, with initial data representing a smooth step from 2 down to 0.

First, test the effect of the parameter ϵ by setting it equal to the values 2, 1, .5. How does the qualitative form of the solution change as ϵ is reduced? What would the limit of solutions look like if ϵ was taken to zero?

(2) Now measure the speed of the moving front (note that initially it’s position is at about x = 10) for different epsilon values. Compare this to the speed given by the Rankine-Hugoniot condition. Is there a trend?

(3) Now change the initial condition to represent a random pulse (uncomment this line in the code) and change the boundary conditions to both be zero. Describe the solution and its asymptotic behavior, and investigate what occurs in changing ϵ.