Convergence of numerical solutions to weak solutions of conservation laws

This worksheet follows the discussion in [LeVeque, 1992] Chapter 12. Solutions of hyperbolic PDEs develop discontinuities, so we have to consider weak solutions that satisfy an integral conservation law. We'll consider the inviscid Burgers' equation:

$$q_t + \left(\frac{1}{2} q^2\right)_x = 0$$

discretized by two simple finite difference methods. The first is simply the upwind method applied to Burgers' equation:

$$Q^{n+1}_j = Q^n_j - \frac{\Delta t}{2\Delta x} \left((Q^n_j)^2 - (Q^n_{j-1})^2\right)$$

The second discretization is obtained by writing Burgers' equation in quasilinear form:

$$u_t + uu_x = 0$$

and using an upwind difference:

$$Q^{n+1}_j = Q^n_j - \frac{\Delta t}{\Delta x} Q^n_j\left(Q^n_j - Q^n_{j-1}\right)$$

The first method is conservative, since the fluxes will cancel out if we sum over the whole grid. The second method is not conservative. Let's implement them and see what happens.

In [14]:
import numpy as np
import matplotlib.pyplot as plt
import time, sys
import matplotlib
from matplotlib import animation
from JSAnimation import IPython_display
In [15]:
def take_step(q,dt,dx,method):
    if method == 'conservative':
        q[1:-1] = q[1:-1] - 0.5*dt/dx * ( q[1:-1]**2 - q[:-2]**2 )
    elif method == 'non-conservative':
        q[1:-1] = q[1:-1] - dt/dx * q[1:-1] * (q[1:-1] - q_nc[0:-2]) 
    return q
In [17]:
N=100;  # Number of grid points
plot_interval = N/10  # Draw a plot every n time steps

nghost = 1; N2=N+2*nghost;  # Ghost cells used for implementing boundary conditions
t=0.; T=0.5; # Initial and final time

q0 = np.sin(2*np.pi*x) + 1.5   # Smooth initial data
q_c = q0.copy();
q_nc= q0.copy();


# Set up plotting
fig = plt.figure()
axes = fig.add_subplot(111)
line1, = axes.plot(x,q0,lw=2)
line2, = axes.plot(x,q0,'r--',lw=2)

frames1 = [q0.copy()]
frames2 = [q0.copy()]

i = 0
while t<T:
    if dt>T-t: dt=T-t
    q_c  = take_step(q_c ,dt,dx,'conservative')
    q_nc = take_step(q_nc,dt,dx,'non-conservative')

    for q in [q_c, q_nc]:  # Periodic boundary conditions
        q[0:nghost]  = q[-2*nghost:-nghost] 
        q[-nghost:]  = q[nghost:2*nghost]
    t = t + dt
    i = i + 1
    if (i % plot_interval) == 0:
def plot_frame(i):

matplotlib.animation.FuncAnimation(fig, plot_frame, frames=len(frames1), interval=200)

Once Loop Reflect