The Crank-Nicholson scheme

The Crank-Nicholson implicit scheme for solving the diffusion equation (see Sect. 6.6) can be adapted to solve the advection equation. Thus, taking the average of the right-hand side of Eq. (234) between the beginning and end of the time-step, we obtain the differencing scheme written below:

$$u_i^{n+1} + \frac{C}{4}\,(u_{i+1}^{n+1}-u_{i-1}^{n+1}) = u_i^n - \frac{C}{4}\,(u_{i+1}^{n}-u_{i-1}^{n}).$$ (246)

A von Neumann stability analysis of the above scheme yields the following expression for the amplification factor:

$$A = \frac{1 - {\rm i} \,(C/2)\,\sin (k \,\delta x)}{1 + {\rm i} \,(C/2)\,\sin (k \,\delta x)}.$$ (247)

Note that $\vert A\vert = 1$ for all values of $k$, irrespective of the value of $C$. This implies that the Crank-Nicholson implicit scheme is not subject to the CFL constraint, $C<1$ (since there is no value of $k$ for which $\vert A\vert>1$). Moreover, there is no spurious decay in the Fourier harmonics of the solution (since $\vert A\vert = 1$). Hence, unlike the Lax scheme, we would not expect the Crank-Nicholson scheme to introduce strong numerical dispersion into the advection problem.

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

// Advect1D.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 Crank-Nicholson scheme.

#include <blitz/array.h>

using namespace blitz;

void Tridiagonal (Array<double,1> a, Array<double,1> b, Array<double,1> c, 
                  Array<double,1> w, Array<double,1>& u);

void Advect1D (Array<double,1>& u, double C)
{  
  // Find N. Declare local arrays.
  int N = u.extent(0) - 2;
  Array<double,1> a(N+2), b(N+2), c(N+2), w(N+2);

  // Initialize tridiagonal matrix
  for (int i = 2; i <= N;   i++) a(i) = - 0.25 * C;
  for (int i = 1; i <= N; i++) b(i) = 1.;
  for (int i = 1; i <= N-1; i++) c(i) = + 0.25 * C;

  // Initialize right-hand side vector
  for (int i = 1; i <= N; i++)
    w(i) = u(i) - 0.25 * C * (u(i+1) - u(i-1));

  // Invert tridiagonal matrix equation
  Tridiagonal (a, b, c, w, u);

  // Calculate i=0 and i=N+1 values
  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 as specified in Eq. (244), and the calculation is performed with $v=1$, $\delta t = 9.95\times 10^{-3}$, and $N=200$. Note that $C=2.0$ for these parameters. It can be seen that the pulse propagates at (almost) the correct speed, and maintains approximately the same shape. Clearly, the performance of the Crank-Nicholson scheme is vastly superior to that of the Lax scheme.

Figure 77: Advection of a 1-d Gaussian pulse. Numerical calculation performed using $v=1$, $\delta t = 9.95\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$.