# Finite Difference Method Stability

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?

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

• It is not correct to refer to CFL numbers in this context -- the CFL number refers to a parameter in discretization of hyperbolic PDEs. The use here of the term "CFL condition" is also not correct. Aug 25, 2014 at 18:38
• Also, if $\alpha<0$ the problem itself is unstable -- this has nothing to do with the question. Aug 25, 2014 at 18:40
• @DavidKetcheson Here is what I suggest : I will now explain my view, please read it and tell me exactly on what aspects you disagree with me. Aug 28, 2014 at 7:15
• Point #1 : If I understand you correctly, you are implying that CFL number only exists for advective terms, is that correct ? This is perhaps the original definition of the CFL number. However, in the litterature, it is not uncommon to see it's definition broadened to cover diffusive terms. To me, it is a fair thing to do since in either cases the idea is to have a coefficient indicating by which factor the $\sum a_k q_k$ of a term will have an effect on the time evolution of the numerical solution. So why not call it a "CFL number for advection" or "CFL number for diffusion" ? Aug 28, 2014 at 7:16
• Point #2 : If $\alpha < 0$, it is true that the solution is unbounded. In my mind, the analytical solution is what it is, and cannot be considered either stable or unstable. The numerical method however may or may not be compliant with the analytical solution. If not, then the numerical solution is said not to converge i.e. it is either non-consistent or unstable or both. Aug 28, 2014 at 7:17

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.

• @DavidKetcheson, I am by no-means an expert in this topic but according to the references below that criterion is right: * Burden & Faires. Numerical Analysis (2010). Section 12.2, p 729 $$\alpha^2 \frac{k}{h^2} \leq \frac{1}{2}$$ * Leveque. Finite Difference Methods for Ordinary and Partial Differential Equations: Steady State and Time Dependent Problems (2007). Section 9.3, p 186 $$\frac{k}{h^2} \leq \frac{1}{2}$$. Aug 25, 2014 at 21:32
• The criterion is right. Calling it a CFL condition is wrong. The actual CFL condition is the one I gave, but (as in many cases) it is not sufficient for convergence. Aug 26, 2014 at 12:50
• Thank your for your clarification. Call that dimensionless number Courant number is also wrong? Aug 26, 2014 at 14:20
• Yes. So if you just delete the last sentence, your answer is great. Aug 28, 2014 at 9:51

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.

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.