Finite Difference Method Stability

Question

The diffusion equation is:

$ \frac{\partial T}{\partial t} = \alpha \left( \frac{\partial^2 T}{\partial x^2} \right) $

An explicit finite difference approach can be used to solve this, forward in time and central differences in space. Approximating the diffusion equation at a node i, yields,

$\frac{T_i^{n+1}-T_i^n}{\Delta t} = \alpha\frac{T_{i+1}^n-2T_i^n+T_{i-1}^n}{\Delta x^2}$

which gives

$T_i^{n+1} = T_i^n(1-\omega)+\omega(0.5T_{i+1}^n+0.5T_{i-1}^n)$

where $\omega = \frac{2\alpha \Delta t}{\Delta x^2}$

The book states that for stability condition, the coefficients of the right-hand side terms must be positive which implies that

$\Delta t \lt \frac{\Delta x^2}{2\alpha}$

My question is why should the coefficients of the right-hand side terms be positive for stability?

Answer 1

Analytical problem :

what you are expecting is positive diffusion : you want the $T_i$ values to spread over your domain as time passes to eventually reach $T_i(t\rightarrow \infty) = cte$

if $\alpha$ was a negative number, you would have what is called negative diffusion : you'll have exactly the opposite i.e. the gradients will get greater through time.

The sign of $\alpha$ hence dictates the behaviour of your analytical solution.

$\alpha < 0$ case :

Ideally, the numerical solution should have the same behaviour as the analytical solution. However, the finite difference theory assumes the solution to be smooth : if the solution features gradients that are too sharp, then your numerical method will not be able to handle them.

We have just said that in the case where $\alpha < 0$, the gradients grow greater with time. The error generated by the simulation will not be smeared out, as would be the case with positive diffusion $\alpha > 0$, but instead will be amplified. For that reason, if α<0, you know for sure your simulation is going to blow up at some point.

$\alpha > 0$ case :

If $\alpha > 0$, you are however not safe. If your time step is too large, your simulation will not be stable either. The stability condition $\Delta t < \frac{\Delta x ^2}{2 \alpha}$ indicates whether your numerical method has a chance of being stable or not. Note that it is a necessary condition for your numerical method to be stable, not a sufficient condition. Yet in practice, it turns out to be a very powerful tool.

Also, the mesh Fourier number for a diffusive term can be defined as $\alpha \frac{\Delta t}{\Delta x^2}$. In practice it is more convenient to write the stability condition in terms of the mesh Fourier $\alpha \frac{\Delta t}{\Delta x^2} < \frac{1}{2}$

This way you can see that the parameters of your simulation $\Delta t$, $\Delta x$ and $\alpha$ are all on the left hand side and $\frac{1}{2}$ is the critical value that must not be exceeded for the simulation to have a chance of being stable. In practice, the value for $\alpha$ is given by your problem and you will have chosen $\Delta x$ already. Hence, $\Delta t$ is the only parameter you can play with so that the stability condition on diffusion is observed.

The value of the critical mesh Fourier number depends on the space and time discretisation you have chosen. Some time integrators have broader stability regions than others, hence they will allow larger mesh Fourier numbers. Practically speaking, this means you'd be able to choose larger time steps while still having a stable numerical method.

To summarise :

• if $\alpha < 0$, you will have negative diffusion and your simulation will in any case not be stable.
• if $\alpha > 0$, your simulation might be stable... or it may not !
• the stability condition on diffusion (and the mesh Fourier number) helps you choose the time step $\Delta t$ for your numerical method to be stable.

I recommend you make a dummy simulation and play with the parameters to see what happens. No need to waste time into programming something : a spreadsheet software is enough for your particular case.

---

Edit: partial rewrite of my answer to make it clearer

Answer 2

For stability you need $\omega \leq 1$, as stated here. If your diffusivity $\alpha$ is positive (that should be based on a physical basis), $\omega$ should also be positive since $\Delta x$ and $\Delta t$ are, by definition, positive numbers.

Answer 3

This is one of the canonical examples used to explain stability theory for PDEs, so there are many good explanations in textbooks. My favorite is in Chapter 8 of LeVeque.

The other answers that have been given contain some incorrect terms and some information that is not directly relevant. Most importantly, they use the term CFL condition in an incorrect way. From LeVeque's text, p. 217, the CFL condition is:

A numerical method can be convergent only if its numerical domain of dependence contains the true domain of dependence of the PDE, at least in the limit as $\Delta t$ and $\Delta x$ go to zero.

In the context of the heat equations, this implies merely that $\Delta t \propto \Delta x^{1+\epsilon}$ for some $\epsilon>0$. Of course, this condition is not sufficient for stability; the CFL condition is only necessary.

Answer 4

A "hand-waving" argument may be possible to explain the need for a positive coefficient for the stability of the time integration scheme. I don't know what that is.

The same condition can be derived rigorously through von Neumann stability analysis. von Neumann Stability analysis. This link examines a few different schemes.