Note: this post was originally written by David Ketcheson.
In their celebrated 1928 paper, Courant, Friedrichs, and Lewy proved a geometric condition that must be satisfied by a convergent partial differential equation discretization – the famous CFL condition. Briefly, the CFL theorem says that the numerical method must transport information at least as quickly as information travels in the true PDE solution. The proof is geometric and is conveyed through numerous diagrams.
Exactly fifty years later, in 1978, Rolf Jeltsch and Olavi Nevanlinna published a theorem [JN] that deals with bounding the modulus of a polynomial \(\psi(z)\) over a disk of the form \[D_r = {z \in \{\mathbb C} : |z-r|\le r\}.\] Their theorem says that if \(\psi(z) = 1 + z + a_2 z^2 + \cdots + a_s z^s\) and \(|\psi(z)|\le 1\) for all \(z\) in such a disk \(D_r\), then the disk radius \(r\) is at most \(s\). The proof of this result is, of course, purely algebraic.
These results apparently have nothing to do with one another. And yet it turns out that they are equivalent statements! That is, the CFL theorem can be proved using the JN disk theorem. And the JN disk theorem can be proved using the CFL condition (and no algebraic techniques). This was explained in a beautiful paper of Sanz-Serna and Spijker [SS] in 1986, and the result deserves to be much more well known.
First order upwinding
Consider the problem of approximating the value \(u(x_i,t_n)\) for the advection equation \[u_t + u_x = 0.\] The exact solution can be obtained by characteristics from the previous time level: \[u(x_i,t_n) = u(x_i-k,t_{n-1}),\] where \(k\) is the time step size. The CFL theorem says that the stencil used for approximating \(u(x_i,t_n)\) must enclose the point \(x_i-k\).
Let’s discretize the advection equation in space using upwind differences: \[U_i'(t) = -\left(U_i-U_{i-1}\right).\] Here for simplicity we’ve assumed a spatial mesh width of 1. Taking periodic boundary conditions, this semi-discretization is a system of ODEs of the form \(U'=LU\) where \(L\) is the circulant matrix \[ \begin{pmatrix} -1 & & & 1 \\ 1 & -1 & & \\ & \ddots & \ddots \\ & & 1 & -1 \\ \end{pmatrix}\] (as usual, all the omitted entries are zero). The eigenvalues of this matrix all lie on the boundary of the disk of radius one centered at \(z=-1\), which we denote by \(D_1\). Here are the eigenvalues of a 50-point discretization:
If we discretize in time with Euler’s method, we get the scheme \[U^n_i = U^{n-1}_i - k\left(U_i-U_{i-1}\right).\] This scheme computes the solution at \((x_i,t_n)\) using values at \((x_{i-1},t_{n-1})\) and \((x_i,t_{n-1})\), so the CFL theorem says it can be convergent only if \(x_i-k\) lies in the interval \((x_{i-1},x_i)\). Since \(x_{i-1} = x_i - 1\), this holds iff the step size \(k\) is smaller than 1.
This result – that the first-order upwind method is stable and convergent only for CFL number at most one – is well known, and can also be derived using basic method of lines stability theory. The stability function for Euler’s method is \(\psi(z) = 1 + z\), so it is stable only if \(z=k\lambda\) lies in the disk \(\{z : |1+z|\le 1\} = D_1\) for each eigenvalue \(\lambda\) of \(L\). What we have seen in the foregoing is that this stability condition can be derived directly from the CFL condition, without considering the eigenvalues of \(L\) or the stability region of Euler’s method.
Proving the JN disk theorem via the CFL theorem
For higher order discretizations, the CFL condition is necessary but not generally sufficient for stability. Nevertheless, we can use it to derive the JN disk theorem. I’ll restrict the explanation here to Runge-Kutta methods, but the extension to multistep methods is very simple. Suppose that we discretize in time using a Runge-Kutta method with \(s\) stages. In each stage, one point further to the left is used, so typically the stencil for computing \(u(x_i,t_n)\) includes the values from the previous step at \(x_{i-s}, x_{i-s+1}, \dots, x_i\). Thus the CFL theorem says the method cannot be convergent unless \(x_i-k\) lies in the interval \((x_{i-s},x_i)\); i.e., unless \(k\le s\). Meanwhile, the stability function \(\psi(z)\) of the Runge-Kutta method is a polynomial of degree at most \(s\). Method of lines analysis tells us that the full discretization is stable if \(kD_1\) lies inside the region \(\{z : |\psi(z)|\le 1\}.\) Since we know it is unstable for \(k>s\), this implies that if \(|\psi(z)|\le 1\) over the disk \(D_k\), then \(k \le s\).
Recap
- An \(s\)-stage upwind discretization has stencil width \(s\).
- The CFL condition implies that this discretization cannot be convergent for Courant numbers larger than \(s\).
- The spectrum of the semi-discretization is the boundary of the disk \(D_1\).
- Stability analysis implies that the full discretization is convergent if the scaled spectrum \(kD_1 = D_k\) lies inside the stability region of the time discretization.
- Thus no \(s\)-stage time discretization can have a stability region including the disk larger than \(D_s\) (this is the content of the JN disk theorem).
Ellipses
Of course, we didn’t have to choose first-order upwinding in space; we could have taken any spatial discretization. For instance, if we use centered differences: \[U_i'(t) = -\left(U_{i+1}-U_{i-1}\right)\] then the spectrum of the semi-discretization lies on the imaginary axis in the interval \([-i,i]\). Then the same line of reasoning then tells us that the largest imaginary-axis interval of stability for an \(s\)-stage method is \([-is,is]\). By considering convex combinations of upwind and centered differences, we get similar results for a family of ellipses; this is the content of Theorem 5 of [SS].
Parabolic problems
It’s well known that the largest interval of stability of a consistent \(s\)-stage method on the negative real axis has length \(s^2\); the corresponding polynomials are (shifted) Chebyshev polynomials. You might hope that this could also be deduced by considering a centered difference semi-discretization of the heat equation and applying the CFL theorem. That would be very neat, since it would provide a connection between PDE stability theory and the optimality of Chebyshev polynomials.
Indeed, explicit time discretizations generally lead to step size restrictions depending on the square of the spatial mesh width when paired with the usual centered spatial discretization. But the CFL theorem is not sharp for these discretizations; it only tells us that \(k\) must vanish vanish more quickly than the spatial mesh width. So no deduction along these lines seems possible.
#spnetwork #recommend doi:/10.1007/BF01389633
#discusses doi:10.1147/rd.112.0215 #discusses doi:10.1007/BF01932030