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

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?

-
Welcome to SciComp! The scope of your question fits this site very well. To get good answers, however, you should be more specific. Please indicate, what in particular you don't understand. Your code looks well structured and documented, but for answering your question, maybe a smaller code snippet will do. –  Jan May 26 '13 at 12:11
It might make your debugging simpler if you use a test case where your input vector has, say, 5 elements, and you step through the code with a debugger like gdb and ddd. Doing so might help narrow down the source of the error. I think that most code debugging questions don't work very well here because the majority of the work involves figuring out where the bug is in the first place. Once you find it, the explanation is frequently (but not always) straightforward. Can you run a unit test to figure out if there are any bugs in Tridiagonal that might be causing this behavior? –  Geoff Oxberry May 26 '13 at 18:37
Have a look at this example in the Wikipedia entry on the Crank-Nicolson method. If you set $k$ and $D_x$ to zero, it will turn into a plain advection problem. It remains to incorporate the boundary conditions... –  Jan May 28 '13 at 12:23

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

"""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()

-
Actually I saw your earlier question too but it was so general I couldn't answer (when you posted a page of code). In my experience, when you ask good question on this site the people are very helpful. Bon voyage! –  boyfarrell May 28 '13 at 13:57
I'm just kidding. –  lost_with_coding May 28 '13 at 13:58
Thank your for your well crafted answer! –  vanCompute May 28 '13 at 20:02
@boyfarrel Is there any chance you have a C++/C version of this. Its fine if no. I don't use matlab much and I don't feel like learning it. Even fortran would be better. –  lost_with_coding May 28 '13 at 20:36
<deleted inappropriate comments> –  Paul May 30 '13 at 3:49

Sorry guys, but the Crank-Nicolson method is totally inappropriate for an advection problem. Stability and accuracy of the local differential approximation unfortunately do not guarantee consistency. You will get the solution of a different problem. In some trivial cases you may get lucky but the advection equation is in general unforgiving. Check the numerical analysis textbooks. The Crank-Nicolson scheme is never used. People still have to adopt it in some engineering problems, because of the stability, but there they have to control the solution as it evolves and have heuristic to test error bounds. That is, they know there can be big errors, but because of the nature of the problem they can live with it. Interestingly, with a wave problem like Schroedinger equation, Crank-Nicholson works like a charm, but that is because the equation is a combination of two advection equations for real and imaginary parts which feed into each other at each time step.

-
Welcome to scicomp. I am afraid that people might take your answer serious. That's why I give a -1: This answer is rather a rant and does not provide any reliable evidence for what it claims to be the case. –  Jan Dec 17 '14 at 13:16