Is Crank-Nicolson a stable discretization scheme for Reaction-Diffusion-Advection (convection) equation?

Question

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?

Answer 1

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.

Answer 2

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.