In [3]:

```
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
```

Out[3]:

$$
\newcommand{Dx}{\Delta x}
\newcommand{Dt}{\Delta t}
\newcommand{imh}{{i-1/2}}
\newcommand{iph}{{i+1/2}}
\newcommand{\DQ}{\Delta Q}
\DeclareMathOperator{\sgn}{sgn}
$$

The methods we have used so far (the *upwind method* and the *Lax-Friedrichs method*) are both dissipative. This means that over time they artificially smear out the solution -- especially shocks. Furthermore, both of these methods are only *first order accurate*, meaning that if we reduce the values of $\Dt$ and $\Dx$ by a factor of two, the overall error decreases (only) by a factor of two. We can do better.

The first step in improving our accuracy is to use a more accurate representation of $q(x)$ at each step. Instead of assuming that $q$ is piecewise-constant, we will now approximate it by a piecewise-linear function:

$$q(x) = Q^n_i + \sigma^n_i (x-x_i).$$Here $\sigma_i$ represents the slope in cell $i$. The most obvious choice to ensure that this results in a second-order accurate approximation is to take the centered approximation

$$\sigma^n_i = \frac{Q^n_{i+1} - Q^n_{i-1}}{2\Dx}.$$We use this to obtain values at the cell interfaces:

\begin{align} q^+_\imh & = Q_i - \sigma_i \frac{\Dx}{2} \\ q^-_\iph & = Q_i + \sigma_i \frac{\Dx}{2}. \end{align}We'll use these interface values to approximate the flux, based on the **Lax-Friedrichs flux**:

This provides second-order accuracy in space. We also need to make the method second-order accurate in time. We can do so by using a second-order Runge--Kutta method. Then the full method is

\begin{align} Q^*_i & = Q^n_i - \frac{\Dt}{\Dx}\left( F^n_\iph - F^n_\imh \right) \\ Q^{n+1}_i & = \frac{1}{2} Q^n_i + \frac{1}{2}\left( Q^*_i - \frac{\Dt}{\Dx}\left( F^*_\iph - F^*_\imh \right) \right) \\ \end{align}In [4]:

```
import sys
sys.path.append('./util')
from ianimate import ianimate
import numpy as np
```

In [5]:

```
def f(q):
return q*(1.0-q)
m = 100 # number of points
dx = 1./m # Size of 1 grid cell
x = np.arange(-3*dx/2, 1.+5*dx/2, dx)
t = 0. # Initial time
T = 0.5 # Final time
dt = 0.9 * dx # Time step
Q = 0.9*np.exp(-100*(x-0.5)**2)
Qnew = np.zeros(m+4)
Qstar = np.zeros(m+4)
sigma = np.zeros(m+4)
F = np.zeros(m+4)
QQ = [Q]
while t < T:
sigma[1:-1] = (Q[2:] - Q[:-2])/(2.0*dx)
qplus = Q[1:-1] - sigma[1:-1] * dx/2.0 # q^+_{i-1/2}
qminus = Q[:-2] + sigma[:-2] * dx/2.0 # q^-_{i-1/2}
F[1:-1] = 0.5*(f(qplus)+f(qminus) - dx/dt*(qplus-qminus) )# F_{i-1/2}
Qstar[2:-2] = Q[2:-2] - dt/dx*(F[3:-1]-F[2:-2])
Qstar[0:2] = Qstar[2]
Qstar[-2:] = Qstar[-3]
sigma[1:-1] = (Qstar[2:] - Qstar[:-2])/(2.0*dx)
qplus = Qstar[1:-1] - sigma[1:-1] * dx/2.0 # q^+_{i-1/2}
qminus = Qstar[:-2] + sigma[:-2] * dx/2.0 # q^-_{i-1/2}
F[1:-1] = 0.5*(f(qplus)+f(qminus) - dx/dt*(qplus-qminus) )# F_{i-1/2}
Qnew[2:-2] = 0.5*Q[2:-2] + 0.5*(Qstar[2:-2] - dt/dx*(F[3:-1]-F[2:-2]))
Q = Qnew.copy()
Q[0:2] = Q[2]
Q[-2:] = Q[-3]
t = t + dt
QQ.append(Q)
ianimate(x,QQ)
```

Out[5]: