How to discretize the advection equation using the Crank-Nicolson method?

Question

The advection equation needs to be discretized in order to be used for the Crank-Nicolson method. Can someone show me how to do that?

Answer 1

Starting with the advection equation is conservative form,

$$ \frac{\partial u}{\partial t} = -\frac{\partial (\boldsymbol{v} u)}{\partial x} + s(x,t) $$

The Crank-Nicolson method consists of a time averaged centered difference.

$$\frac{u_{j}^{n+1} - u_{j}^{n}}{\Delta t} = -\boldsymbol{v} \left[ \frac{1-\beta}{2\Delta x} \left( u_{j+1}^{n} - u_{j-1}^{n} \right) + \frac{\beta}{2\Delta x} \left( u_{j+1}^{n+1} - u_{j-1}^{n+1} \right) \right] + s(x,t)$$

Regarding the notation, subscripts are for points in space, and the superscripts are for points in time.

The points at $n+1$ are in the future: they are unknowns. We now need to rearrange the above equation so that all knowns are on the r.h.s and unknowns are on the l.h.s..

Making the substitution,

$$ r = \frac{\boldsymbol{v}}{2}\frac{\Delta t}{\Delta x} $$

gives,

$$-\beta r\phi_{j-1}^{n+1} + \phi_{j}^{n+1} + \beta r\phi_{j+1}^{n+1} = (1-\beta)r\phi_{j-1}^{n} + \phi_{j}^{n} - (1-\beta)r\phi_{j+1}^{n}$$

This is the advection equation discretised using the Crank-Nicolson method. You can write it as a matrix equation,

$$ \begin{pmatrix} 1 & \beta r & & & 0 \\ -\beta r & 1 & \beta r & & \\ & \ddots & \ddots & \ddots & \\ & & -\beta r & 1 & \beta r \\ 0 & & & -\beta r & 1 \\ \end{pmatrix} \begin{pmatrix} u_1^{n+1} \\ u_2^{n+1} \\ \vdots \\ u_{J-1}^{n+1} \\ u_{J}^{n+1} \\ \end{pmatrix} = \begin{pmatrix} 1 & -(1 - \beta)r & & & 0 \\ (1 - \beta)r & 1 & -(1 - \beta)r & & \\ & \ddots & \ddots & \ddots & \\ & & (1 - \beta)r & 1 & -(1 - \beta)r \\ 0 & & &(1 - \beta)r & 1 \\ \end{pmatrix} \begin{pmatrix} u_1^{n} \\ u_2^{n} \\ \vdots \\ u_{J-1}^{n} \\ u_{J}^{n} \\ \end{pmatrix} $$ Setting $\beta=1/2$ will give you trapezoidal integration in time, so for Crank-Nicolson this is what you want.

A few words of warning. This is basic solution you wanted, but you will need to include some sort of boundary condition for a well-posed problem. Also, Crank-Nicolson is not necessarily the best method for the advection equation. It is second order accurate and unconditionally stable, which is fantastic. However it will generate (as with all centered difference stencils) spurious oscillation if you have very sharp peaked solutions or initial conditions.

I wrote the following code for you in Python, it should get you started. The code solves the advection equation for an initial Gaussian curve moving to the right with constant velocity.

from __future__ import division
from scipy.sparse import spdiags
from scipy.sparse.linalg import spsolve
import numpy as np
import pylab

def make_advection_matrices(z, r):
    """Return matrices A and M for advection equations"""
    ones = np.ones(len(z))
    A = spdiags( [-beta*r, ones, beta*r], (-1,0,1), len(z), len(z) )
    M = spdiags( [(1-beta) * r, ones, -(1-beta) * r], (-1,0,1), len(z), len(z) )
    return A.tocsr(), M.tocsr()

def plot_iteration(z, u, iteration):
    """Plot the solver progress"""
    pylab.plot(z, u, label="Iteration %d" % iteration)

# Set up basic constants
beta = 0.5
J = 200 # total number of mesh points
z = np.linspace(-10,10,J) # vertices
dz = abs(z[1]-z[0]) # space step
dt = 0.2    # time step
v = 2 * np.ones(len(z)) # velocity field (constant)
r = v / 2 * dt / dz

# Initial conditions (peak function)
gaussian = lambda z, height, position, hwhm: height * np.exp(-np.log(2) * ((z - position)/hwhm)**2)
u_init = gaussian(z, 1, -3, 2)

A, M = make_advection_matrices(z, r)
u = u_init
for i in range(10):
    u = spsolve(A, M * u)
    plot_iteration(z, u, i)

pylab.legend()
pylab.show()

Answer 2

User 03161 asserts that the Crank Nicolson method is not appropriate for advection problems, but boyfarrell provides a working code with results visualized in a movie. In fact they are both correct, but neither gives the full perspective.

At the beginning of boyfarrell's answer the correct C-N formula for linear advection is written out. In that formula are two important dimensionless numbers. One is explicitly pointed out; it is $\beta$, the proportion of the spatial derivative that is estimated from the advanced time level. The other is denoted by $r$, which is half of what, in a more usual notation, is known as the Courant number, $\nu=\frac{v\Delta t }{\Delta x}$. This is the ratio between the true propagation speed $v$ and the speed $\Delta x/\Delta t $ with which information propagates through the grid. For explicit methods taking the spatial derivative only from the early time, the Courant number must be less than 1.0, for reasons of causality.

When the later time is taken into account also, there is no obviously necessary limit on $\nu$, which is why people are willing to go to such methods, that are called implicit and are very costly, especially for non-linear problems in more than one dimension. Nevertheless, there are theoretical attractions.

The correct outcome from a Fourier analysis is that if $\beta$ is chosen greater than 0.5, there really is no limit on the time step. If beta is chosen exactly equal to 0.5 we find neutral stability, which should be ideal. Unfortunately, waveforms get very distorted and large oscillations appear. The results are only acceptable in practice if $\beta$ is made greater than 0.5, unless the data is very smooth. (And in that case you are paying more than you have to)

Use of the C-N method has been very rare (in any application that I know of) for the past fifty years. However, it still makes an instructive homework problem. But there are a lot of much better advection schemes.