Understanding the Courant–Friedrichs–Lewy condition

Question

I understand these equations in particular can be solved easily without use of computational methods. Although right now I am concerned with trying to solve these equations using numerical integration strictly.

$$ \frac {\partial ^2 E}{\partial t ^2} = c^2 \frac {\partial ^2 E}{\partial z ^2} \tag{1}$$

$$E(t=0) = \sin(z), \frac {\partial E}{\partial t}(t=0) = -c \cos(z)$$

where c $= 3 \times 10^8 \frac{\text{m}}{\text{s}}$

Using the method of lines, I am trying to reduce this to a system of first-order differential equations and then apply a 4th-order Runge-Kutta method.

Although a necessary condition is the CFL condition which states (for explicit methods):

$$ \frac {\text{c} \Delta t}{\Delta z} < \text{C}_{max} = 1 $$

It makes sense that if I'm considering a domain of $ 0 < z < 100 \text{ m}$ where $\Delta z \sim 1 \text{m}$, then $\Delta t$ should be on the order of $10^{-9} \text{s} $. Although is the CFL condition assuming my $\Delta z$ is at an appropriate spatial width to account for changes in $E$ over space? Since if I made $\Delta z$ arbitrarily large (e.g. $10^9 \text{m}$) to consider a larger domain in $z$ now, then I could set $\Delta t$ to $1\text{s}$ and satisfy the CFL condition. Although based on (1) and the initial conditions, we know $E = \sin(z-ct)$, and would need to be sampled at a spatial frequency of $1\text{m}^{-1}$ and time frequency of GHz to be evaluated correctly.

To raise my question: it appears the CFL condition is only a necessary but not sufficient condition. But in a more general case where one does not necessarily know how rapidly $E$ will vary with $z$ (e.g. a more complicated hyperbolic PDE where $E = \text{A}(z,t)\text{B}(z,t)\text{C}(z,t)$ where $\text{A, B,}$ and $\text{C}$ are variable coefficients that vary with $t$ and $z$ on unknown scales scales (e.g. stiff)), how can we define $\Delta z$ a priori when applying the method of lines?

For example, if I was to advance in $z$ as opposed to $t$ (i.e. discretizing $t$) which is the contrary to the typical method of lines, I would have a condition reading:

$$ \frac {\Delta z}{\text{c} \Delta t} < \text{C}^{z}_{max} = 1 $$

which would thus place a limit on how small $\Delta z$ should be based on $\Delta t$ (which I define ahead of time) if I want to model the spatial frequency of $E$ correctly. But once again, this assumes I have chosen an appropriate $\Delta t$ for the problem—it is the discretized variable I am not quite understanding how to define its step since we can always implement an adaptive step for the other variable. But how can one implement an adaptive step for both independent variables when these equations involve derivatives of the other variable?

Also, is there a well-defined value for $\text{C}_{max}$ when using implicit methods (e.g. BDF)?

Answer 1

I have two extra points I would like to add to Wolfgang's answer. A formulation of the CFL condition that I find more useful than the classic formula is this:

A necessary condition for the stability of a numerical scheme is that the numerical domain of dependence bounds the physical domain of dependence.

This is exactly what good old $$ \dfrac{\Delta t}{c \Delta x} \leq 1 $$

says for simple schemes. The quote I wrote above has a bit more nuance, though, that can be used to answer the questions about implicit methods and higher order algorithms (my two points).

My points below are for the first order wave equation just to make things easy: I hope that everyone can see how the same ideas apply to the second order equation.

Higher order algorithms

It is true that for most time stepping methods increasing the order imposes stricter stability constraints but a notable exception to this is explicit Runge Kutta methods (which include some part of the imaginary axis in their stability regions; that is a separate stability concern). Consider the first order wave equation $$ u_t + c u_x = 0 $$ with periodic boundary conditions and second order finite differences $$ u_x(x_j) = \dfrac{u(x_j + \Delta x) - u(x_j - \Delta x)}{2 \Delta x} + O(\Delta x^2) $$

(to make things as easy as possible) in space. One can write the semidiscrete system as $$ \dfrac{d\vec{U}}{dt} = A \vec{U}. $$ If we use RK3 then we have $$ \begin{align} k_1 &= A \vec{U}\,^n \\ k_2 &= A \left(\vec{U}\,^n + \dfrac{\Delta t}{2} k_1\right) \\ k_3 &= A \left(\vec{U}\,^n + \dfrac{3 \Delta t}{4} k_2\right) \\ \vec{U}\,^{n + 1} &= \vec{U}\,^n + \dfrac{\Delta t}{9} (2 k_1 + 3 k_2 + 4 k_3) \\ \end{align} $$ so, by definition of (A):

• the $j$th entry of $k_1$ uses $U_{j - 1}^n$ and $U_{j + 1}^n$
• the $j$th entry of $k_2$ uses $U_{j - 2}^n$, $U_{j - 1}^n$, $U_j^n$, $U_{j + 1}^n$, and $U_{j + 2}^n$
• the $j$th entry of $k_3$ uses $U_{j - 3}^n$ through $U_{j + 3}^n$

since for each stage we use the same space operator on the current solution which effectively widens the stencil. Hence doing one step with this method has the CFL constraint $$ \dfrac{c \Delta t}{\Delta x} \leq 3 $$ since calculating $U_j^{n + 1}$ requires data from three grid points to the left and three grid points to the right. The same argument holds for higher order methods: a method with more stages tends to have a wider stencil and therefore a larger Courant number.

Warning: due to the spatial operator we used this scheme is not stable, even if it satisfies the CFL constraint. The correct Courant number to use (which comes from the stability region of RK3) is (roughly)

$$ \dfrac{c \Delta t}{\Delta x} \leq 1.75. $$

Implicit methods

The answer to your question about ''well-defined values for $C_{max}$ when using implicit methods'' is that $C_{max} = \infty$ by the same plain-text version of the CFL condition. To see why this is true consider the same ODE system with backward Euler: $$ \vec{U}\,^{n + 1} = \vec{U}\,^{n} + \Delta t A \vec{U}\,^{n + 1} \Rightarrow (I - \Delta t A) \vec{U}\,^{n + 1} = \vec{U}\,^{n} $$ so solving for the solution at $t^{n + 1}$ requires knowledge of all values at $t^n$, meaning that no matter what time step we use we are pulling in the entire physical domain of dependence. Implicit methods for wave equations are an unusual choice for a variety of (unrelated to the original question) reasons but they do satisfy the CFL condition for any choice of $\Delta t > 0$.

Edit:

@Mathews24 asked some more closely related questions whose answers won't really fit in a comment so I am appending them to my answer here. I hope that that is an acceptable way to answer questions (please correct me if there is a better place for this).

Why is the centered difference scheme unstable past a Courant number of $1.75$?

As mentioned previously the CFL condition is a necessary, but not sufficient, condition for the stability of a numerical scheme. The other usual stability condition comes from the ODE solution method (which depends on the eigenvalues of $A$, which in turn depends on the spatial discretization: this is why I say that the choice of spatial discretization impacts the stability of the time discretization). Ignoring the physical meaning of things for a moment, we say that an ODE solver for a scalar problem $$ y' = \lambda y + f $$ is stable if the numerical approximation of the unforced problem, $y' = \lambda y$, remains bounded for all time. I like to think of this thermodynamically: if one considers an unforced wave in a closed system then that wave will not grow in energy over time (but we allow it to shrink for the purposes of our numerical scheme, since that is better than exponential growth). An easy example of this is forward Euler: one can check that if we do $$ w^n \approx y(t^n), \,w^{n + 1} = w^n + \lambda \Delta t w^n \Rightarrow w^{n + 1} = (1 + \lambda \Delta t) w^n $$

that the scheme is stable when (for complex $\lambda$) $|\lambda \Delta t + 1| \leq 1$. The value on the left is called the amplification factor: the solution will grow exponentially if it is higher than one. The same idea holds for all other ODE schemes (in particular Runge-Kutta methods): one can calculate a stability region for the three stage scheme I presented above with a restriction on $\lambda \Delta t$ that guarantees that the solution does not grow. The stability plot is the blue curve under ''Bogacki-Shampine-4-2-3'' at

http://runge.math.smu.edu/arkode_dev/doc/guide/build/html/Butcher.html

if anyone wants to see it.

The rule of thumb that people use (that works for all non-pathological cases; there is a nice book by Trefethen and Reddy about how to do this more generally) to examine stability of a system of ODEs $y' = A y$ is to do 1D stability analysis with all of the eigenvalues of $A$. This is rigorous if the matrix is diagonalizable and usually works even if it is not. If one looks at the stability region of the RK3 scheme then one will see that it stays (barely!) to the right of the imaginary axis and intersects at about $1.75$. For the matrix I gave above the eigenvalues are purely imaginary and are bounded (in modulus) by

$$ |\text{eig}(A)| \leq \dfrac{c}{\Delta x} $$

which is exactly why things blow up for a Courant number greater than $1.75$.

Why are implicit methods unusual for wave equations?

Running this problem with RK3 near the stability limit is nice because the energy in the system is very nearly conserved: things decay very slightly over time since all of the amplification factors are very close to $1$ (the product $\lambda_i \Delta t$ is always near the boundary of the stable region). If we used an implicit method then we won't have to worry about the CFL constraint (or, for the right choice of implicit method, the amplification factors) but one can check that as a result of this the amplification factors will be significantly less than one, leading to a lot of extra dissipation in the scheme: this is usually not what you want with waves (usually some sort of conservation is desired, even if it is only approximate). I recommend calculating the amplification factors of backward Euler along the imaginary axis to convince oneself of why this is a concern.

Answer 2

You are correct: If you satisfy the CFL condition, then all that guarantees is that your scheme is stable, i.e., the numerical solution does not go to infinity. But the CFL condition says nothing about how accurate the numerical solution is. For that, indeed $\Delta z$ and $\Delta t$ must also be small enough compared to the features of the exact solution.

I'm going to note that $C_\text{max}$ may not be 1, but can also be smaller (sometimes much smaller) depending on what time stepping scheme you use. Generally, higher order schemes will yield smaller values for $C_\text{max}$.