Burgers' equation

In [1]:
%matplotlib inline
In [2]:
%config InlineBackend.figure_format = 'svg'
from utils.snapshot_widgets import interact
from ipywidgets import widgets
from ipywidgets import FloatSlider, fixed
from exact_solvers import burgers
from exact_solvers import burgers_demos
from IPython.display import HTML

Will create static figures with single value of parameters

In this chapter, we study a simple scalar nonlinear conservation law: Burgers' equation. Burgers' equation models momentum transport in a fluid of uniform density and pressure, and it is the simplest equation that captures some key features of gas dynamics or water waves.

To examine the Python code for this chapter, see:

Burgers' equation has been used extensively for developing both theory and numerical methods, and it will allow us to explore the Riemann problem for a nonlinear conservation law. Burgers' equation is a scalar conservation law with flux $f(q)=q^2/2$: \begin{align} q_t + \left(\frac{1}{2}q^2\right)_x = 0. \label{burgers0} \end{align} The quasilinear form is obtained by applying the chain rule to the flux term:
\begin{align*} q_t + qq_x = 0. \end{align*}
This equation looks very similar to the advection equation, with the difference that the advection speed at each point is given by the solution $q$. Burgers' equation is often viewed as a simplified version of equations in fluid dynamics or water waves that also have nonlinear fluxes. Our study of the dynamics of equation (\ref{burgers0}) is a step toward understanding these more complex nonlinear systems, in Shallow_water and Euler. In Traffic_flow we consider another scalar nonlinear conservation law that has a similar structure to Burgers' equation and a simpler physical interpretation.

The characteristic speed for Burgers' equation is $f'(q) = q$. As long as the solution is smooth, the solution is constant along characteristic curves $X(t)$ satisfying $X'(t) = f'(q(X(t),t)$ since then $$ \frac{d}{dt} q(X(t),t) = q_x(X(t),t)X'(t) + q_t(X(t),t) = 0. $$ Because of this, the characteristics are straight lines. However, since the charactistic speed depends on the solution, these lines are not parallel and characterstics may converge or spread out.

Shock formation

In the figure below we consider Burgers' equation with a Gaussian hump as the initial data. Since the characteristic speed in Burgers' equation is given by $q$ itself, the peak of the hump travels faster than the rest, and characteristics are converging at the front of the traveling wave (where $f'(q)$ decreases with $x$) while they are spreading out behind the peak (where $f'(q)$ increases with $x$). The dashed line shows the initial condition while the solid lines show the solution at later times.

In [3]:
%matplotlib inline
In [4]:
%config InlineBackend.figure_format = 'svg'
from utils.snapshot_widgets import interact
from ipywidgets import widgets
from ipywidgets import FloatSlider, fixed
from exact_solvers import burgers
from exact_solvers import burgers_demos
In [5]:
interact(burgers_demos.multivalued_solution, 
         t=FloatSlider(min=0.,max=9.,value=7.75), fig=fixed(0));

Notice that at first $q$ remains single-valued for every $x$. However, after some time the crest of the wave overtakes the leading edge. After this time, we obtain a triple-valued solution for certain values of $x$. The first time this overtaking happens is referred to as the breaking time -- a reference to waves breaking on the beach. It is also the point where the conservation law, in its differential form, breaks down and where the characteristics cross each other for the first time. When characteristics cross, a shock wave, or discontinuity, forms. Mathematically, replacing the triple-valued region with a discontinuity will avoid the problem of the solution being multivalued. Where should the shock be located?

If we replace part of the multivalued solution interval with a shock, some mass will be removed (area $A_1$ in the figure below) and some mass will be added (area $A_2$ in the figure below). In the live notebook, the figure below shows possible locations for the shock at a given time; in order to maintain conservation of the integral of $q$, the shock must be placed so that these two areas are equal.

In [6]:
interact(burgers_demos.shock_location, 
         xshock = FloatSlider(min=6,max=10,step=0.25,value=7.75),
         fig=fixed(0));

This geometric reasoning provides a nice intuition for the shock location, but is cumbersome in practice. To determine the location of the shock we can use the Rankine-Hugoniot jump condition, which we will derive in Traffic_flow, and which requires that the jump in flux across a propagating shock must be related to the jump in $q$ by $$ s(q_r-q_l) = f(q_r) - f(q_l), $$ at each instant in time, where $s$ is the shock speed at this time.

Since the flux for Burgers' equation is $f(q) = q^2/2$, this gives \begin{align*} s(q_r-q_l) & = \frac{1}{2} (q_r^2 - q_l^2) \\ \Rightarrow \ \ \ \ s &= \frac{1}{2}(q_l + q_r). \end{align*}
In general (as in the image above) the states $q_l$ and $q_r$ just to the left and right of the shock will not be constant and the speed of a shock will change in time.

Shock solution

For Burgers' equation, (or any scalar hyperbolic conservation law with a convex flux function, such as the traffic flow problem considered in Traffic_flow), the solution of the Riemann problem consists of a single wave, which may be either a shock or a rarefaction, separating regions of constant value $q_l$ and $q_r$. Here convexity of the flux $f(q)$ means that $f'(q)$ is either monotonically increasing or monotonically decreasing as $q$ varies. If $f(q)$ is not convex then the solution can be more complicated; see the section on Convexity below and in Nonconvex_scalar.

The analysis above already yields the solution to the Riemann problem in the case that the resulting wave is a shock. The entropy condition, as we will explain below, indicates that a shock will form when $f'(q_l) > f'(q_r)$; i.e. when $q_l>q_r$. In this case, the solution is \begin{align*} q(x,t) = \begin{cases} q_l \quad \text{if} \ \ x/t<s \\ q_r \quad \text{if} \ \ x/t>s. \end{cases} \end{align*}
Below, we plot the solution of Burgers' equation for an initial condition that leads to a shock (i.e. with $q_l>q_r$).

In [7]:
interact(burgers_demos.shock(), 
         t=widgets.FloatSlider(min=0,max=1.0,value=0.5),
         which_char=widgets.Checkbox(value=True,
                                     description='Show characteristics'));

Rarefaction wave

In the previous figure with the hump as the initial condition, we observed that a shock formed on the right side of the hump. However, on the left side, the characteristics spread out and will never cross. This part of the solution is called a rarefaction wave. This is the kind of behavior we will observe in the solution of the Riemann problem when $q_l<q_r$.

In the next figure, we consider such a Riemann problem. Although the initial data is discontinuous at $x=0$, we can think of smearing it out slightly to a continuous function so that all values between $q_l$ and $q_r$ are taken over a very narrow region around $x=0$. As time evolves, each value of $q$ propagates according to the advection speed $q$ given by the quasi-linear equation. Therefore, after time $t$, each value $q$ must propagate a distance $x=qt$, so the solution along the rarefaction is then $q = x/t$. As $q_l<q_r$, the smallest and largest displacements are given by $q_l t$ and $q_r t$, respectively.

In [8]:
interact(burgers_demos.rarefaction_figure, 
         t = FloatSlider(min=0,max=9,value=5));

The full rarefaction solution for the Burgers' Riemann problem is then simply given by \begin{align*} q(x,t) = \begin{cases} q_l, \quad \text{for} \ \ x<q_l t \\ x/t, \quad \text{for} \ \ q_l t \le x \le q_r t \\ q_r, \quad \text{for} \ \ x>q_r t. \end{cases} \end{align*}

As we will see in Traffic_flow, the rarefaction solution is always a self-similar solution. This means that it can be expressed as a function of the ratio between position and time $q(x,t) = \tilde{q}(x/t)$, so it remains the same when rescaling both $x$ and $t$ by the same factor. In Burgers' equation, the form of the rarefaction is particularly simple since the advection speed is simply $q$ and so the rarefaction wave is linear in $x$ with slope $1/t$ at time $t$.

Rarefaction solution

In the figure below, we plot a solution of the Riemann problem with $q_l<q_r$, exhibiting a rarefaction.

In [9]:
interact(burgers_demos.rarefaction(), 
         t=widgets.FloatSlider(min=0,max=1.0,value=0.5),
         which_char=widgets.Checkbox(value=True,
                                     description='Show characteristics'));

Weak solutions

As we mentioned before, the differential form of the equation breaks down in the presence of shocks/discontinuities. However, the integral form of the conservation law remains valid. Let's integrate the general conservation law $q_t+f(q)_x=0$ from $x=x_1$ to $x=x_2$ and $t=t_1$ to $t=t_2$: \begin{align*} \int_{t_1}^{t_2}\int_{x_1}^{x_2} [q_t+f(q)_x] dx dt = 0. \end{align*} This integral can be rewritten in terms of an indicator function $\phi(x,t)$ in $[x_1,x_2]\times[t_1,t_2]$ (defined to be 1 in this region, zero elsewhere): \begin{align} \int_{0}^{\infty}\int_{-\infty}^{\infty} [q_t+f(q)_x]\phi(x,t) dx dt = 0. \label{eq:Burgersintclaw2} \end{align}

We can further replace $\phi(x,t)$ by a smooth function with compact support on some region of the $x-t$ plane (zero outside a closed and bounded region). Assuming $t=0$ is in the support of $\phi(x,t)$, integration by parts yields \begin{align} \int_{0}^{\infty}\int_{-\infty}^{\infty} [q\phi_t+f(q)\phi_x] dx dt = -\int_{-\infty}^{\infty}q(x,0)\phi(x,0)dx, \label{eq:Burgersintclaw3} \end{align} where now the derivatives are on $\phi(x,t)$ and not on $q$ or $f(q)$, so the equation still makes sense for discontinuous $q$. Note that we only obtain one boundary term along $t=0$ since $\phi(x,t)$ vanishes at infinity. The function $q(x,t)$ is called a weak solution of the conservation law if (\ref{eq:Burgersintclaw3}) holds for all functions $\phi(x,t)$ that are continuously differentiable and have compact support (bump functions). However, the function $\phi(x,t)$ chosen in (\ref{eq:Burgersintclaw2}) does not satisfy these conditions since it is not smooth. Nonetheless, we can approximate this function arbitrarily well by a smooth function. Note that any weak solution is a solution of the integral conservation law and vice versa.

Weak solutions are thus allowed to include discontinuities. The shock solution presented above, for instance, is a weak solution of Burgers' equation. After characteristics cross, there is no strong solution of the conservation law, and we must resort to weak solutions.

An unphysical weak solution

A potential problem with the notion of weak solutions is that in some cases the weak solution is not unique. For example, consider the Riemann problem for Burgers' equation with $q_l < q_r$. As mentioned already we expect a rarefaction to be the correct solution. This is because if we smeared out the initial data just slightly then there would be smoothly varying characteristics emerging from each point $x$ and these characteristics would spread out. However, for exactly discontinuous initial data, there exists another weak solution in which the initial discontinuity propagates as a shock wave with the Rankine-Hugoniot speed $(q_l + q_r)/2$: \begin{align*} q(x,t) = \begin{cases} q_l \quad \text{if} \ \ x/t<s \\ q_r \quad \text{if} \ \ x/t>s. \end{cases} \end{align*} This unphysical solution is also referred to as an expansion shock, and it is also a weak solution to Burgers' equation. In the next figure, we plot this weak solution.

In [10]:
interact(burgers_demos.unphysical(), 
         t=widgets.FloatSlider(min=0,max=1.0,value=0.5),
         which_char=widgets.Checkbox(value=True,
                                     description='Show characteristics'));

Note the behavior of the characteristics with respect to the shock. The fact that the characteristics are spreading away from the shock rather than converging on it indicates that the correct solution should instead be a rarefaction. In order to be able to specify which of the weak solutions is physically correct, we need to derive a mathematical condition from our physical intuition gained from observing the behavior of the characteristics. This condition is referred to as the entropy condition.

The entropy condition

A given initial value problem for a hyperbolic PDE may have many weak solutions. A condition that selects a unique physically correct solution out of these weak solutions is called an admissibility condition or more often an entropy condition. This name comes from gas dynamics, where the physical entropy must increase in the gas passing through a shock, according to the second law of thermodynamics. A discontinuity that violates this is non-physical. A mathematical "entropy function" with a similar property can often be defined for other conservation laws. More discussion of this and several other formulations of admissibility conditions, with more detailed explanations, are available in the literature, for example in many of the books cited in the Preface.

In the context of the Riemann problem, the entropy condition allows us to determine if the physical solution should involve a shock or a rarefaction. In the case of scalar conservation laws, the entropy condition is quite simple and can be formulated in terms of the flux function of the conservation law as follows: the solution of a scalar Riemann problem will consist of a shock only if $$ f'(q_l) > f'(q_r). $$ In other words, the solution is a shock only if nearby characteristics from the left and right approach the shock as time progresses. This is often called the Lax Entropy Condition. If nearby characteristics are spreading out, as in the last example, the correct solution is instead a rarefaction. In the case of Burgers' equation, the flux function is $f(q)=q^2/2$, so the correct solution is a shock only if $$ q_l > q_r, $$ which can be clearly observed in the interactive solution below (in the notebook). In later chapters, we will see how this condition generalizes to systems of conservation laws.

Convexity

The flux $f(q) = q^2/2$ for Burgers' equation is a convex function, since $f''(q) = 1 > 0$ for all values of $q$. This means that as $q$ varies between $q_l$ and $q_r$ the the characteristic speed $f'(q) = q$ is either monotonically increasing (if $q_l < q_r$) or monotonically decreasing (if $q_l > q_r$). Hence for any Riemann problem the characteristic speeds for intermediate states are either purely converging or purely diverging, and the Riemann solution is always either a single rarefaction wave or a single shock wave.

The flux $f(\rho) = \rho(1-\rho)$ of the simple traffic flow model that we will consider in Traffic_flow also has the property that $f'(\rho)$ varies monotonically with $\rho$ (in the context of hyperbolic PDEs, a flux function is referred to as convex if either $f''(q)\ge0$ for all $q$ or $f''(q)\le 0$ for all $q$). Since $f''$ is always negative for the traffic flow model, the correct Riemann solution is a shock if $\rho_l < \rho_r$, or a rarefaction wave if $\rho_l > \rho_r$.

The solution to the Riemann problem can be much more complicated if $f'(q)$ is not monotonically varying between $q_l$ and $q_r$, i.e. if $f''(q)$ changes sign. In this case the Riemann solution can consist of multiple shock and rarefaction waves. The nonconvex case is explored further in Nonconvex_scalar.

Interactive solution and examples

In the live notebook, the figure below is an interactive solution of the Riemann problem for Burgers' equation. The values of the initial conditions and the time can be modified to observe their effect on the characteristic structure and the solution.

In [11]:
burgers.riemann_solution_interact()

Animations of more complex solutions

We can now explore a few more examples that are representative of phenomena we will observe in more complicated systems, with more interesting initial data that the single jump discontinuity of a Riemann problem. For the animations shown below, the solution is computed numerically using PyClaw.

Bump initial conditions

The animation below in the notebook, or on this webpage, shows the solution of Burgers' equation for an initial Gaussian hump, as discussed at the beginning of this chapter. A shock forms on the leading edge and the trailing edge spreads out as a rarefaction.

In [12]:
anim = burgers_demos.bump_animation(numframes = 50)
HTML(anim)
Out[12]:


Once Loop Reflect

Three-state data with merging shocks

For an initial condition that is piecewise-constant with three values, one can solve Burgers' equation by solving a Riemann problem locally for each discontinuity and then considering how the resulting waves interact. For instance, consider the initial condition \begin{align*} q_0(x) = \begin{cases} q_l \quad \text{if} \ \ x < -1 \\ q_m \quad \text{if} \ \ -1\le x \le 1 \\ q_r \quad \text{if} \ \ x > 1. \end{cases} \end{align*}
We can decompose this into two Riemann problems: one at $x=-1$ and another at $x=1$. Each of these two Riemann problems will yield either a shock or a rarefaction. The interesting part is what happens when the generated shocks or rarefactions interact. In the scenario shown below, let us assume that both Riemann problems produce a right-propagating shock. However, the shock generated at $x=-1$ propagates faster, so it will eventually reach the shock originally generated at $x=1$. At the point in time when the shocks collide, one can simply restate the problem again as a Riemann problem and solve the whole problem analytically. In this animation, the plot on the left shows the solution $q(x)$ evolving through time, while the one on the right shows the characteristic structure in the $x-t$ plane, with the shocks shown as wide lines and the characteristics as thin lines, and with time is marked by a horizontal dashed line. In this animation, one can clearly see the two shocks merge to become a single shock at later times. This animation can be viewed in the live notebook on or this webpage.

In [13]:
anim = burgers_demos.triplestate_animation(ql = 4.0, qm = 2.0, 
                                           qr = 0.0, numframes = 50)
HTML(anim)
Out[13]:


Once Loop Reflect

Three-state data with interacting shock and rarefaction

As one would expect, this can become more complicated when one or both of the waves generated are given by rarefactions. The animation below, or on this webpage, shows how a right-propagating shock can overtake a rarefaction. The shock speed changes as it passes through the rarefaction, since the state just to the left of the shock is changing. This is seen clearly in the $x-t$ plane, where the slope of the shock is changed after interacting with the rarefaction. Note the shock is again marked as a wider line.

In [14]:
anim = burgers_demos.triplestate_animation(ql = 4.0, qm = -1.5, 
                                           qr = 0.5, numframes = 50)
HTML(anim)
Out[14]:


Once Loop Reflect

Analogously, one can also have a rarefaction overtaking a shock. As seen in the animation below, or this webpage, this can even change the direction in which the shock is propagating. This happens because the velocity of the shock depends on the slopes of the impinging characteristics. When a shock interacts with a rarefaction, the varying slopes of the characteristic curves produce a constant change of the shock propagation velocity.

In the notebook, try modifying the values of $q_l$, $q_m$ and $q_r$ in order to see what other behavior can be observed. Is it possible to observe two interacting rarefactions?

In [15]:
anim = burgers_demos.triplestate_animation(ql = -1, qm = 3.0, 
                                           qr = -2, numframes = 50)
HTML(anim)
Out[15]:


Once Loop Reflect