The Lax scheme

The instability in the differencing scheme (237) can be fixed by replacing $u_i^n$ on the right-hand side by the spatial average of $u$ taken over the neighbouring grid points. Thus, we obtain

$$u_i^{n+1} = \frac{1}{2}\,(u_{i+1}^n+u_{i-1}^n) -\frac{C}{2}\,(u_{i+1}^n-u_{i-1}^n),$$ (240)

which is known as the Lax scheme. A von Neumann stability analysis of the Lax scheme yields the following expression for the amplification factor:

$$A = \cos (k\,\delta x) - {\rm i} \,C\,\sin (k \,\delta x).$$ (241)

Now

$$\vert A\vert^2 = 1 - (1-C^2)\,\sin^2(k\,\delta x).$$ (242)

It follows that the Lax scheme is unconditionally stable (i.e., $\vert A\vert<1$ for all $k$), provided that $C<1$. From the definition of $C$, the inequality $C<1$ can also be written

$$\delta t < \frac{\delta x}{v}.$$ (243)

This is the famous Courant-Friedrichs-Lewy (or CFL) stability criterion. In fact, all stable explicit differencing schemes for solving the advection equation are subject to the CFL constraint, which determines the maximum allowable time-step.

Listed below is a routine which solves the 1-d advection equation via the Lax method.

// Lax1D.cpp

// Function to evolve advection equation in 1-d:

//  du / dt + v du / dx = 0  for  xl <= x <= xh

//  u = 0  at x=xl and x=xh

// Array u assumed to be of extent N+2.

// Now, ith element of array corresponds to

//  x_i = xl + i * dx    i=0,N+1

// Here, dx = (xh - xl) / (N+1) is grid spacing.

// Function evolves u by single time-step.

// C = v dt / dx, where dt is time-step.

// Uses Lax scheme.

#include <blitz/array.h>

using namespace blitz;

void Lax1D (Array<double,1>& u, double C)
{
  // Set N. Declare local array.
  int N = u.extent(0) - 2;
  Array<double,1> u0(N+2);

  // Evolve u
  u0 = u;
  for (int i = 1; i <= N; i++)
    u(i) = 0.5 * (u0(i+1) + u0(i-1)) - 0.5 * C * (u0(i+1) - u0(i-1));

  // Set boundary conditions
  u(0) = 0.;
  u(N+1) = 0.;
}

Figure 75 shows an example calculation which uses the above routine to advect a Gaussian pulse. The initial condition is

$$u(x,0) = \exp[-100\,(x-0.5)^2],$$ (244)

and the calculation is performed with $v=1$, $\delta t = 2.49\times 10^{-3}$, and $N=200$. Furthermore, $x_l = v\,t$ and $x_h = 1 + v\,t$. Note that $C=0.5$ for these parameters. It can be seen that the pulse is advected at the correct speed: i.e., the pulse appears approximately stationary when plotted versus $x-v\,t$. Unfortunately, the pulse does not remain the same shape (as it should). Instead, the pulse becomes gradually lower and wider as it propagates, and eventually diffuses away entirely.

Figure 75: Advection of a 1-d Gaussian pulse. Numerical calculation performed using $v=1$, $\delta t = 2.49\times 10^{-3}$, and $N=200$. The solid curve shows the initial condition at $t=0$, the short-dashed curve the numerical solution at $t=1$, the long-dashed curve the numerical solution at $t=2$, and the dot-dashed curve the numerical solution at $t=3$.

It is clear, from the above calculation, that the Lax scheme introduces a spurious dispersion effect into the advection problem. We can understand the origin of this effect by attempting a Fourier solution, $u(x,t)=\hat{u}_k(t)\,\exp(\,{\rm i}\,k\,x)$, of Eq. (234). We easily obtain

$$\hat{u}_k(t) = \hat{u}_k(0)\,\exp(-{\rm i}\,k\,v\,t).$$ (245)

Note that $\vert\hat{u}_k\vert$ is constant in time for all values of $k$. In other words, the amplitudes of the Fourier harmonics of a true solution of the advection equation remain constant in time--it is the phases of the harmonics which evolve. Let us now examine Eq. (242). It can be seen that, provided the CFL condition $C<1$ is satisfied, the magnitude of the amplification factor, $\vert A\vert$, is less than unity for all Fourier harmonics. In other words, the Lax differencing scheme causes the Fourier harmonics to decay in time. It is this unphysical attenuation of the Fourier harmonics which gives rise to the strong dispersion effect illustrated in Fig. 75.

Figure 76: Advection of a 1-d Gaussian pulse. Numerical calculation performed using $v=1$, $\delta t = 9.95\times 10^{-3}$, and $N=200$. The dotted curve shows the initial condition at $t=0.00$, the short-dashed curve the numerical solution at $t=0.15$, the long-dashed curve the numerical solution at $t=0.30$, and the solid curve the numerical solution at $t=0.37$.

Figure 76 shows a calculation made using the Lax scheme in which the CFL condition is violated. This calculation is identical to the one discussed previously, except that the time-step has been increased to $\delta t = 9.95\times 10^{-3}$, yielding a CFL parameter, $C=2.0$, which exceeds unity. It can be seen that the pulse grows in amplitude, and eventually starts to break up due to a short-wavelength instability.