Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

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?

share|improve this question
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
add comment

1 Answer

up vote 14 down vote accepted

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} $$


$$-\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.

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)

share|improve this answer
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
add comment

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.