8
$\begingroup$

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?

$\endgroup$
4
  • 1
    $\begingroup$ 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. $\endgroup$
    – Jan
    May 26 '13 at 12:11
  • $\begingroup$ 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? $\endgroup$ May 26 '13 at 18:37
  • $\begingroup$ 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... $\endgroup$
    – Jan
    May 28 '13 at 12:23
  • 1
    $\begingroup$ 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 erro $\endgroup$
    – user13061
    Dec 12 '14 at 16:18
20
+100
$\begingroup$

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.

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()
$\endgroup$
4
  • 3
    $\begingroup$ 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! $\endgroup$
    – boyfarrell
    May 28 '13 at 13:57
  • $\begingroup$ I'm just kidding. $\endgroup$ May 28 '13 at 13:58
  • $\begingroup$ @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. $\endgroup$ May 28 '13 at 20:36
  • 2
    $\begingroup$ <deleted inappropriate comments> $\endgroup$
    – Paul
    May 30 '13 at 3:49
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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