41  Problem class on PDEs and the CFL condition

During the introduction to numerical methods on PDEs, it was mentioned that the first-order explicit Euler scheme is not well-suited for certain problems, such as the heat equation. In this chapter, we provide a brief note on deriving the famous CFL stability condition, which informs us of the timestep size required for a stable solution.

The following analysis is known as von Neumann analysis; it was originally derived in 1947 by Crank and Nicolson, and then later made rigorous by Charney et al. in 1950 (of which the third author was John von Neumann). You can read about the procedure here.

Consider the finite-difference formulation of the heat equation, as applied to a regularly spaced grid with spatial discretisation size \(\Delta x\) and temporal discretisation size \(\Delta t\): \[ \frac{u^{n+1}_j - u^n_j}{\Delta t} = D \frac{u^n_{j+1} - 2u^n_j + u^n_{j-1}}{\Delta x^2}. \] For ease of symbol manipulation, we have replaced the usual thermal diffusivity parameter with \(D\). Manipulation then gives \[ u^{n+1}_j = u^n_j + \frac{Dk}{h^2} (u^n_{j+1} - 2u^n_j + u^n_{j-1}), \] where we have replaced \(\Delta t\) with \(k\) and \(\Delta x\) with \(h\).

The basic idea of von Neumann analysis is to examine the growth rate when the solution is assumed to consist of sinusoidals. Let us assume that the solution at a mesh point is given by \[ u_k^n = G^n e^{i k x_j} = G^n e^{i\xi jh} \] where for simplicity we have assumed that the spatial points are given by \(x_j = jh\). The above is related to a sinusoidal wave with a wavenumber of \(\xi\) (assumed unspecified for now). The theory of Fourier series tells us that many functions can be considered as sums of the above sinusoidal waves, over many different values of \(\xi\). So our task is to understand the effect of one such wave, and to examine the growth rate \(G\).

Substitution into the finite-difference equation gives \[ G^{n+1} e^{i\xi jh} = G^n e^{i\xi jh} + \frac{Dk}{h^2} (G^n e^{i\xi (j+1)h} - 2G^n e^{i\xi jh} + G^n e^{i\xi (j-1)h}), \] or alternatively \[ \frac{G^{n+1}}{G^n} = 1 + \frac{Dk}{h^2} (e^{i\xi h} - 2 + e^{-i\xi h}). \] However, we can use the fact that \[ \cos(\xi h) = \frac{e^{i\xi h} + e^{-i\xi h}}{2}, \] and therefore \[ \frac{G^{n+1}}{G^n} = 1 + 2\frac{Dk}{h^2} (\cos(\xi h) - 1). \] Moreover, using another trigonometric identity, we can write \[ \frac{G^{n+1}}{G^n} = 1 - 4\frac{Dk}{h^2} \sin^2(\xi h/2) \equiv g(\xi, h, k). \] This is essentially a typical recurrence relation for \(G^n\) with a growth factor \(g(\xi, h, k)\). In order for the scheme to stable over all possible wavenumbers \(\xi\), we must then require that \(|g| < 1\) (otherwise \(G^n \to \infty\) as \(n \to \infty\)). Thus we require that \[ |1 - 4\frac{Dk}{h^2}| < 1. \] or \[ k < \frac{h^2}{2D} \Longrightarrow \Delta t < \frac{\Delta x^2}{2D}. \]