Background

The Burgers equation is a nonlinear partial differential equation of second order. It is used in fluid dynamics teaching and in engineering as a simplified model for turbulence, boundary layer behavior, shock wave formation, and mass transport. It has been studied and applied for many decades. Many different closed-form, series approximation, and numerical solutions are known for particular sets of boundary conditions.

Burgers equations is an example of a non-linear scalar equation that produce shocks even starting from smooth initial data. It can be written as:

$$\dot q + q \frac{\partial q}{\partial x} = 0.$$ (30)

The method of solution that we are going to use relays on specific form of the equations, they need to be cast in conservation law form, i.e.:

$$\frac{\partial \,{\bf q}}{\partial t} + \frac{\partial \,{\bf f}({\bf q})}{\partial x^i} = {\mbox{\boldmath$\psi$}},$$ (31)

where $\bf q$ is a vector formed by the so called the conservation variables and $\bf {f}({\bf q})$ is a vector of fluxes which depend on the conservation variables (in general not explicitly). Note that $x^i, \ \ i=1,2,3$ stands for the spatial coordinates. ${\mbox{\boldmath$\psi$}}$ is a source term (when we have no source term we call this equation the inviscid Burgers equation).

Figure 1: A gaussian profile and Burgers equation characteristic lines, shown in the applet we implemented for this homework (see web site)

Figure 2: A shock profile and its characteristic lines.

Figure 3: Now a rare fraction.

It is easy to cast Burgers equation in conservation form:

$$\dot q + \left(\frac{1}{2} q^2 \right)^\prime = 0.$$ (32)

Using the notation introduced previously:

$${\bf q} = q, \ \ \ \ \ \ \ \ {\bf f} = \frac{1}{2} q^2, \ \ \ \ \ \ \ \ {\mbox{\boldmath$\psi$}}=0.$$ (33)

In order to solve this equation we use a finite volumes discretization with a Roe solver. The finite volume discretization of this equation can be written as:

$$\frac{\bar q^{n+1}_{i}-\bar q^n_i}{dt} + \frac{F^{n}_{i+1/2}-F^{n}_{i-1/2}}{dx}=0$$ (34)

where $\bar q^n_i$ is the spatial average [10] and $F^{n}_{i+1/2}$is the numerical flux.

If the initial data is smooth, then the method of characteristics can be used to determine the solution for small enough t such that the characteristics do not intersect. For larger t, after characteristics have interested, the PDE will fail to have a classical solution (a classical solution is a single valued solution or a function), as the information obtained by following the characteristics will produce a multi-valued solution, or possibly, no solution at all. To overcome this lack of existence of a classical solution, we must introduce a broader notion of a solution, a weak solution. Roughly speaking, a weak solution may contain discontinuities, may not be differentiable, and will require less smoothness to be considered a solution than a classical solution. Working with the weak solution of a PDE usually requires that the PDE be reformulated in an integral form. If a classical solution to the problem exists, it will also satisfy the definition of a weak solution.

Figure 4: Godunov's scheme: local solutions of Riemann problems. At every interface,, and, a local Riemann problem is set up as a result of the discretization process (bottom panel), when approximating the numerical solution by piecewise constant data. At time these discontinuities decay into three elementary waves, which propagate the solution forward to the next time level $t^{n+1}$ (top panel). The time step of the numerical scheme must satisfy the Courant-Friedrichs-Lewy condition, being small enough to prevent the waves from advancing more than $\Delta x/2$ in $\Delta t$ (From Numerical Hydrodynamics in General Relativity Jose A. Font, http://www.livingreviews.org/lrr-2003-4)