You are given the following linear diffusion advection equat

You are given the following linear diffusion advection equation to solve with c, D as constants. Write down the Crank Nicholson discretization for this PDE with parameter beta. Make sure that your discretization has an error of no more than O(delta t, delta x^2) Use von Neumann analysis to prove that this scheme is unconditionally stable.

Solution

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 (uj+1uj1)/2t(uj+1uj1)/2t and switch to one-sided differences (ujuj1)/t(ujuj1)/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 ei(kjxnt)ei(kjxnt), where kk is the wave number and 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 neikjxneikjxand we want to make sure that ||1||1.

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

un+1junjt=Dun+1j12un+1j+un+1j+1x2ujn+1ujnt=Duj1n+12ujn+1+uj+1n+1x2

If we substitute in a plane wave, then divide by nn and eikjxeikjx, we get the equation

1t=Deikx2+eikxx21t=Deikx2+eikxx2

Clean this up a little bit now and we get:

=11+2Dtx2(1coskx)=11+2Dtx2(1coskx).

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

un+1junjt=vunj1unj+12xujn+1ujnt=vuj1nuj+1n2x

and see what you get. (It\'ll have an imaginary part this time.) You\'ll find that ||2>1||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 ujuj1ujuj1 if v>0v>0 and ujuj+1ujuj+1 if v<0v<0 and pick a timestep so that Vt/x1Vt/x1. (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.