In [3]:
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
Out[3]:
Content provided under a Creative Commons Attribution license, CC-BY 4.0; code under MIT License. (c)2014 David I. Ketcheson
version 0.1 - May 2014

High-resolution methods

$$ \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.

Reducing dissipation

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:

$$F_\imh = \frac{1}{2} \left( f(q^-_\imh) + f(q^+_\imh) - \frac{\Dt}{\Dx} (q^+_\imh - q^-_\imh)\right)$$

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]:


Once Loop Reflect