I am not very familiar with the common discretization schemes for PDEs. I know that Crank-Nicolson is popular scheme for discretizing the diffusion equation. Is also a good choice for the advection term?

I am interesting in solving the Reaction-Diffusion-Advection equation,

$\frac{\partial u}{\partial t} + \nabla \cdot \left( \boldsymbol{v} u - D\nabla u \right) = f$

where $D$ is the diffusion coefficient of substance $u$ and $\boldsymbol{v}$ is the velocity.

For my specific application the equation can be written in the form,

$\frac{\partial u}{\partial t} = \underbrace{D\frac{\partial^2 u}{\partial x^2}}_{\textrm{Diffusion}} + \underbrace{\boldsymbol{v}\frac{\partial u}{\partial x}}_{\textrm{Advection (convection)}} + \underbrace{f(x,t)}_{\textrm{Reaction}}$

Here is the Crank-Nicolson scheme I have applied,

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

Notice the $\alpha$ and the $\beta$ terms. This enables scheme to move between:

  • $\beta=\alpha=1/2$ Crank-Niscolson,
  • $\beta=\alpha=1$ it is fully implicit
  • $\beta=\alpha=0$ it is fully explicit

The values can be different, which allows the diffusion term to be Crank-Nicolson and the advection term to be something else. What is the most stable approach, what would you recommend?


This is a well-framed question and a very useful thing to understand. Korrok is correct to refer you to von Neumann analysis and LeVeque's book. I can add a bit more to that. I'd like to write a detailed answer, but at the moment I only have time for a short one:

With $\alpha=\beta=1/2$, you get a method that is absolutely stable for arbitrarily large step sizes, as well as second-order accurate. However, the method is not L-stable, so very high frequencies will not be damped, which is unphysical.

With $\alpha=\beta=1$, you get a method that is also unconditionally stable, but only 1st-order accurate. This method is very dissipative. It is L-stable.

If you take $\alpha\ne\beta$, your method can be understood as applying an additive Runge-Kutta method to the centered-difference semi-discretization. The stability and accuracy analysis for such methods is considerably more complicated. A very nice paper on such methods is here.

Which approach to recommend depends strongly on the magnitude of $D$, the kind of initial data you deal with, and the accuracy you seek. If very low accuracy is acceptable, then $\alpha=\beta=1$ is a very robust approach. If $D$ is moderate or large, then the problem is diffusion-dominated and very stiff; typically $\alpha=\beta=1/2$ will give good results. If $D$ is very small, then it may be advantageous to use an explicit method and higher-order upwinding for the convective terms.

  • $\begingroup$ A very insightful answers, thank you! Is there a way to define the different regimes of being diffusion dominated and advection dominated? Other than comparing the magnitude of the terms? For example, by comparing coefficients only? What is the meaning of the technical term L-stability. Everybody recommends this book, I must buy it! $\endgroup$ – boyfarrell Mar 3 '13 at 23:52
  • $\begingroup$ The criterion I gave you only involves the coefficients. Briefly, L-stability means that high frequencies will be strongly damped. $\endgroup$ – David Ketcheson Mar 4 '13 at 0:41
  • $\begingroup$ So when $u(x)$ is a smooth function (as in the sense, it doesn't have high frequency Fourier components) Crank-Nicolson is a good choice. If however, $u(x)$ has sharp edges then $\beta=1$ is a good choice. $\endgroup$ – boyfarrell Mar 4 '13 at 14:35
  • $\begingroup$ That's a reasonable, though very rough, generalization. Those choices will at least work if you don't need very much accuracy. $\endgroup$ – David Ketcheson Mar 4 '13 at 17:26
  • $\begingroup$ Question: is a method a Crank-Nicholson method only if spatial derivatives (first and second) are centered? Or can one use any spatial derivative? $\endgroup$ – Millemila Apr 15 at 19:40

Generally speaking, you'll want to use an implicit method for parabolic equations (the diffusion part) -- explicit schemes for parabolic PDE need to have a very short timestep to be stable. Conversely, for the hyperbolic part (advection) you'll want an explicit method as it's cheaper and doesn't disrupt the symmetry of the linear system you have to solve by using an implicit scheme for diffusion. In that case, you want to avoid centered differences like $(u_{j+1}-u_{j-1})/2\Delta t$ and switch to one-sided differences $(u_j-u_{j-1})/\Delta t$ for reasons of stability.

I'd suggest you look at Randy Leveque's book or Dale Durran's book for "von Neumann stability analysis". It's a general approach to ascertaining the stability of your discretization scheme, provided you have periodic boundary conditions. (There's also a good wiki article here.)

The basic idea is to assume that your discrete approximation can be written a sum of plane waves $e^{i(k j\Delta x-\omega n\Delta t)}$, where $k$ is the wave number and $\omega$ the frequency. You cram a plane wave into your approximation to the PDE and pray it doesn't blow up. We can rewrite the plane wave as $\xi^n e^{ikj\Delta x}$ and we want to make sure that $|\xi|\le 1$.

By way of illustration, consider the ordinary diffusion equation with fully implicit differencing:

$\frac{u^{n+1}_j-u^n_j}{\Delta t} = D\frac{u^{n+1}_{j-1}-2u^{n+1}_j+u^{n+1}_{j+1}}{\Delta x^2}$

If we substitute in a plane wave, then divide by $\xi^n$ and $e^{ikj\Delta x}$, we get the equation

$ \frac{\xi-1}{\Delta t} = D\frac{e^{-ik\Delta x}-2+e^{ik\Delta x}}{\Delta x^2}\xi $

Clean this up a little bit now and we get:

$\xi = \frac{1}{1+\frac{2D\Delta t}{\Delta x^2}\left(1-\cos k\Delta x\right)}$.

This is always less than one, so you're in the clear. Try applying this for the explicit, centered scheme for the advection equation:

$ \frac{u_j^{n+1}-u_j^n}{\Delta t} = v\frac{u_{j-1}^n-u_{j+1}^n}{2\Delta x}$

and see what $\xi$ you get. (It'll have an imaginary part this time.) You'll find that $|\xi|^2 > 1$, which is sad times. Hence my admonition that you don't use it. If you can do that, then you shouldn't have much trouble finding a stable scheme for the full advection-diffusion equation.

That said, I'd use a fully implicit scheme for the diffusion part. Change the differencing in the advection part to $u_j-u_{j-1}$ if $v > 0$ and $u_j-u_{j+1}$ if $v < 0$ and pick a timestep so that $V\Delta t/\Delta x \le 1$. (This is the Courant-Friedrichs-Lewy condition.) It's only first-order accurate, so you may want to look up higher order discretization schemes if that concerns you.

  • $\begingroup$ That's a really detailed answer, thank you. $\endgroup$ – boyfarrell Mar 2 '13 at 5:32
  • $\begingroup$ This answer considers only discretizations based on the forward and backward Euler methods in time. The question is about Crank-Nicholson. $\endgroup$ – David Ketcheson Mar 3 '13 at 3:56

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.