# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# + [markdown] slideshow={"slide_type": "skip"}
# <table>
#  <tr align=left><td><img align=left src="./images/CC-BY.png">
#  <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) <NAME></td>
# </table>

# + slideshow={"slide_type": "skip"}
from __future__ import print_function

# %matplotlib inline

import os
import numpy
import matplotlib.pyplot as plt
import matplotlib.animation

from IPython.display import HTML
from clawpack import pyclaw
from clawpack import riemann


# + [markdown] slideshow={"slide_type": "slide"}
# # Nonlinear Scalar Conservation Laws

# + [markdown] slideshow={"slide_type": "subslide"}
# Up until now we have mostly assumed that our PDE
# $$
#     q_t + f(q)_x = 0
# $$
# has been linear, i.e. $f(q) = A q$.  Part of this was to establish the basic theory while we know what the true solution should look like.  We now will investigate how nonlinear problems do not simply translate but can deform the solution in interesting ways.  Note that we will assume that $f(q)$ is convex/concave, i.e. $f''(q)$ does not change sign.

# + [markdown] slideshow={"slide_type": "subslide"}
# As a test problem we will consider basic traffic flow as our PDE of choice.  For this we assume that we can model traffic, or the density of cars, as a continuum as we do for fluids.  This leads us to the conservation law
# $$
#     q_t + (U(q) q)_x = 0
# $$
# where $q$ represents the density of cars and $U(q)$ the speed of those cars.  We then need to specify the relationship between the speed at which the cars travel and the local density of cars.

# + [markdown] slideshow={"slide_type": "subslide"}
# The Lighthill, Whitham, and Richards (LWR) equation uses a simple function for the speed that is
# $$
#     U(q) = u_{\text{max}}(1 - q) \text{   for } 0 \leq q \leq 1
# $$
# where $u_{\text{max}}$ is the speed limit and assumed constant.  Adding this to our previous equation leads us to a flux defined as
# $$
#     f(q) = u_{\text{max}}(1 - q) q
# $$
# that is concave with respect to $q$.

# + slideshow={"slide_type": "skip"}
def traffic_animation(init_condition="gaussian"):
    solver = pyclaw.ClawSolver1D(riemann.traffic_1D)
    solver.bc_lower[0] = pyclaw.BC.extrap
    solver.bc_upper[0] = pyclaw.BC.extrap

    x = pyclaw.Dimension(-30.0,30.0,500,name='x')
    domain = pyclaw.Domain(x)
    num_eqn = 1
    state = pyclaw.State(domain,num_eqn)

    grid = state.grid
    xc=grid.p_centers[0]

    if init_condition.lower() == "shock":
        state.q[0, :] = 0.75 * (xc > 0)
    elif init_condition.lower() == "rarefaction":
        state.q[0, :] = 0.75 * (xc < 0)
    elif init_condition.lower() == "gaussian":
        state.q[0, :] = 1.0 * numpy.exp(-xc**2 / 10.0**2)
    else:
        raise ValueError("Unknown initial condition requested.")

    state.problem_data['efix']=True
    state.problem_data['umax']=1.

    claw = pyclaw.Controller()
    claw.tfinal = 25.0
    claw.num_output_times = 20
    claw.solution = pyclaw.Solution(state,domain)
    claw.solver = solver
    claw.keep_copy = True
    
    claw.run()
    x = claw.frames[0].grid.dimensions[0].centers
    
    fig = plt.figure()
    axes = plt.subplot(1, 1, 1)
    axes.set_xlim((x[0], x[-1]))
    axes.set_ylim((-0.1, 1.1))
    axes.set_title("Traffic Flow")
    
    def init():
        axes.set_xlim((x[0], x[-1]))
        axes.set_ylim((-0.1, 1.1))
        computed_line, = axes.plot(x[0], claw.frames[0].q[0, :][0], 'bo-')
        return (computed_line, )
    
    computed_line, = init()
    
    def fplot(n):
        computed_line.set_data([x,], [claw.frames[n].q[0, :]])
        return (computed_line, )

    frames_to_plot = range(0, len(claw.frames))
    plt.close(fig)
    return matplotlib.animation.FuncAnimation(fig, fplot, frames=frames_to_plot, interval=100,
                                   blit=True, init_func=init, repeat=False)
HTML(traffic_animation().to_jshtml())

# + [markdown] slideshow={"slide_type": "subslide"}
# Let us first examine how the characteristics behave in this simple case where the initial condition is smooth but leads to a discontinuous solution.  If we think of the $k$th car whose trajectory is $X_k(t)$.  If we assume that the $k$th car makes a decision on how fast she goes $U_(q_k(t))$ then we have
# $$
#     X'_k(t) = U_k(q_k(t)) = U([X_{k+1}(t) - X_k(t)]^{-1}) \quad \text{for} ~ k=1,\ldots,m
# $$
# with the boundaries appropriately specified.

# + [markdown] slideshow={"slide_type": "subslide"}
# Note that for these test cases that the characteristic speed agrees with the speed at which the particles are moving.  This is a general property of nonlinear hyperbolic PDEs and arises as $f'(q) \neq u$ as it is in the linear case.

# + [markdown] slideshow={"slide_type": "slide"}
# ## Quasilinear Forms and Characteristics
#
# One useful way we can analyze what may go on given a nonlinear flux is to study the quasilinear form of the PDE defined as
# $$
#     q_t + f'(q)q_x = 0.
# $$
# If $q(x,t)$ is smooth then we would have
# $$
#     X'(t) = f'(q(X(t), t))
# $$
# and
# $$
#     \frac{\text{d}}{\text{d}t} q(X(t), t) = X'(t) q_x + q_t = 0
# $$
# and therefore $q(x,t)$ is constant on the characteristic curves defined by $X(t)$.

# + [markdown] slideshow={"slide_type": "subslide"}
# As we know once we have the characteristics and know that $q$ needs to be constant on them that we need to simply trace back the characteristics to the initial condition (or boundary condition) to find out what the solution is.  Often times we represent this as
# $$
#     q(x,t) = q_0(\xi)
# $$
# where
# $$
#     x = \xi + f'(q_0(\xi)) t.
# $$

# + [markdown] slideshow={"slide_type": "slide"}
# ## Burgers Equation
#
# Another of the often used basic equations that are solved when considering nonlinear conservation laws is Burgers equation, defined as
# $$
#     u_t + \left( \frac{1}{2} u^2 \right)_x = 0.
# $$
# Sometimes this is also referred to as the inviscid Burgers equation as the Burgers equation is meant to be a simplified model of fluid flow which usually also includes viscosity in
# $$
#     u_t + \left( \frac{1}{2} u^2 \right)_x = \epsilon u_{xx}.
# $$
# We will use this eventually to define the set of solutions that are "physical".

# + [markdown] slideshow={"slide_type": "subslide"}
# The Burgers equation also can be written in non-conservative form as
# $$
#     u_t + u u_x = 0
# $$
# but as we will see for discontinuous solutions that we will come to different answers than the conservative form.  We will discuss this later again in relation to entropy conditions.

# + slideshow={"slide_type": "skip"}
path = os.path.join(os.environ.get("CLAW", os.getcwd()), "pyclaw", "fvmbook", "chap11", "burgers")
os.chdir(path)
import burgers1D

def burgers_animation():
   
    # compute the solution with the method define above:
    claw = burgers1D.burgers()
    claw.keep_copy = True
    claw.run()
    x = claw.frames[0].grid.dimensions[0].centers
    
    fig = plt.figure()
    axes = plt.subplot(1, 1, 1)
    axes.set_xlim((x[0], x[-1]))
    axes.set_ylim((-3.5, 6))
    axes.set_title("Burgers Equation")
    
    def init():
        axes.set_xlim((x[0], x[-1]))
        axes.set_ylim((-3.5, 6))
        computed_line, = axes.plot(x[0], claw.frames[0].q[0, :][0], 'bo-')
        return (computed_line, )
    
    computed_line, = init()
    
    def fplot(n):
        computed_line.set_data([x,], [claw.frames[n].q[0, :]])
        return (computed_line, )

    frames_to_plot = range(0, len(claw.frames))
    plt.close(fig)
    return matplotlib.animation.FuncAnimation(fig, fplot, frames=frames_to_plot, interval=100,
                                   blit=True, init_func=init, repeat=False)
HTML(burgers_animation().to_jshtml())

# + [markdown] slideshow={"slide_type": "slide"}
# ## Rarefaction Waves
#
# The first type of nonlinaer wave we will consider is ararefaction wave.  In traffic flow this occurs if $q_x(x,0) < 0$.  The characteristic speed is
# $$
#     f'(q) = U(q) + q U'(q) = u_{\text{max}} ( 1- 2q).
# $$

# + slideshow={"slide_type": "skip"}
HTML(traffic_animation(init_condition="rarefaction").to_jshtml())

# + [markdown] slideshow={"slide_type": "subslide"}
# From this example we see that the cars are being "rarefied", spreading out, or becoming less dense and that the solution is smooth.  In this case as well we sometimes call this a **centered rarefaction wave** as it is centered at zero and moving in both directions.

# + [markdown] slideshow={"slide_type": "slide"}
# ## Compression Waves
#
# The original solution of traffic flow is an example of a **compression wave**.  These types of waves usually are composed of a shock and then a rarefaction wave and are called compression waves as they represent a compression of density of cars.

# + slideshow={"slide_type": "skip"}
HTML(traffic_animation(init_condition="gaussian").to_jshtml())

# + [markdown] slideshow={"slide_type": "slide"}
# ## Vanishing Viscosity Limit
#
# As we know, once a shock form the strong form of the PDE breaks down and only the integral form, or **weak form**, is valid.
#
# One way to avoid this situation is to add a vanishingly small amount of viscosity to the solution such as
# $$
#     q_t + f(q)_x = \epsilon q_{xx}.
# $$
# This keeps the shock from completely forming and is also motivated by fluid mechanics where often times the viscosity is small with respect to the other terms in the equation.  

# + [markdown] slideshow={"slide_type": "subslide"}
# Unfortunately this formulation is now a parabolic equation so often we look at this as a means for picking a limiting solution as we will see later.  As $\epsilon \rightarrow 0$ we recover the **vanishing viscosity limit** and re-obtain the original hyperbolic formulation.

# + [markdown] slideshow={"slide_type": "slide"}
# ## Equal-Area Rule
#
# The equal-area rule concerns how to take a solution that is no longer a function, similar to a breaking wave, and placing a shock such that the triple-valued function divides the overlap into equal areas under the curve.  This in reality is a consequence of conservation, we do not want to loose the value of the integral of $q(x,t)$ and so place the shock accordingly.

# + [markdown] slideshow={"slide_type": "slide"}
# ## Shock Speed
#
# Although the equal-area rule works we often desire a m ore fundamental way to know the location and speed of a shock.  Say that we have a shock traveling at a speed $s(t)$.  The integral conservation law provides us a way to find this speed by integrating over the rectangle $[x_0, x_0 + \Delta x] \times [t_0, t_0 + \Delta t]$ that the shock will exactly split from corner to corner.  If we set up this integral we then have
# $$
#     \int^{x_0 + \Delta x}_{x_0} q(x, t_0 + \Delta t) dx - \int^{x_0 + \Delta x}_{x_0} q(x, t_0) dx = \int^{t_0 + \Delta t}_{t_0} f(q(x_0, t)) dt - \int^{t_0 + \Delta t}_{t_0} f(q(x_0 + \Delta x, t)) dt.
# $$

# + [markdown] slideshow={"slide_type": "subslide"}
# By design the values along the rectangle are constant and the integral conservation law can be rewritten as
# $$\begin{aligned}
#     \int^{x_0 + \Delta x}_{x_0} q(x, t_0 + \Delta t) dx - \int^{x_0 + \Delta x}_{x_0} q(x, t_0) dx &= \int^{t_0 + \Delta t}_{t_0} f(q(x_0, t)) dt - \int^{t_0 + \Delta t}_{t_0} f(q(x_0 + \Delta x, t)) dt \\
#     \Delta x q_r - \Delta x q_\ell &= \Delta t f(q_\ell) - \Delta t f(q_r) + \mathcal{O}(\Delta t^2).
# \end{aligned}$$
# Noting that the speed of the shock should approximately be
# $$
#     s(t) = \frac{\Delta x}{\Delta t} + \mathcal{O}(\Delta t^2)
# $$
# as this rectangle shrinks we have
# $$
#     s(t) (q_r - q_\ell) = f(q_r) - f(q_\ell).
# $$
# This expression is often called the **Rankine-Hugoniot jump condition**.

# + [markdown] slideshow={"slide_type": "subslide"}
# If we have a scalar conservation law then this condition is often rewritten as
# $$
#     s = \frac{f(q_r) - f(q_\ell)}{q_r - q_\ell}.
# $$

# + [markdown] slideshow={"slide_type": "subslide"}
# Note that we have written the jump condition in terms of $q_r$ and $q_\ell$.  In general these values are changing in time and could cause the shock speed to change in time.

# + [markdown] slideshow={"slide_type": "subslide"}
# In the case where $q_r \approx q_\ell$ we see that this is starting to approach a derivative of $f(q)$ with respect to $q$ and the problem starts to become linear.  We often call this type of shock a **weak shock**.

# + [markdown] slideshow={"slide_type": "subslide"}
# For traffic flow we have
# $$
#     s = u_{\text{max}} [1 - (q_\ell + q_r)] = \frac{1}{2} [ f'(q_\ell) + f'(q_r)].
# $$ and for Burgers equation we have
# $$
#     s = \frac{1}{2} (u_\ell + u_r).
# $$
#
# This averaging is not coincidental.  For quadratic flux functions the shock speed is the average of the characteristic speeds on either side of the shock.  This is also indicative of the equal-area rule.

# + [markdown] slideshow={"slide_type": "slide"}
# ## Rankine-Hugoniot Conditions for Systems
#
# In the case of systems we can of course not divide through by the difference in states to isolate the shock speed.  Instead we must in general solve a nonlinear system to find the shock speed.

# + [markdown] slideshow={"slide_type": "subslide"}
# In the case of a linear hyperbolic PDE we can start to see though how we may accomplish this.  If $f(q) = A q$ then we would have
# $$
#     A (q_r - q_\ell) = s(q_r - q_\ell).
# $$
# This is of course a eigenproblem where the shock speed $s$ corresponds to the eigenvalues and the jump in $q$ the eigenvectors.

# + [markdown] slideshow={"slide_type": "slide"}
# ## Similarity Solutions and Centered Rarefactions
#
# So far we have not specifically cared about Riemann problem data.  In this case a conservation law in general will be a similarity solution.  In the case of our form of equation then the solution is a function
# $$
#     q(x,t) = \widetilde{q}(\xi) \quad \text{where} \quad \xi = x/t.
# $$

# + [markdown] slideshow={"slide_type": "subslide"}
# Using this fact in combination with our Riemann problem initial condition we can compute
# $$
#     q_t(x,t) = -\frac{x}{t^2} \widetilde{q}'(\xi) \quad \text{and} \quad f(q)_x = \frac{1}{t} f'(\widetilde{q}(\xi)) \widetilde{q}'(\xi)
# $$
# leading to
# $$\begin{aligned}
#     f(q)_x &= -q_t \\
#     \frac{1}{t} f'(\widetilde{q}(\xi)) \widetilde{q}'(\xi) &= \frac{x}{t^2} \widetilde{q}'(\xi) \\
#     f'(\widetilde{q}(\xi)) \widetilde{q}'(\xi) &= \xi \widetilde{q}'(\xi)
# \end{aligned}$$
# For the scalar case either $\widetilde{q}(\xi) = 0$ or 
# $$
#     f'(\widetilde{q}(\xi)) = \xi
# $$
# allowing us to calculate the solution of the centered rarefaction directly.

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Example:  Traffic Flow
#
# We know that if $q_\ell > q_r$ that we have a rarefaction.  This implies that
# $$
#     \widetilde{q}(\xi) = \left \{ \begin{aligned}
#         q_\ell && \text{for } \xi \leq f'(q_\ell) \\
#         q_r && \text{for } \xi \geq f'(q_r) 
# \end{aligned} \right .$$
# This leaves a wedge, the rarefaction wave to still compute.  From our statement before we know
# $$\begin{aligned}
#     f'(\widetilde{q}(\xi)) &= \xi \\
#     u_\text{max} [ 1 - 2 \widetilde{q}(\xi)] &= \xi \\
#     \widetilde{q}(\xi) &= \frac{1}{2} \left[1 - \frac{\xi}{u_\text{max}} \right] \quad \text{for} \quad f'(q_\ell) \leq \xi \leq f'(q_r)
# \end{aligned}$$
# Note that this satisfies the condition that the solution is continuous.

# + [markdown] slideshow={"slide_type": "slide"}
# ## Weak Solutions
#
# We will now derive a different from of the conservation law, more similar to the weak form of PDEs that are common in finite element analysis.  If we instead start from the strong form and integrate over space-time arbitrarily we have
# $$
#     \int^{t_1}_{t_0} \int^{x_1}_{x_0} [q_t + f(q)_x] dx dt = 0.
# $$
# Introduce a **test function** $\phi(x,t)$ and rewrite the above integral as
# $$
#     \int^{\infty}_{0} \int^{\infty}_{-\infty} [q_t + f(q)_x] \phi(x, t) dx dt = 0.
# $$
# where we require that $\phi(x,t)$ has compact support.  In particular we would get back our original integral if
# $$
#     \phi(x,t) = \left \{ \begin{aligned}
#         1 &\quad\text{if} ~ (x, t) \in [x_0, x_1] \times [t_0, t_1] \\
#         0 &\quad \text{otherwise}
#     \end{aligned} \right .
# $$

# + [markdown] slideshow={"slide_type": "subslide"}
# Integrating by parts then gives us
# $$
#     \int^{\infty}_{0} \int^{\infty}_{-\infty} [q \phi_t + f(q) \phi_x] dx dt = - \int^\infty_0 q(x, 0) \phi(x, 0) dx.
# $$
# Note that now the derivatives are no longer on $q$ or $f(q)$ but on the test functions.

# + [markdown] slideshow={"slide_type": "subslide"}
# The function $q(x,t)$ is called a **weak solution** of the conservation law with given initial data $q(x,0)$ if
# $$
#     \int^{\infty}_{0} \int^{\infty}_{-\infty} [q \phi_t + f(q) \phi_x] dx dt = - \int^\infty_0 q(x, 0) \phi(x, 0) dx
# $$
# holds for all function $\phi \in C^1_0$.

# + [markdown] slideshow={"slide_type": "slide"}
# ## Why Should You Be Careful When Manipulating Conservation Laws?
#
# Transforming a conservation law into different differential forms has a consequence that the weak solutions one might find may not be equivalent.  This is best shown by an example.

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Example:  Burgers' Equation
#
# The canonical form of Burgers' equation is
# $$
#     u_t + \left( \frac{1}{2} u^2 \right )_x = 0
# $$
# with a quasilinear form defined by
# $$
#     u_t + u u_x = 0.
# $$
#
# If we multiply this latter form by $2u$ we get
# $$\begin{aligned}
#     2 u u_t + 2 u^2 u_x &= 0 \\
#     (u^2)_t + \left(\frac{2}{3} u^3 \right)_x &= 0.
# \end{aligned}$$
# This form again is now a conservation law for $u^2$.  For smooth solutions, i.e. strong solutions, the first form and last form have equivalent solutions.  However for weak solutions we find differences.  If we have shocks in our solution, i.e. $u_\ell > u_r$, then we have shock speeds
# $$
#     s_1 = \frac{1}{2} (u_\ell + u_r)
# $$
# as before and
# $$
#     s_2 = \frac{2}{3} \left(\frac{u_r^3 - u_\ell^3}{u^2_r - u_\ell^2} \right).
# $$
# Clearly $s_1 \neq s_2$ in this case.  What we have done requires that the solution $u$ must be smooth and therefore leads to a different set of weak solutions.

# + [markdown] slideshow={"slide_type": "slide"}
# ## Nonuniqueness, Admissibility , and Entropy Conditions
#
# Unfortunately when we are computing weak solutions to a PDE we often find that nonuniqueness can be an issue.  In this last section we consider ways to pick out the relevant solution.

# + [markdown] slideshow={"slide_type": "subslide"}
# **Basic Problem:** is that if we consider the weak form of the equation for a Riemann problem (or others for that matter) that both a rarefaction and a shock are solutions.  We know that there really is only one solution but the weak form allows both so we need a means of picking the solution that is correct.

# + [markdown] slideshow={"slide_type": "subslide"}
# One idea we have already mentioned, the vanishing viscosity limit, can be used to pick out a solution.  We know that in the case where we have added the elliptic term $\epsilon u_{xx}$ that the equation now has a unique, smooth solution.

# + [markdown] slideshow={"slide_type": "subslide"}
# Often these types of conditions are called **admissibility conditions** or **entropy conditions**.  The latter often includes a function $\eta(q)$ called an **entropy function** that helps to determine the solution that we want.

# + [markdown] slideshow={"slide_type": "subslide"}
# The first entropy condition we will define is used only for scalar conservation laws but is useful none the less.
#
# **Entropy Condition - Lax:**  For a convex scalar conservation law, a discontinuity propagating with speed $s$ given by the Rankine-Hugoniot condition satisfies the Lax entropy condition if
# $$
#     f'(q_\ell) > s > f'(q_r).
# $$
#
# In this condition $f'(q)$ is the characteristic speed so this condition suggests that the characteristics on either side should be *impinging* on the shock.

# + [markdown] slideshow={"slide_type": "subslide"}
# If the flux function is non-convex we do have a more general entropy condition as well.
#
# **Entropy Condition - Oleinik:**  $q(x,t)$ is the entropy solution to a scalar conservation law $q_t + f(q)_x = 0$ with $f''(q) > 0$ if there exists a constant $E > 0$ s.t. $\forall a > 0, t > 0$, and $x \in \mathbb R$,
# $$
#     \frac{q(x + a, t) - q(x,t)}{a} < \frac{E}{t}.
# $$

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Entropy Functions
#
# As mentioned an alternative approach to admissibility conditions is defining an entropy function.  Usually these are derived from some physical principle as the namesake suggests.
#
# In general an entropy should be conserved when $q(x,t)$ is smooth and decrease or increase at a shock.

# + [markdown] slideshow={"slide_type": "subslide"}
# One example is the entropy of gas.  By the laws of thermodynamics we know that entropy must be produced at an admissible shock but is reduced across an inadmissible shock.

# + [markdown] slideshow={"slide_type": "subslide"}
# Turning to the mathematical definitions we need an entropy function $\eta(q)$ and an entropy flux $\psi(q)$ that satisfies a new integral conservation law defined by
# $$
#     \int_{x_1}^{x_2} \eta(q(x,t_2)) dx = \int^{x_2}_{x_1} \eta(q(x,t_1)) dx + \int^{t_2}_{t_1} \psi(q(x_1, t)) dt - \int^{t_2}_{t_1} \psi(q(x_2, t)) dt.
# $$
# This holds whenever $q(x,t)$ is smooth.  When $q(x,t)$ is not smooth we expect this equation not to hold but instead turns into an inequality statement
# $$
#     \int_{x_1}^{x_2} \eta(q(x,t_2)) dx \leq \int^{x_2}_{x_1} \eta(q(x,t_1)) dx + \int^{t_2}_{t_1} \psi(q(x_1, t)) dt - \int^{t_2}_{t_1} \psi(q(x_2, t)) dt.
# $$

# + [markdown] slideshow={"slide_type": "subslide"}
# As perhaps expected if $q(x,t)$ is indeed smooth it can be manipulated as we have done before to derive
# $$
#     \eta(q)_t + \psi(q)_x = 0.
# $$
# Note also that if $\eta$ and $\psi$ are smooth functions we can write
# $$
#     \eta'(q) q_t + \psi'(q) q_x = 0,
# $$
# which comparing this to the original quasilinear form 
# $$
#     q_t + f'(q) q_x = 0
# $$
# from before allows us to identify
# $$
#     \psi'(q) = \eta'(q) f'(q)
# $$
# that should hold.

# + [markdown] slideshow={"slide_type": "subslide"}
# Interestingly for systems of equations we still require only scalar functions for $\eta$ and $\psi$.  This makes sense physically as entropy is only a scalar function.  Now however when we take a derivative of $\eta$ for instance we have
# $$
#     \eta'(q) = \begin{bmatrix} \frac{\partial \eta}{\partial q_1} & \frac{\partial \eta}{\partial q_2} & \cdots & \frac{\partial \eta}{\partial q_1} \end{bmatrix}.
# $$

# + [markdown] slideshow={"slide_type": "subslide"}
# We can also relate the idea of vanishing viscosity solution to the idea of entropy functions.  Define the related viscous equation as
# $$
#     q^\epsilon_t + f(q^\epsilon)_x = \epsilon q^\epsilon_{xx}
# $$
# and the associated entropy equation as
# $$
#     \eta(q^\epsilon)_t + \psi(q^\epsilon)_x = \epsilon \eta'(q^\epsilon) q^\epsilon_{xx}
# $$
# we can write (for smooth solutions) the expression derived before as
# $$
#     \eta(q^\epsilon)_t + \psi(q^\epsilon)_x = \epsilon(\eta'(q^\epsilon) q^\epsilon_x)_x - \epsilon \eta''(q^\epsilon) (q^\epsilon_x)^2.
# $$

# + [markdown] slideshow={"slide_type": "subslide"}
# Integrating this equation over our usual rectangle leads to
# $$\begin{aligned}
#     \int^{x_2}_{x_1} \eta(q^\epsilon(x, t_2)) dx &= \int^{x_2}_{x_1} \eta(q^\epsilon(x,t_1) dx - \left ( \int^{t_2}_{t_1} \psi(q^\epsilon(x_2, t)) dt - \int^{t_2}_{t_1} \psi(q^\epsilon(x_1, t)) dt \right) \\
#     &+ \epsilon \int^{t_2}_{t_1} [\eta'(q^\epsilon(x_2,t)) q^\epsilon_x(x_2, t) - \eta'(q^\epsilon(x_1,t)) q^\epsilon_x(x_1, t)] dt \\
#     &- \epsilon \int^{t_2}_{t_1} \int^{x_2}_{x_1} \eta''(q^\epsilon) (q^\epsilon_x)^2 dx dt
# \end{aligned}$$
# The usual flux differences are here but there are two additional terms that are multiplied by $\epsilon$.  As $\epsilon \rightarrow 0$ one can show that an entropy integral relation results.

# + [markdown] slideshow={"slide_type": "subslide"}
# We can also derive a true weak form of the entropy condition as was done before.  A weak solution $q$ satisfied the **weak form of the entropy inequality** if
# $$
#     \int^\infty_0 \int^\infty_{-\infty} [\phi_t \eta(q) + \phi_x \psi(q)] dx dt + \int^\infty_{-\infty} \phi(x, 0) \eta(q(x,0)) dx \geq 0
# $$
# for all $\phi in C^1_0(\mathbb R \times \mathbb R)$ with $\phi(x,t) \geq 0$ $\forall x,t$.

# + [markdown] slideshow={"slide_type": "subslide"}
# Another useful form of the entropy conditions can be defined similar to the Rankine-Hugoniot conditions and are derived similar to as before where now they take the form
# $$
#     s(\eta(q_r) - \eta(q_\ell)) \leq \psi(q_r) - \psi(q_\ell).
# $$

# + [markdown] slideshow={"slide_type": "subslide"}
# #### Example:  Burgers' Equation
#
# As was discussed previously we can use 
# $$
#     \eta(u) = u^2
# $$ 
# with 
# $$
#     \psi(u) = \frac{2}{3} u^3.
# $$
# If we have a shock propagating at speed
# $$
#     s = \frac{1}{2} (u_\ell + u_r)
# $$
# then the entropy condition requires
# $$
#     \frac{1}{2} (u_\ell + u_r) (u^2_r - u^2_\ell) \geq \frac{2}{3} (u^3_r - u^3_\ell).
# $$
# This then leads to the expression
# $$
#     \frac{1}{6} (u_\ell - u_r)^3 \geq 0
# $$
# implying that the entropy condition is only satisfied if $u_\ell > u_r$ we have already noticed due to the Lax entropy condition.

# + [markdown] slideshow={"slide_type": "subslide"}
# ### The Kruzkov Entropies
#
# The previous analysis assumes that the entropy function should be strictly convex with $\eta'' > 0$ at all $q$.  We can relax this condition by utilizing an idea due to Kruzkov who suggested a family of entropy functions rather than just one.  This is generalized to
# $$
#     \eta_k(q) = |q - k|, \quad \psi_k(q) = \text{sign}(q-k)[f(q) - f(k)]
# $$
# where $k \in \mathbb R$.  Note that ever $\eta_k$ is a piecewise defined, linear function of $q$ but that the higher derivatives are not neccesarily continuous.  We also see that $\eta''(q) = \delta(q - k)$ because of this.

# + [markdown] slideshow={"slide_type": "subslide"}
# ### Long-Term Behavior and N-Wave Decay
#
# We have seen before that for longer term behavior of solutions to both Burgers' equation and traffic-flow that we see an N-wave like solution that appears to be decaying.  This is indicative of the idea that entropy increases over a shock and information is "lost".

StarcoderdataJupyterScriptsDedupFiltered