Numerical solutions to diffusion equations

[Signifier]

Hello, the work I am doing now requires numerical solution of several reaction-diffusion equations such as the standard advection diffusion equation (ADE) or IE Fisher's equation. So far I have been using the only method I am familiar with, the simplest one I believe: to express the time derivative as a simple Euler finite difference approximation, and to express the Laplacian using the following approximation:

\frac {d^2C} {dx^2} = \frac{C(x + 1) - 2C(x) + C(x - 1)}{dx^2}}

where C is the density function (the variable diffusing), which is here a 1d matrix and dx is the step size. I approximated dC/dx with a backward difference approximation, and I approximated \frac {d^2C} {dx^2} with a forward difference approximation, then combined these to get the approximation above.

These techniques have worked for what I'm doing fairly well until now, but I would like to try using more accurate numerical methods. I am wondering what a simple, but higher order or more accurate method would be. Can anyone give me a name of a method? Can anyone recommend any good numerical methods books that I might be able to get through my university library or an interlibrary loan?

Also for a more accurate method would I have to keep track for more than just the state of the 1d matrix at the time before? IE would I have to keep track of the last 2 iterations, etc.

Thank you.

 

[jim mcnamara]

This is not my area. Apologies of this is of no help. That said, I once had a complimentary copy, according to my little card, of:

Numerical Methods for Advection-Diffusion Problems:
Notes on Numerical Fluid Mechanics
C. B. Vreugdenhil and B. Koren 1993

I do not how well that meets your needs, but university library will have a copy.

 

[nmtim]

Strikwerda's "Finite Difference Schemes and Partial Differential Equations" is a general introduction to finite difference methods, with a lot of coverage of convergence, stability, etc.

 

[Beeza]

What is the full equation?

 

[Dr Transport]

4-th order Runge-Kutta would be able to work for you in this case. You would need to write the 2nd order equaiton as a coupled 1st order set of equations. Try www.nr.com for the code and google it for the general method.

 

[AlephZero]

I'm curious as to how you would use R-K for a partial differential equation - I haven't seen that before.

A standard method for the diffusion equation is the Crank-Nicholson method. Basically it's the OP's discretization in space, but a central difference in time, not an Euler forward difference in time.

See http://farside.ph.utexas.edu/teaching/329/lectures/node80.html (and lots more Google hits).

C-N is unconditionally stable and 2nd order accurate, but if the time steps are large the solution has numerical oscillations in time between successive steps.

You can kill the oscillations (and reduce the accuracy) by using a backwards difference in time, or something in between. A general difference scheme is

T_i^{n+1} - T_i^n = \theta C\left(T_{i-1}^{n+1}-2\,T_i^{n+1}+T_{i+1}^{n+1}\right) + (1-\theta)C\left(T_{i-1}^n-2\,T_i^n+T_{i+1}^n\right).

with \theta = 1/2 = 1/2 for C-N,

\theta = 1 for backward difference (good for severely nonlinear problems)

\theta = about 0.6 or 2/3 is a reasonable compromise between accuracy and damping the oscillations - can work well for mildly nonlnear problems.

\theta = 0 (Euler method in time) is neither accurate nor stable. Its only "advantage" is it is an explicit method, i.e. you don't need to solve a set of equations at each time step. But for 1 space dimension, the equations are tridiagonal, therefore quick and easy to solve, so that "advantage" isn't worth much.

For 2 or 3 space dimensions, there are more efficient solution methods - e.g. alternating-direction implicit (ADI) methods, methods based on Fourier transforms, etc. But for one space dimension, C-N is pretty good.