Notebook

Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 L.A. Barba, C.D. Cooper, G.F. Forsyth.¶

Riding the wave¶

Numerical schemes for hyperbolic PDEs¶

Welcome back! This is the second notebook of Riding the wave: Convection problems, the third module of "Practical Numerical Methods with Python".

The first notebook of this module discussed conservation laws and developed the non-linear traffic equation. We learned about the effect of the wave speed on the stability of the numerical method, and on the CFL number. We also realized that the forward-time/backward-space difference scheme really has many limitations: it cannot deal with wave speeds that move in more than one direction. It is also first-order accurate in space and time, which often is just not good enough. This notebook will introduce some new numerical schemes for conservation laws, continuing with the traffic-flow problem as motivation.

Red light!¶

Let's explore the behavior of different numerical schemes for a moving shock wave. In the context of the traffic-flow model of the previous notebook, imagine a very busy road and a red light at $x=4$ . Cars accumulate quickly in the front, where we have the maximum allowed density of cars between $x=3$ and $x=4$ , and there is an incoming traffic of 50% the maximum allowed density $(\rho = 0.5\rho_{\rm max})$ .

Mathematically, this is:

$\begin{equation} \rho(x,0) = \left\{ \begin{array}{cc} 0.5 \rho_{\rm max} & 0 \leq x < 3 \\ \rho_{\rm max} & 3 \leq x \leq 4 \\ \end{array} \right. \end{equation}$

Let's find out what the initial condition looks like.

In [1]:

import numpy
from matplotlib import pyplot
%matplotlib inline

In [2]:

# Set the font family and size to use for Matplotlib figures.
pyplot.rcParams['font.family'] = 'serif'
pyplot.rcParams['font.size'] = 16

In [3]:

def rho_red_light(x, rho_max):
    """
    Computes the "red light" initial condition with shock.
    
    Parameters
    ----------
    x : numpy.ndaray
        Locations on the road as a 1D array of floats.
    rho_max : float
        The maximum traffic density allowed.
    
    Returns
    -------
    rho : numpy.ndarray
        The initial car density along the road
        as a 1D array of floats.
    """
    rho = rho_max * numpy.ones_like(x)
    mask = numpy.where(x < 3.0)
    rho[mask] = 0.5 * rho_max
    return rho

In [4]:

# Set parameters.
nx = 81  # number of locations on the road
L = 4.0  # length of the road
dx = L / (nx - 1)  # distance between two consecutive locations
nt = 40  # number of time steps to compute
rho_max = 10.0  # maximum taffic density allowed
u_max = 1.0  # maximum speed traffic

# Get the road locations.
x = numpy.linspace(0.0, L, num=nx)

# Compute the initial traffic density.
rho0 = rho_red_light(x, rho_max)

In [5]:

# Plot the initial traffic density.
fig = pyplot.figure(figsize=(6.0, 4.0))
pyplot.xlabel(r'$x$')
pyplot.ylabel(r'$\rho$')
pyplot.grid()
line = pyplot.plot(x, rho0,
                   color='C0', linestyle='-', linewidth=2)[0]
pyplot.xlim(0.0, L)
pyplot.ylim(4.0, 11.0)
pyplot.tight_layout()

The question we would like to answer is: How will cars accumulate at the red light?

We will solve this problem using different numerical schemes, to see how they perform. These schemes are:

Lax-Friedrichs
Lax-Wendroff
MacCormack

Before we do any coding, let's think about the equation a little bit. The wave speed $u_{\rm wave}$ is $-1$ for $\rho = \rho_{\rm max}$ and $\rho \leq \rho_{\rm max}/2$ , making all velocities negative. We should see a solution moving left, maintaining the shock geometry.

squarewave

Figure 1. The exact solution is a shock wave moving left.¶

Now to some coding! First, let's define some useful functions and prepare to make some nice animations later.

In [6]:

def flux(rho, u_max, rho_max):
    """
    Computes the traffic flux F = V * rho.
    
    Parameters
    ----------
    rho : numpy.ndarray
       Traffic density along the road as a 1D array of floats.
    u_max : float
        Maximum speed allowed on the road.
    rho_max : float
        Maximum car density allowed on the road.
    
    Returns
    -------
    F : numpy.ndarray
        The traffic flux along the road as a 1D array of floats.
    """
    F = rho * u_max * (1.0 - rho / rho_max)
    return F

Before we investigate different schemes, let's create the function to update the Matplotlib figure during the animation.

In [7]:

from matplotlib import animation
from IPython.display import HTML

In [8]:

def update_plot(n, rho_hist):
    """
    Update the line y-data of the Matplotlib figure.
    
    Parameters
    ----------
    n : integer
        The time-step index.
    rho_hist : list of numpy.ndarray objects
        The history of the numerical solution.
    """
    fig.suptitle('Time step {:0>2}'.format(n))
    line.set_ydata(rho_hist[n])

Lax-Friedrichs scheme¶

Recall the conservation law for vehicle traffic, resulting in the following equation for the traffic density:

$\begin{equation} \frac{\partial \rho}{\partial t} + \frac{\partial F}{\partial x} = 0 \end{equation}$

$F$ is the traffic flux, which in the linear traffic-speed model is given by:

$\begin{equation} F = \rho u_{\rm max} \left(1-\frac{\rho}{\rho_{\rm max}}\right) \end{equation}$

In the time variable, the natural choice for discretization is always a forward-difference formula; time invariably moves forward!

$\begin{equation} \frac{\partial \rho}{\partial t}\approx \frac{1}{\Delta t}( \rho_i^{n+1}-\rho_i^n ) \end{equation}$

As is usual, the discrete locations on the 1D spatial grid are denoted by indices $i$ and the discrete time instants are denoted by indices $n$ .

In a convection problem, using first-order discretization in space leads to excessive numerical diffusion (as you probably observed in Lesson 1 of Module 2). The simplest approach to get second-order accuracy in space is to use a central difference:

$\begin{equation} \frac{\partial F}{\partial x} \approx \frac{1}{2\Delta x}( F_{i+1}-F_{i-1}) \end{equation}$

But combining these two choices for time and space discretization in the convection equation has catastrophic results! The "forward-time, central scheme" (FTCS) is unstable. (Go on: try it; you know you want to!)

The Lax-Friedrichs scheme was proposed by Lax (1954) as a clever trick to stabilize the forward-time, central scheme. The idea was to replace the solution value at $\rho^n_i$ by the average of the values at the neighboring grid points. If we do that replacement, we get the following discretized equation:

$\begin{equation} \frac{\rho_i^{n+1}-\frac{1}{2}(\rho^n_{i+1}+\rho^n_{i-1})}{\Delta t} = -\frac{F^n_{i+1}-F^n_{i-1}}{2 \Delta x} \end{equation}$

Take a careful look: the difference formula no longer uses the value at $\rho^n_i$ to obtain $\rho^{n+1}_i$ . The stencil of the Lax-Friedrichs scheme is slightly different than that for the forward-time, central scheme.

Stencil of the forward-time central scheme

Figure 2. Stencil of the forward-time/central scheme.¶

Stencil of the Lax-Friedrichs scheme

Figure 3. Stencil of the Lax-Friedrichs scheme.¶

This numerical discretization is stable. Unfortunately, substituting $\rho^n_i$ by the average of its neighbors introduces a first-order error. Nice try, Lax!

To implement the scheme in code, we need to isolate the value at the next time step, $\rho^{n+1}_i$ , so we can write a time-stepping loop:

$\begin{equation} \rho_i^{n+1} = \frac{1}{2}(\rho^n_{i+1}+\rho^n_{i-1}) - \frac{\Delta t}{2 \Delta x}(F^n_{i+1}-F^n_{i-1}) \end{equation}$

The function below implements Lax-Friedrichs for our traffic model. All the schemes in this notebook are wrapped in their own functions to help with displaying animations of the results. This is also good practice for developing modular, reusable code.

In order to display animations, we're going to hold the results of each time step in the variable rho, a 2D array. The resulting array rho_n has nt rows and nx columns.

In [9]:

def lax_friedrichs(rho0, nt, dt, dx, bc_values, *args):
    """
    Computes the traffic density on the road 
    at a certain time given the initial traffic density.
    Integration using Lax-Friedrichs scheme.
    
    Parameters
    ----------
    rho0 : numpy.ndarray
        The initial traffic density along the road
        as a 1D array of floats.
    nt : integer
        The number of time steps to compute.
    dt : float
        The time-step size to integrate.
    dx : float
        The distance between two consecutive locations.
    bc_values : 2-tuple of floats
        The value of the density at the first and last locations.
    args : list or tuple
        Positional arguments to be passed to the flux function.
    
    Returns
    -------
    rho_hist : list of numpy.ndarray objects
        The history of the car density along the road.
    """
    rho_hist = [rho0.copy()]
    rho = rho0.copy()
    for n in range(nt):
        # Compute the flux.
        F = flux(rho, *args)
        # Advance in time using Lax-Friedrichs scheme.
        rho[1:-1] = (0.5 * (rho[:-2] + rho[2:]) -
                     dt / (2.0 * dx) * (F[2:] - F[:-2]))
        # Set the value at the first location.
        rho[0] = bc_values[0]
        # Set the value at the last location.
        rho[-1] = bc_values[1]
        # Record the time-step solution.
        rho_hist.append(rho.copy())
    return rho_hist

Lax-Friedrichs with $\frac{\Delta t}{\Delta x}=1$ ¶

We are now all set to run! First, let's try with CFL=1

In [10]:

# Set the time-step size based on CFL limit.
sigma = 1.0
dt = sigma * dx / u_max  # time-step size

# Compute the traffic density at all time steps.
rho_hist = lax_friedrichs(rho0, nt, dt, dx, (rho0[0], rho0[-1]),
                          u_max, rho_max)

In [11]:

# Create an animation of the traffic density.
anim = animation.FuncAnimation(fig, update_plot,
                               frames=nt, fargs=(rho_hist,),
                               interval=100)
# Display the video.
HTML(anim.to_html5_video())

Out[11]:

Think¶

What do you see in the animation above? How does the numerical solution compare with the exact solution (a left-traveling shock wave)?
What types of errors do you think we see?
What do you think of the Lax-Friedrichs scheme, so far?

Lax-Friedrichs with $\frac{\Delta t}{\Delta x} = 0.5$ ¶

Would the solution improve if we use smaller time steps? Let's check that!

In [12]:

# Set the time-step size based on CFL limit.
sigma = 0.5
dt = sigma * dx / u_max  # time-step size

# Compute the traffic density at all time steps.
rho_hist = lax_friedrichs(rho0, nt, dt, dx, (rho0[0], rho0[-1]),
                          u_max, rho_max)

In [13]:

# Create an animation of the traffic density.
anim = animation.FuncAnimation(fig, update_plot,
                               frames=nt, fargs=(rho_hist,),
                               interval=100)
# Display the video.
HTML(anim.to_html5_video())

Out[13]:

Dig deeper¶

Notice the strange "staircase" behavior on the leading edge of the wave? You may be interested to learn more about this: a feature typical of what is sometimes called "odd-even decoupling." Last year we published a collection of lessons in Computational Fluid Dynamics, called CFD Python, where we discuss odd-even decoupling.

How does this solution compare with the previous one, where the Courant number was $\frac{\Delta t}{\Delta x}=1$ ?

Lax-Wendroff scheme¶

The Lax-Friedrichs method uses a clever trick to stabilize the central difference in space for convection, but loses an order of accuracy in doing so. First-order methods are just not good enough for convection problems, especially when you have sharp gradients (shocks).

The Lax-Wendroff (1960) method was the first scheme ever to achieve second-order accuracy in both space and time. It is therefore a landmark in the history of computational fluid dynamics.

To develop the Lax-Wendroff scheme, we need to do a bit of work. Sit down, grab a notebook and grit your teeth. We want you to follow this derivation in your own hand. It's good for you! Start with the Taylor series expansion (in the time variable) about $\rho^{n+1}$ :

$\begin{equation} \rho^{n+1} = \rho^n + \frac{\partial\rho^n}{\partial t} \Delta t + \frac{(\Delta t)^2}{2}\frac{\partial^2\rho^n}{\partial t^2} + \ldots \end{equation}$

For the conservation law with $F=F(\rho)$ , and using our beloved chain rule, we can write:

$\begin{equation} \frac{\partial \rho}{\partial t} = -\frac{\partial F}{\partial x} = -\frac{\partial F}{\partial \rho} \frac{\partial \rho}{\partial x} = -J \frac{\partial \rho}{\partial x} \end{equation}$

where

$\begin{equation} J = \frac{\partial F}{\partial \rho} = u _{\rm max} \left(1-2\frac{\rho}{\rho_{\rm max}} \right) \end{equation}$

is the Jacobian for the traffic model. Next, we can do a little trickery:

$\begin{equation} \frac{\partial F}{\partial t} = \frac{\partial F}{\partial \rho} \frac{\partial \rho}{\partial t} = J \frac{\partial \rho}{\partial t} = -J \frac{\partial F}{\partial x} \end{equation}$

In the last step above, we used the differential equation of the traffic model to replace the time derivative by a spatial derivative. These equivalences imply that

$\begin{equation} \frac{\partial^2\rho}{\partial t^2} = \frac{\partial}{\partial x} \left( J \frac{\partial F}{\partial x} \right) \end{equation}$

Let's use all this in the Taylor expansion:

$\begin{equation} \rho^{n+1} = \rho^n - \frac{\partial F^n}{\partial x} \Delta t + \frac{(\Delta t)^2}{2} \frac{\partial}{\partial x} \left(J\frac{\partial F^n}{\partial x} \right)+ \ldots \end{equation}$

We can now reorganize this and discretize the spatial derivatives with central differences to get the following discrete equation:

$\begin{equation} \frac{\rho_i^{n+1} - \rho_i^n}{\Delta t} = -\frac{F^n_{i+1}-F^n_{i-1}}{2 \Delta x} + \frac{\Delta t}{2} \left(\frac{(J \frac{\partial F}{\partial x})^n_{i+\frac{1}{2}}-(J \frac{\partial F}{\partial x})^n_{i-\frac{1}{2}}}{\Delta x}\right) \end{equation}$

Now, approximate the rightmost term (inside the parenthesis) in the above equation as follows:

$\begin{equation} \frac{J^n_{i+\frac{1}{2}}\left(\frac{F^n_{i+1}-F^n_{i}}{\Delta x}\right)-J^n_{i-\frac{1}{2}}\left(\frac{F^n_i-F^n_{i-1}}{\Delta x}\right)}{\Delta x}\end{equation}$

Then evaluate the Jacobian at the midpoints by using averages of the points on either side:

$\begin{equation}\frac{\frac{1}{2 \Delta x}(J^n_{i+1}+J^n_i)(F^n_{i+1}-F^n_i)-\frac{1}{2 \Delta x}(J^n_i+J^n_{i-1})(F^n_i-F^n_{i-1})}{\Delta x}.\end{equation}$

Our equation now reads:

$\begin{align} &\frac{\rho_i^{n+1} - \rho_i^n}{\Delta t} = -\frac{F^n_{i+1}-F^n_{i-1}}{2 \Delta x} + \cdots \\ \nonumber &+ \frac{\Delta t}{4 \Delta x^2} \left( (J^n_{i+1}+J^n_i)(F^n_{i+1}-F^n_i)-(J^n_i+J^n_{i-1})(F^n_i-F^n_{i-1})\right) \end{align}$

Solving for $\rho_i^{n+1}$ :

$\begin{align} &\rho_i^{n+1} = \rho_i^n - \frac{\Delta t}{2 \Delta x} \left(F^n_{i+1}-F^n_{i-1}\right) + \cdots \\ \nonumber &+ \frac{(\Delta t)^2}{4(\Delta x)^2} \left[ (J^n_{i+1}+J^n_i)(F^n_{i+1}-F^n_i)-(J^n_i+J^n_{i-1})(F^n_i-F^n_{i-1})\right] \end{align}$

with

$\begin{equation}J^n_i = \frac{\partial F}{\partial \rho} = u_{\rm max} \left(1-2\frac{\rho^n_i}{\rho_{\rm max}} \right).\end{equation}$

Lax-Wendroff is a little bit long. Remember that you can use \ slashes to split up a statement across several lines. This can help make code easier to parse (and also easier to debug!).

In [14]:

def jacobian(rho, u_max, rho_max):
    """
    Computes the Jacobian for our traffic model.
    
    Parameters
    ----------
    rho : numpy.ndarray
       Traffic density along the road as a 1D array of floats.
    u_max : float
        Maximum speed allowed on the road.
    rho_max : float
        Maximum car density allowed on the road.
    
    Returns
    -------
    J : numpy.ndarray
        The Jacobian as a 1D array of floats.
    """
    J = u_max * (1.0 - 2.0 * rho / rho_max)
    return J

In [15]:

def lax_wendroff(rho0, nt, dt, dx, bc_values, *args):
    """
    Computes the traffic density on the road 
    at a certain time given the initial traffic density.
    Integration using Lax-Wendroff scheme.
    
    Parameters
    ----------
    rho0 : numpy.ndarray
        The initial traffic density along the road
        as a 1D array of floats.
    nt : integer
        The number of time steps to compute.
    dt : float
        The time-step size to integrate.
    dx : float
        The distance between two consecutive locations.
    bc_values : 2-tuple of floats
        The value of the density at the first and last locations.
    args : list or tuple
        Positional arguments to be passed to the
        flux and Jacobien functions.
    
    Returns
    -------
    rho_hist : list of numpy.ndarray objects
        The history of the car density along the road.
    """
    rho_hist = [rho0.copy()]
    rho = rho0.copy()
    for n in range(nt):
        # Compute the flux.
        F = flux(rho, *args)
        # Compute the Jacobian.
        J = jacobian(rho, *args)
        # Advance in time using Lax-Wendroff scheme.
        rho[1:-1] = (rho[1:-1] -
                     dt / (2.0 * dx) * (F[2:] - F[:-2]) +
                     dt**2 / (4.0 * dx**2) *
                     ((J[1:-1] + J[2:]) * (F[2:] - F[1:-1]) -
                      (J[:-2] + J[1:-1]) * (F[1:-1] - F[:-2])))
        # Set the value at the first location.
        rho[0] = bc_values[0]
        # Set the value at the last location.
        rho[-1] = bc_values[1]
        # Record the time-step solution.
        rho_hist.append(rho.copy())
    return rho_hist

Now that's we've defined a function for the Lax-Wendroff scheme, we can use the same procedure as above to animate and view our results.

Lax-Wendroff with $\frac{\Delta t}{\Delta x}=1$ ¶

In [16]:

# Set the time-step size based on CFL limit.
sigma = 1.0
dt = sigma * dx / u_max  # time-step size

# Compute the traffic density at all time steps.
rho_hist = lax_wendroff(rho0, nt, dt, dx, (rho0[0], rho0[-1]),
                        u_max, rho_max)

In [17]:

# Create an animation of the traffic density.
anim = animation.FuncAnimation(fig, update_plot,
                               frames=nt, fargs=(rho_hist,),
                               interval=100)
# Display the video.
HTML(anim.to_html5_video())

Out[17]:

Interesting! The Lax-Wendroff method captures the sharpness of the shock much better than the Lax-Friedrichs scheme, but there is a new problem: a strange wiggle appears right at the tail of the shock. This is typical of many second-order methods: they introduce numerical oscillations where the solution is not smooth. Bummer.

Lax-Wendroff with $\frac{\Delta t}{\Delta x} =0.5$ ¶

How do the oscillations at the shock front vary with changes to the CFL condition? You might think that the solution will improve if you make the time step smaller ... let's see.

In [18]:

# Set the time-step size based on CFL limit.
sigma = 0.5
dt = sigma * dx / u_max  # time-step size

# Compute the traffic density at all time steps.
rho_hist = lax_wendroff(rho0, nt, dt, dx, (rho0[0], rho0[-1]),
                        u_max, rho_max)

In [19]:

# Create an animation of the traffic density.
anim = animation.FuncAnimation(fig, update_plot,
                               frames=nt, fargs=(rho_hist,),
                               interval=100)
# Display the video.
HTML(anim.to_html5_video())

Out[19]:

Eek! The numerical oscillations got worse. Double bummer!

Why do we observe oscillations with second-order methods? This is a question of fundamental importance!

MacCormack Scheme¶

The numerical oscillations that you observed with the Lax-Wendroff method on the traffic model can become severe in some problems. But actually the main drawback of the Lax-Wendroff method is having to calculate the Jacobian in every time step. With more complicated equations (like the Euler equations), calculating the Jacobian is a large computational expense.

Robert W. MacCormack introduced the first version of his now-famous method at the 1969 AIAA Hypervelocity Impact Conference, held in Cincinnati, Ohio, but the paper did not at first catch the attention of the aeronautics community. The next year, however, he presented at the 2nd International Conference on Numerical Methods in Fluid Dynamics at Berkeley. His paper there (MacCormack, 1971) was a landslide. MacCormack got a promotion and continued to work on applications of his method to the compressible Navier-Stokes equations. In 1973, NASA gave him the prestigious H. Julian Allen award for his work.

The MacCormack scheme is a two-step method, in which the first step is called a predictor and the second step is called a corrector. It achieves second-order accuracy in both space and time. One version is as follows:

$\begin{equation} \rho^*_i = \rho^n_i - \frac{\Delta t}{\Delta x} (F^n_{i+1}-F^n_{i}) \ \ \ \ \ \ \text{(predictor)} \end{equation}$

$\begin{equation} \rho^{n+1}_i = \frac{1}{2} (\rho^n_i + \rho^*_i - \frac{\Delta t}{\Delta x} (F^*_i - F^{*}_{i-1})) \ \ \ \ \ \ \text{(corrector)} \end{equation}$

If you look closely, it appears like the first step is a forward-time/forward-space scheme, and the second step is like a forward-time/backward-space scheme (these can also be reversed), averaged with the first result. What is so cool about this? You can compute problems with left-running waves and right-running waves, and the MacCormack scheme gives you a stable method (subject to the CFL condition). Nice! Let's try it.

In [20]:

def maccormack(rho0, nt, dt, dx, bc_values, *args):
    """
    Computes the traffic density on the road 
    at a certain time given the initial traffic density.
    Integration using MacCormack scheme.
    
    Parameters
    ----------
    rho0 : numpy.ndarray
        The initial traffic density along the road
        as a 1D array of floats.
    nt : integer
        The number of time steps to compute.
    dt : float
        The time-step size to integrate.
    dx : float
        The distance between two consecutive locations.
    bc_values : 2-tuple of floats
        The value of the density at the first and last locations.
    args : list or tuple
        Positional arguments to be passed to the flux function.
    
    Returns
    -------
    rho_hist : list of numpy.ndarray objects
        The history of the car density along the road.
    """
    rho_hist = [rho0.copy()]
    rho = rho0.copy()
    rho_star = rho.copy()
    for n in range(nt):
        # Compute the flux.
        F = flux(rho, *args)
        # Predictor step of the MacCormack scheme.
        rho_star[1:-1] = (rho[1:-1] -
                          dt / dx * (F[2:] - F[1:-1]))
        # Compute the flux.
        F = flux(rho_star, *args)
        # Corrector step of the MacCormack scheme.
        rho[1:-1] = 0.5 * (rho[1:-1] + rho_star[1:-1] -
                           dt / dx * (F[1:-1] - F[:-2]))
        # Set the value at the first location.
        rho[0] = bc_values[0]
        # Set the value at the last location.
        rho[-1] = bc_values[1]
        # Record the time-step solution.
        rho_hist.append(rho.copy())
    return rho_hist

MacCormack with $\frac{\Delta t}{\Delta x} = 1$ ¶

In [21]:

# Set the time-step size based on CFL limit.
sigma = 1.0
dt = sigma * dx / u_max  # time-step size

# Compute the traffic density at all time steps.
rho_hist = maccormack(rho0, nt, dt, dx, (rho0[0], rho0[-1]),
                      u_max, rho_max)

In [22]:

# Create an animation of the traffic density.
anim = animation.FuncAnimation(fig, update_plot,
                               frames=nt, fargs=(rho_hist,),
                               interval=100)
# Display the video.
HTML(anim.to_html5_video())

Out[22]:

MacCormack with $\frac{\Delta t}{\Delta x}= 0.5$ ¶

Once again, we ask: how does the CFL number affect the errors? Which one gives better results? You just have to try it.

In [23]:

# Set the time-step size based on CFL limit.
sigma = 0.5
dt = sigma * dx / u_max  # time-step size

# Compute the traffic density at all time steps.
rho_hist = maccormack(rho0, nt, dt, dx, (rho0[0], rho0[-1]),
                      u_max, rho_max)

In [24]:

# Create an animation of the traffic density.
anim = animation.FuncAnimation(fig, update_plot,
                               frames=nt, fargs=(rho_hist,),
                               interval=100)
# Display the video.
HTML(anim.to_html5_video())

Out[24]:

Dig Deeper¶

You can also obtain a MacCormack scheme by reversing the predictor and corrector steps. For shocks, the best resolution will occur when the difference in the predictor step is in the direction of propagation. Try it out! Was our choice here the ideal one? In which case is the shock better resolved?

Challenge task¶

In the red light problem, $\rho \geq \rho_{\rm max}/2$ , making the wave speed negative at all points . You might be wondering why we introduced these new methods; couldn't we have just used a forward-time/forward-space scheme? But, what if $\rho_{\rm in} < \rho_{\rm max}/2$ ? Now, a whole region has negative wave speeds and forward-time/backward-space is unstable.

How do Lax-Friedrichs, Lax-Wendroff and MacCormack behave in this case? Try it out!
As you decrease $\rho_{\rm in}$ , what happens to the velocity of the shock? Why do you think that happens?

References¶

Peter D. Lax (1954), "Weak solutions of nonlinear hyperbolic equations and their numerical computation," Commun. Pure and Appl. Math., Vol. 7, pp. 159–193.
Peter D. Lax and Burton Wendroff (1960), "Systems of conservation laws," Commun. Pure and Appl. Math., Vol. 13, pp. 217–237.
R. W. MacCormack (1969), "The effect of viscosity in hypervelocity impact cratering," AIAA paper 69-354. Reprinted on Journal of Spacecraft and Rockets, Vol. 40, pp. 757–763 (2003). Also on Frontiers of Computational Fluid Dynamics, edited by D. A. Caughey, M. M. Hafez (2002), chapter 2: read on Google Books.
R. W. MacCormack (1971), "Numerical solution of the interaction of a shock wave with a laminar boundary layer," Proceedings of the 2nd Int. Conf. on Numerical Methods in Fluid Dynamics, Lecture Notes in Physics, Vol. 8, Springer, Berlin, pp. 151–163.

The cell below loads the style of the notebook.¶

In [25]:

from IPython.core.display import HTML
css_file = '../../styles/numericalmoocstyle.css'
HTML(open(css_file, 'r').read())

Out[25]: