• «
• 1
• ...
• 65
• 66
• 67
• 68
• 69
• 70
• 71
• 72
• 73
• 74
• 75
• 76
• 77
• 78
• 79
• »

8.7 Lax-Wendroff Schemes

These schemes were proposed in 1960 by P.D. Lax and B. Wendroff [15] for solving, approximately, systems of hyperbolic conservation laws on the generic form given in (8.3).

A large class of numerical methods for solving (8.3) are the so-called conservative methods: $$ \begin{equation} u_j^{n+1} = u_j^{n} + \frac{\Delta t}{\Delta x} \left ( F_{j-1/2} - F_{j+1/2} \right ) \tag{8.41} \end{equation} $$

Linear advection.

The Lax–Wendroff method belongs to the class of conservative schemes (8.3) and can be derived in various ways. For simplicity, we will derive the method by using a simple model equation for (8.3), namely the linear advection equation with \( F(u)=a\,u \) as in (8.1), where \( a \) is a constant propagation velocity. The Lax-Wendroff outset is a Taylor approximation of \( u_j^{n+1} \): $$ \begin{equation} u_j^{n+1}=u_j^n+\Delta t \frac{\partial u}{\partial t}\Bigg|_j^n + \frac{(\Delta t)}{2}\frac{\partial^2u}{\partial t^2}\Bigg|_j^n+\cdots \tag{8.42} \end{equation} $$ From the differential equation (8.3) we get by differentiation $$ \begin{equation} \begin{array}{c} \dfrac{\partial u}{\partial t}\Bigg|_j^n = -a_0\dfrac{\partial u}{\partial x}\Bigg|_j^n \qquad \text{ and} \qquad \dfrac{\partial^2u}{\partial t^2}\Bigg|_j^n = a_0^2\dfrac{\partial^2u}{\partial x^2}\Bigg|_j^n \end{array} \tag{8.43} \end{equation} $$ Before substitution of (8.43) in the Taylor expansion (8.42) we approximate the spatial derivatives by central differences: $$ \begin{equation} \dfrac{\partial u}{\partial x}\Bigg|_j^n \approx \dfrac{u_{j+1}^n-u_{j-1}^n}{(2\Delta x)} \qquad \text{ and } \qquad \dfrac{\partial^2u}{\partial x^2}\Bigg|_j^n \approx \dfrac{u_{j+1}^n -2u_j^n + u_{j-1}^n}{(\Delta x)^2} \tag{8.44} \end{equation} $$ and then the Lax-Wendroff scheme follows by substitution: $$ \begin{equation} \tag{8.45} u_j^{n+1}=u_j^n-\frac{C}{2}\left( u_{j+1}^n-u_{j-1}^n \right) + \frac{C^2}{2}\left( u_{j+1}^n - 2u_j^n+u_{j-1}^n \right) \end{equation} $$ with the local truncation error \( T_j^n \): $$ \begin{equation} T_j^n = \frac{1}{6}\cdot \left[ (\Delta t)^2 \frac{\partial^3u}{\partial t^3} + a_0(\Delta x)^2 \frac{\partial^3u}{\partial x^3}\right]_j^n = O[(\Delta t)^2,(\Delta x)^2] \label{} \end{equation} $$ The resulting difference equation in (8.45) may also be formulated as: $$ \begin{equation} \tag{8.46} u_j^{n+1} = \frac{C}{2}(1+C)u_{j-1}^n + (1-C^2)u_j^n - \frac{C}{2}(1-C)u_{j+1}^n \end{equation} $$

The explicit Lax-Wendroff stencil is illustrated in Figure 112

An example of how to implement the Lax-Wendroff scheme is given as follows:

def lax_wendroff(u): 
    u[1:-1] = c/2.0*(1+c)*u[:-2] + (1-c**2)*u[1:-1] - c/2.0*(1-c)*u[2:]
    return u[1:-1]

8.7.1 Lax-Wendroff for non-linear systems of hyperbolic PDEs

Richtmyer Scheme. One of the earliest extensions of the scheme is the Richtmyer two-step Lax–Wendroff method, which is on the conservative form (8.41) with the numerical fluxes computed as follows: $$ \begin{align} &u_{j+1/2}^{n+1/2} = \frac{1}{2} \left (u_{j}^{n} + u_{j+1}^{n} \right ) + \frac{1}{2} \frac{\Delta t}{\Delta x} \left (F_j^n - F_{j+1}^n \right ) \tag{8.47}\\ & F_{j+1/2} = F(u_{j+1/2}^{n+1/2}) \tag{8.48} \end{align} $$

Lax-Wendroff two step.

A Lax-Wendroff two step method is outlined in the following. In the first step \( u(x, t) \) is evaluated at half time steps \( n + 1/2 \) and half grid points \( j+1/2 \). In the second step values at the next time step \( n + 1 \) are calculated using the data for \( n \) and \( n + 1/2 \).

First step: $$ \begin{align} & u_{j+1/2}^{n+1/2} = \frac{1}{2} \left (u_{j+1}^{n} + u_{j}^{n} \right ) - \frac{\Delta t}{2 \Delta x} \left (F(u_{j+1}^{n}) - F(u_{j}^{n}) \right ) \tag{8.49}\\ & u_{j-1/2}^{n+1/2} = \frac{1}{2} \left (u_{j}^{n} + u_{j-1}^{n} \right ) - \frac{\Delta t}{2 \Delta x} \left (F(u_{j}^{n}) - F(u_{j-1}^{n}) \right ) \tag{8.50} \end{align} $$

Second step: $$ \begin{equation} u_{j}^{n+1} = u_{j}^{n} - \frac{\Delta t}{\Delta x} \left ( F(u_{j+1/2}^{n+1/2}) - F(u_{j-1/2}^{n+1/2}) \right ) \tag{8.51} \end{equation} $$

Notice that for a linear flux \( F=a_0 \, u \), the two-step Lax-Wendroff method ((8.50) and (8.51)) may be shown to reduce to the one-step Lax-Wendroff method outlined in (8.45) or (8.46).

MacCormack Scheme.

A simpler and popular extension/variant of Lax-Wendroff schemes like in the previous section, is the MacCormack scheme [16]: $$ \begin{align} & u_j^p =u_j^n + \frac{\Delta t}{\Delta x} \left (F_{j}^{n} - F_{j+1}^{n} \right ) \nonumber \\ \tag{8.52} \\ & u_j^{n+1} =\frac{1}{2} \left (u_{j}^{n} + u_{j}^{p} \right ) + \frac{1}{2} \frac{\Delta t}{\Delta x} \left (F_{j-1}^{p} - F_{j}^{p} \right ) \nonumber \end{align} $$

where we have introduced the convention \( F_{j}^{p}=F(u_{j}^{p}) \).

Note that in the predictor step we employ the conservative formula (8.41) for a time \( \Delta t \) with forward differencing, i.e. . \( F_{j+1/2} = F^n_{j+1}=F(u^n_{j+1}) \). The corrector step may be interpreted as using (8.41) for a time \( \Delta t/2 \) with initial condition \( \frac{1}{2} \left (u_{j}^{n} + u_{j+1}^{p} \right ) \) and backward differencing.

Another MacCormack scheme may be obtained by reversing the predictor and corrector steps. Note that the MacCormack scheme (8.52) is not written in conservative form (8.41). However, it easy to express the scheme in conservative form by expressing the flux in (8.41) as: $$ \begin{equation} F_{j+1}^{m} = \frac{1}{2} \left (F_j^p + F_{j+1}^n \right ) \label{} \end{equation} $$

For a linear flux \( F(u) = a_0 \, u \), one may show that the MacCormack scheme in (8.52) reduces to a two-step scheme: $$ \begin{align} & u_j^p = u_j^n + C \left ( u_j^n - u_{j+1}^n \right ) \tag{8.53}\\ & u_j^{n+1} = \frac{1}{2} \left ( u_j^n + u_j^p \right )+ \frac{C}{2} \left ( u_{j-1}^p - u_{j}^p \right ) \tag{8.54} \end{align} $$

and substitution of (8.53) into (8.54) shows that the MacCormack scheme is identical to the Lax-Wendroff scheme (8.46) for the linear advection flux. A python implementation is given by:

def macCormack(u):
    up = u.copy()
    up[:-1] = u[:-1] - c*(u[1:]-u[:-1])
    u[1:] = .5*(u[1:]+up[1:] -  c*(up[1:]-up[:-1]))
    return u[1:-1] 


# Constants and parameters
a = 1.0 # wave speed
tmin, tmax = 0.0, 1.0 # start and stop time of simulation
xmin, xmax = 0.0, 2.0 # start and end of spatial domain
Nx = 80 # number of spatial points
c = 0.9 # courant number, need c<=1 for stability

8.7.2 Code example for various schemes for the advection equation

# src-ch6/advection_schemes.py

import numpy as np
from matplotlib import animation
from scipy import interpolate
from numpy import where
from math import sin

import matplotlib; matplotlib.use('Qt4Agg')
import matplotlib.pylab as plt
plt.get_current_fig_manager().window.raise_()


LNWDT=2; FNT=15
plt.rcParams['lines.linewidth'] = LNWDT; plt.rcParams['font.size'] = FNT


init_func=1   # Select stair case function (0) or sin^2 function (1)

# function defining the initial condition
if (init_func==0):
    def f(x):
        """Assigning a value of 1.0 for values less than 0.1"""
        f = np.zeros_like(x)
        f[np.where(x <= 0.1)] = 1.0
        return f
elif(init_func==1):
    def f(x):
        """A smooth sin^2 function between x_left and x_right"""
        f = np.zeros_like(x)
        x_left = 0.25
        x_right = 0.75
        xm = (x_right-x_left)/2.0
        f = where((x>x_left) & (x<x_right), np.sin(np.pi*(x-x_left)/(x_right-x_left))**4,f) 
        return f

def ftbs(u): # forward time backward space
    u[1:-1] = (1-c)*u[1:-1] + c*u[:-2]
    return u[1:-1]

# Lax-Wendroff
def lax_wendroff(u): 
    u[1:-1] = c/2.0*(1+c)*u[:-2] + (1-c**2)*u[1:-1] - c/2.0*(1-c)*u[2:]
    return u[1:-1]

# Lax-Friedrich Flux formulation
def lax_friedrich_Flux(u):
    u[1:-1] = (u[:-2] +u[2:])/2.0 -  dt*(F(u[2:])-F(u[:-2]))/(2.0*dx)
    return u[1:-1] 

# Lax-Friedrich Advection
def lax_friedrich(u):
    u[1:-1] = (u[:-2] +u[2:])/2.0 -  c*(u[2:] - u[:-2])/2.0
    return u[1:-1] 

# macCormack for advection quation
def macCormack(u):
    up = u.copy()
    up[:-1] = u[:-1] - c*(u[1:]-u[:-1])
    u[1:] = .5*(u[1:]+up[1:] -  c*(up[1:]-up[:-1]))
    return u[1:-1] 


# Constants and parameters
a = 1.0 # wave speed
tmin, tmax = 0.0, 1.0 # start and stop time of simulation
xmin, xmax = 0.0, 2.0 # start and end of spatial domain
Nx = 80 # number of spatial points
c = 0.9 # courant number, need c<=1 for stability


# Discretize
x = np.linspace(xmin, xmax, Nx+1) # discretization of space
dx = float((xmax-xmin)/Nx) # spatial step size
dt = c/a*dx # stable time step calculated from stability requirement
Nt = int((tmax-tmin)/dt) # number of time steps
time = np.linspace(tmin, tmax, Nt) # discretization of time

# solve from tmin to tmax

solvers = [ftbs,lax_wendroff,lax_friedrich,macCormack]
#solvers = [ftbs,lax_wendroff,macCormack]
#solvers = [ftbs,lax_wendroff]
#solvers = [ftbs]

u_solutions=np.zeros((len(solvers),len(time),len(x)))
uanalytical = np.zeros((len(time), len(x))) # holds the analytical solution


    
for k, solver in enumerate(solvers): # Solve for all solvers in list
    u = f(x)
    un = np.zeros((len(time), len(x))) # holds the numerical solution

    for i, t in enumerate(time[1:]):
        
        if k==0:
            uanalytical[i,:] = f(x-a*t) # compute analytical solution for this time step
            
        u_bc = interpolate.interp1d(x[-2:], u[-2:]) # interplate at right bndry
        
        u[1:-1] = solver(u[:]) # calculate numerical solution of interior
        u[-1] = u_bc(x[-1] - a*dt) # interpolate along a characteristic to find the boundary value
        
        un[i,:] = u[:] # storing the solution for plotting
    
    u_solutions[k,:,:] = un



### Animation 
 
# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure()
ax = plt.axes(xlim=(xmin,xmax), ylim=(np.min(un), np.max(un)*1.1))

lines=[]     # list for plot lines for solvers and analytical solutions
legends=[]   # list for legends for solvers and analytical solutions

for solver in solvers:
    line, = ax.plot([], [])
    lines.append(line)
    legends.append(solver.__name__)

line, = ax.plot([], []) #add extra plot line for analytical solution
lines.append(line)
legends.append('Analytical')

plt.xlabel('x-coordinate [-]')
plt.ylabel('Amplitude [-]')
plt.legend(legends, loc=3, frameon=False)
 
# initialization function: plot the background of each frame
def init():
    for line in lines:
        line.set_data([], [])
    return lines,

# animation function.  This is called sequentially
def animate(i):
    for k, line in enumerate(lines):
        if (k==0):
            line.set_data(x, un[i,:])
        else:
            line.set_data(x, uanalytical[i,:])
    return lines,

def animate_alt(i):
    for k, line in enumerate(lines):
        if (k==len(lines)-1):
            line.set_data(x, uanalytical[i,:])
        else:
            line.set_data(x, u_solutions[k,i,:])
    return lines,

 
# call the animator.  blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate_alt, init_func=init, frames=Nt, interval=100, blit=False)
 
 
plt.show()

Result from code example above using a step function as the initial value.

Result from code example above using a sine squared function as the initial value

• «
• 1
• ...
• 65
• 66
• 67
• 68
• 69
• 70
• 71
• 72
• 73
• 74
• 75
• 76
• 77
• 78
• 79
• »