von Neumann stability analysis is a great tool to assess the stabilty of discretizing schemes for PDEs. But too often, imho, the discussion is too convoluted. Here, I try to provide a shortcut.

In a nutshell:

• The standard approach involves Fourier techniques, involving (of course) complex numbers
• The real part of these numbers is analysed, with some trigonometric expression resulting, identifying the troublesome modes
• I claim this mode can be identified in advance, which makes the whole Fourier procedure unnecessary

BTW: it’s pronounced “fon no ee man”

Standard procedure

Starting from the diffusion equation (aka heat equation):

A centered space and forward time discretization throws a explicit scheme:

where

is a diffusion Courant number. The question is what value of this number is acceptable for the simulation to remain stable. In more practical terms, for a given spacial resolution we are asking what time step is acceptable. The shorter, the more accurate (in principle, but for stuff such as roundoff errors), but of course we’d rather have our results in minutes rather than hours or days.

Now, the usual discussion centers around the growth of the error in the solution. Let’s just suppose that the solution has a sine wave shape:

i.e.

Then the discretized equation (1) becomes

Now, the usual discussion relates the part with the complex exponentials to trigonometrical functions:

Now this means

Any law of the form

means that:

• the mode will decay monotonically if real( C ) < 1 ( C is a complex number!) (well, not in this case, but it might be, in general)
• the mode will decay in an oscillatory manner if real( C ) > – 1. In other words, the scheme is stable if | C | < 1
• If | C | > 1 the scheme will fail, because that mode will grow and grow, which it must not in principle (see below for Cahn Hilliard)
• A possible imaginary part is not so important, it being a change of phase of the mode

This applies to every mode (i.e. value of k ), but our main offender is seen to be the one whose real( C ) is the lowest: the one for which the argument of the squared sin is π/2 (for the time being, we’ll pretend this does not mean anything special). For this mode, C = 1- 4 α, and:

Hence, we finally arrive at our result:

• α must be below 1/4 for a monotonic decrease of the worst modes, it suffices it be below 1/2 for stability

Our approach

Now, seriously, the worst mode is the one for which

The wavelength of this mode is λ = 2 Δx. This is precisely the shortest wavelength mode that is allowed by our spacial resolution (that’s about everything we need to know of Fourier analysis).

This may look like a cosine function centered at j and with this wavelength takes vales +1 at j, -1 at j+1, 1 at j+2, et caetera.  It turns out, then, we did not need such a detour. Back to Eq. (1),

It seems clear all along that we would get the worst possible case if the values of the field Φ are staggered, producing a -1 -2 -1= -4 in the parenthesis. Then:

and in this case

$latex  A^{n + 1} = \A^{n} [ 1 – 4 \alpha ] .$

Which recovers the result above in two lines!

This may be extended very simply to other schemes. For example, an implicit scheme ends up in something like this (for the worst mode):

So, C=1/(1 + 4α ), which is always between 1 and 0. The method is then unconditionally stable.

A Crank-Nicolson method yields

So, for this scheme, α must be above 1.

Cahn-Hilliard equation

Let’s apply this to the more complicated Cahn-Hilliard equation

We’ll suppose that γ = 1 / 2 (it can be shown that this happens if the spacial length is set equal to the interfacial length at equilibrium).

Now, a tricky thing about this equation is that it describes segregation, and domain formation. So, the φ = 0 field is actually unstable against fluctuations. The systems should begin to form structures with values of  φ departing from 0, with a well defined wave-length (see Note below). [That’s the featured image, by the way. It started at values of 0.01, and see what it has done in just about 500 reduced time units.]

It makes little sense to use our stability analysis in this case, then.

We have to go close to the equilibrium regime, when φ ~ 1 (there is also φ ~ -1, what follows applies in the same way).

Let us write

The small field is now

and our equation can be approximated by

Now, a standard explicit scheme would yield:

(There’s a handy finite differences calculator for higher order derivatives). In this expression, α is as above, and

Now, the worst case is exactly as in diffusion, an in this case,

Hence, the scheme is stable if 1 – 8 α – 8 β > -1, i.e.  α + β < 1/4.

An implicit scheme leads to

which is unconditionally stable.

A mixed scheme such as the one that I am forced to work with (long story…) produces

hence the criterium for stability is

This means a typical, easy choice of Δt = Δx , in which α = β = 1, would be stable.

Actually (… but this gets progresively less interesting for the general public), my approach involves an explicit cube term, in which we may write, close to saturation

which leads to

hence the criterium for stability is

This means the choice Δt = Δx , in which α = β = 1, would not be stable!

Note on Cahn-Hilliard coarsening

If φ << 1, the cubed term is negibible. Then

Which is linear. Going to Fourier space,

Hence, all modes between k = 0 and √2 will grow, with the one with k=1 growing the most. This is a wavelength of 2 π , in units of the interfacial length.