At the end of the last post I described a problem I did not actually explain any further. During rollout the network eats its own predictions: a slightly wrong output becomes the input for the next step, that step inherits the error and adds a little of its own, and over a long horizon the whole thing (error) can drift away from physical reality or blow up entirely. I wrote that the errors "compound". I did not say why they compound rather than simply pile up, and I did not say what decides whether a given scheme is the kind that amplifies a small perturbation or the kind that quietly kills it.

This is where you stumble upon the CFL condition as the rule that dictates how small a time step a wave solver is allowed to take. This topic is indirectly tied to the von Neumann stability analysis. It turns out to be the machinery that produces those conditions, and it answers exactly the question we've been unable to answer about rollout. This post will work through it from scratch.

The CFL condition has a genuinely intuitive picture behind it, so let me state that first.

\(\frac{c \Delta t}{\Delta x} \leq 1\)0

where \(c\) denotes the wave speed, \(\Delta t\) denotes the size of one time step, and \(\Delta x\) denotes the spacing between two adjacent grid points. The left-hand side is called the Courant number.

In the continuous wave equation a disturbance travels at a fixed speed \(c\). Over one time step of size \(\Delta t\), the physical disturbance therefore covers a distance \(c \Delta t\) per step. The numerical scheme of the finite-difference stencil that approximates the the 2nd spatial derivative, meanwhile, computes each point's new value by looking only at that point and its two immediate neighbours, which sit a distance \(\Delta x\) away as \(\frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2}\).

The set of places whose values can possibly influence a given point after one step is called that point's domain of dependence. For the true PDE it is everything within \(c\Delta t\). For the scheme it is everything within \(\Delta x\) around a given point. If \(c\Delta t > \Delta x\), then real physical information arrives from farther away than the scheme ever looks. The scheme is being asked to produce a number that depends on data it never reads. It cannot possibly get the right answer, as it's lacking information. That is the intuition.

Two questions structure everything below:

1. Does an error present in the field grow or shrink as the scheme runs, and what exactly does that depend on?
2. Turns out, if it grows, it grows exponentially rather than additively — so why is that?

Everything in this post is bookkeeping between three things, and the whole derivation falls apart if they get blurred together. So let me name them explicitly.

\(u_i^n\) denotes the true solution — the exact value the real PDE would have at grid point \(i\) at time step \(n\), where \(i\) denotes the index of a point along the spatial grid and \(n\) denotes how many time steps have been taken so far.

\(U_i^n\) denotes the numerical solution — the value the computer actually holds in memory at that same grid point and time step, produced by running the Euler scheme from the previous posts.

\(e_i^n := U_i^n - u_i^n\)(1)

where \(e_i^n\) denotes the error at grid point \(i\) and time step \(n\), \(U_i^n\) denotes the numerical value there, and \(u_i^n\) denotes the true value there.

For explicity reasons note that the superscript \(n\) here is an index, not a power. Nothing is being raised to anything.

I want a rule that tells me how the error at one time step becomes the error at the next. Feel free to skip this section if you are willing to just accept, that the error obeys the same update rule as the field itself. I will show it anyway, because it is a nice little exercise in bookkeeping and it is worth seeing how the pieces fit together.

Start with the numerical solution, which is the Euler scheme applied to the heat equation (moved the \(\Delta x^2\) divisor into the coefficient as it will be more convenient for later):

\(U_i^{n+1} = U_i^n + \frac{\alpha \Delta t}{\Delta x^2}\left(U_{i+1}^n - 2U_i^n + U_{i-1}^n\right)\)(2)

where \(U_i^{n+1}\) denotes the numerical value at grid point \(i\) at the next time step, \(U_i^n\) denotes its numerical value now, \(\alpha\) denotes the thermal diffusivity from the heat equation, \(\Delta t\) denotes the time-step size, \(\Delta x\) denotes the grid spacing, and \(U_{i+1}^n\) and \(U_{i-1}^n\) denote the numerical values at the right and left neighbours at the current time step. Nothing new here.

Now the true solution. For the moment I am going to pretend that the true solution satisfies the same-shaped formula:

\(u_i^{n+1} = u_i^n + \frac{\alpha \Delta t}{\Delta x^2}\left(u_{i+1}^n - 2u_i^n + u_{i-1}^n\right)\)(3)

where every symbol means the same as above except that each \(U\) has been replaced by the corresponding true value \(u\). This is not quite right: the true solution satisfies the continuous heat equation, and the finite-difference stencil is only an approximation of the continuous second derivative, so there is a small mismatch. But we are ignoring that mismatch for now and will add it back explicitly near the end of the post, once the machinery exists to say something useful about it.

Now subtract the second equation from the first, term by term:

\(U_i^{n+1} - u_i^{n+1} = \left(U_i^n - u_i^n\right) + \frac{\alpha \Delta t}{\Delta x^2}\left(\left(U_{i+1}^n - u_{i+1}^n\right) - 2\left(U_i^n - u_i^n\right) + \left(U_{i-1}^n - u_{i-1}^n\right)\right)\)(4)

where each bracketed difference is a numerical value minus the true value at one specific grid point and time step, and the coefficient \(\frac{\alpha \Delta t}{\Delta x^2}\) denotes the same constant as before, which survives the subtraction untouched because it multiplies both bracketed groups identically.

Every one of those differences is by definition an error. Replacing them gives:

\(e_i^{n+1} = e_i^n + \frac{\alpha \Delta t}{\Delta x^2}\left(e_{i+1}^n - 2e_i^n + e_{i-1}^n\right)\)(5)

where \(e_i^{n+1}\) denotes the error at grid point \(i\) at the next time step, \(e_i^n\) denotes the error there now, \(e_{i+1}^n\) and \(e_{i-1}^n\) denote the errors at the right and left neighbours now, and the coefficient is unchanged.

The error obeys the exact same update rule as the field itself. The structural reason is worth stating plainly, because it is the only reason any of this works: the update rule is linear in the field values. It multiplies values by constants and adds them, and it does nothing else. Pushing two quantities through the same linear rule and then subtracting the results is the same as subtracting first and pushing the difference through once.

So I have a rule for how errors evolve. The obvious next move is to ask whether that rule makes errors bigger or smaller, and the obvious way to answer it is to try some errors and see.

That does not work, and it is worth seeing why it does not, because the failure is what motivates everything else.

A real error vector is some arbitrary shape, say \([0.3,\, -1.2,\, 0.7,\, 0.1,\, -0.4]\), where each entry denotes the error at one grid point. Push it through the update once and you get another arbitrary shape. There is no visible relationship between the input list and the output list, because every output entry is built by mixing three different input entries. The shape does not survive; it turns into a different shape. Doing it again gives a third unrelated shape. Nothing about this tells me whether errors in general grow.

So what I need is a way to break an arbitrary error into pieces that do keep their shape under the update. Then I only have to understand the pieces.

Pull the spatial part of the update out and give it a name:

\(S(e)_i = e_{i+1} - 2e_i + e_{i-1}\)(6)

where \(S\) denotes the stencil operator — the thing that takes a whole vector of values and returns a whole new vector — \(S(e)_i\) denotes entry \(i\) of the result, \(e_i\) denotes the input value at that point, and \(e_{i+1}\) and \(e_{i-1}\) denote its right and left neighbours. This is the same stencil derived from scratch in the PDE post — how much a point differs from the combined total of its two neighbours — just written without the \(\Delta x^2\) divisor, which we can safely ignore for the moment as it is just a constant factor.

The property we need from it is linearity:

\(S(P + Q) = S(P) + S(Q)\)(7)

where \(P\) and \(Q\) denote two arbitrary vectors of grid values, \(P + Q\) denotes their entry-by-entry sum, and \(S\) denotes the stencil operator. In words: applying the stencil to a sum gives the same answer as applying it separately and adding afterwards.

One clarifications, I wanted spelled out.

Linear does not mean an operation has to be applied elementwise to the input vector. In this case for example, output entry \(i\) is built from three input entries, so information genuinely moves sideways across the grid. Linearity says nothing about points staying independent of each other. It only says that whole vectors can be split apart and recombined without loss.

That is the payoff, and it is worth stating as a plan, because the rest of the post is just executing it:

1. Pick a family of simple building-block shapes.
2. Write any messy error vector as an exact weighted sum of them.
3. Work out what the stencil does to each building block on its own.
4. Add the results back together, knowing that linearity guarantees no cross-terms were missed.

An intractable problem — "check every possible error" — becomes a tractable one: "check every building block."

Before choosing building blocks, I need a way of saying that one error is bigger than another. The standard choice is the squared L² norm:

\(\|e^n\|^2 = \sum_i \left(e_i^n\right)^2\)(8)

where \(\|e^n\|^2\) denotes the squared total size of the error field at time step \(n\), \(e_i^n\) denotes the error at grid point \(i\) at that time step, and \(\sum_i\) denotes the sum over every grid point.

The reason for squaring rather than just adding the values up is that raw sums cancel. Take \(e = [5,\, -5]\), where each entry denotes the error at one of two grid points. Summed directly it gives \(0\), which would claim there is no error at all, when in fact both points are badly wrong. Squared and summed it gives \(50\). Squaring removes the sign, so opposite mistakes stop hiding each other.

Any error vector, however messy, can be written as a weighted sum of simple wave shapes — sines and cosines of various spatial frequencies. If each of those shapes comes back out of the stencil rescaled but otherwise unchanged, then I know everything. I can decompose any error into its shapes, multiply each one by its own scaling factor, and add them back up.

The important word is exact. This is not an approximation or a truncation.

Take an 11-point grid, \(i = 0, 1, \ldots, 10\), and put a sine on it:

\(T_i = \sin\!\left(\frac{i\pi}{2}\right)\)(9)

where \(T_i\) denotes the value of the test shape at grid point \(i\), \(i\) denotes the grid index, and \(\frac{\pi}{2}\) denotes how much the sine's angle advances per grid point. Evaluating it gives \(T = [0,\, 1,\, 0,\, -1,\, 0,\, 1,\, 0,\, -1,\, 0,\, 1,\, 0]\).

Now apply the stencil

\(S(T)_i = T_{i+1} - 2T_i + T_{i-1}\)

at a few points and compare each result against the input value at the same point.

At \(i = 1\): \(0 - 2(\sin\!\left(\frac{\pi}{2}\right)\)) + 0 = -2), and \(-2 \cdot T_1 = -2\).

At \(i = 2\): \(-1 - 0 + 1 = 0\), and \(-2 \cdot T_2 = 0\).

At \(i = 3\): \(0 - 2(\sin\!\left(3\frac{\pi}{2}\right)\)) + 0 = 2), and \(-2 \cdot T_3 = 2\).

Every time, the output at a point is the input at that same point multiplied by \(-2\):

\(S(T)_i = -2 \cdot T_i\)(10)

where \(S(T)_i\) denotes the stencil applied to the sine shape and read off at grid point \(i\), \(T_i\) denotes the sine shape's own value at that point, and \(-2\) denotes the constant multiplier that came out. The shape went in, the same shape came out, only scaled. That is exactly the behaviour we want, and it is not an accident of the numbers — the next few sections show where the multiplier comes from and how to get it for any frequency without evaluating anything by hand.

Here you can see what happens if you apply the finite-difference stencil to a Fourier mode. You'll see what that means throughout the rest of the post. But for now, just play around with the controls and note, how a wave shaped error simply gets scaled by the stencil, whereas more complex error shapes seem to get transformed in a more complicated way.

Doing the above with sines works but is clumsy, because shifting a sine by one grid point turns it into a mixture of sine and cosine and the algebra gets messy fast (I honestly don't even know how that would work, so this was just Claude speaking, but we'll give him the benefit of the doubt as it seems like a reasonable argument). There is a form that makes shifting trivial:

\(T_i = e^{j\kappa i}\)(11)

where \(T_i\) denotes the value of the building-block shape at grid point \(i\), \(\kappa\) denotes the wavenumber (some of you mightz also know it as the 'winding frequency' as often mentioned by channels such as 3blue1brown is the Fourier Transform video) — the spatial frequency, meaning how much the wave's angle advances from one grid point to the next — \(i\) denotes the grid index, and \(j\) denotes the imaginary unit, the number whose square is \(-1\).

A word on that \(j\). Mathematicians normally write the imaginary unit as \(i\). I am not this time, because \(i\) is already the grid index in every formula in this post and has been since the PDE post.

Euler's formula unpacks what that exponential actually is:

\(T_i = e^{j\kappa i} = \cos(i\kappa) + j\sin(i\kappa)\)(12)

where \(T_i\) denotes the shape's value at grid point \(i\), \(\cos(i\kappa)\) denotes its real part, \(\sin(i\kappa)\) denotes its imaginary part, \(\kappa\) denotes the wavenumber and \(j\) denotes the imaginary unit. So a complex exponential is not some exotic new object — it is a cosine and a sine carried around together in one symbol.

Check what that gives for \(\kappa = \pi/2\): the values are \(T = [1,\, j,\, -1,\, -j,\, 1,\, \ldots]\), where each entry denotes the shape's complex value at one grid point. Take just the imaginary parts and you get \([0,\, 1,\, 0,\, -1,\, 0,\, \ldots]\) — literally the plain sine from the previous section. The sine case is the imaginary part of the complex case. Nothing has been abandoned by going complex; the sine is still in there.

This is the reason the complex form is worth the trouble, and it rests entirely on the ordinary exponent rule \(e^{a+b} = e^a \cdot e^b\), where \(a\) and \(b\) denote any two exponents. Apply it to the shape shifted one point to the right:

\(T_{i+1} = e^{j\kappa(i+1)} = e^{j\kappa i} \cdot e^{j\kappa} = T_i \cdot e^{j\kappa}\)(13)

where \(T_{i+1}\) denotes the shape's value at the right neighbour, \(T_i\) denotes its value at the point itself, \(\kappa\) denotes the wavenumber, and \(e^{j\kappa}\) denotes a fixed complex number that does not depend on \(i\) at all.

The same argument to the left:

\(T_{i-1} = e^{j\kappa(i-1)} = T_i \cdot e^{-j\kappa}\)(14)

where \(T_{i-1}\) denotes the shape's value at the left neighbour.

Check both with numbers at \(\kappa = \pi/2\), where \(e^{j\pi/2} = j\) and \(e^{-j\pi/2} = -j\). From \(T_1 = j\), the rule predicts \(T_2 = j \cdot j = -1\), which matches the list above, and \(T_0 = j \cdot (-j) = 1\), which also matches.

Now that we summarized all necessary tools and information, let's feed the complex exponential into the stencil and substitute the two shift rules:

\(S(T)_i = T_i e^{j\kappa} - 2T_i + T_i e^{-j\kappa} = T_i\left(e^{j\kappa} - 2 + e^{-j\kappa}\right)\)(15)

where \(S(T)_i\) denotes the stencil applied to the shape and evaluated at grid point \(i\), \(T_i\) denotes the shape's value at that point, and \(e^{j\kappa}\) and \(e^{-j\kappa}\) denote the right-shift and left-shift multipliers.

The bracket has no \(i\) in it, so give it a name:

\(g(\kappa) := e^{j\kappa} - 2 + e^{-j\kappa}\)(16)

where \(g(\kappa)\) denotes the factor by which the stencil rescales the shape of wavenumber \(\kappa\), and the three terms denote the right neighbour's contribution, the point's own contribution and the left neighbour's contribution respectively. It depends only on \(\kappa\) — never on the position \(i\), never on the time step \(n\).

Check it against the hand-computed sine result. At \(\kappa = \pi/2\): \(g(\pi/2) = j - 2 - j = -2\). That is exactly the \(-2\) I got by evaluating the stencil point by point earlier.

The Von Neumann analysis however is not solely bound to the 2nd order-derivative stencil. Take a first-derivative-style stencil \(T_{i+1} - T_i\), where the two terms denote the right neighbour and the point itself, and the identical substitution gives \(T_i(e^{j\kappa} - 1)\). Von Neumann analysis works for any linear, constant-coefficient finite-difference scheme. Only the formula for \(g\) changes.

The complex form is correct but hard to reason about, so simplify it. Euler's formula both ways round gives \(e^{j\kappa} = \cos\kappa + j\sin\kappa\) and \(e^{-j\kappa} = \cos\kappa - j\sin\kappa\), where \(\cos\kappa\) denotes the real part, \(\sin\kappa\) the imaginary part and \(\kappa\) the wavenumber. Adding them, the imaginary parts cancel exactly:

\(e^{j\kappa} + e^{-j\kappa} = 2\cos(\kappa)\)(17)

where the left-hand side denotes the sum of the two shift multipliers and \(2\cos(\kappa)\) denotes what remains once the two \(j\sin\kappa\) terms annihilate each other. Substituting into the definition:

\(g(\kappa) = 2\cos(\kappa) - 2\)(18)

where \(g(\kappa)\) denotes the stencil's rescaling factor for wavenumber \(\kappa\), \(\cos(\kappa)\) denotes the cosine of the wavenumber, and the \(-2\) comes from the point's own term in the stencil. This is now a plain real-valued function of one variable, which is a much better thing to reason about than a sum of complex exponentials.

Check again at \(\kappa = \pi/2\): \(2(0) - 2 = -2\). Consistent with both earlier computations.

Linearity says a sum of two shapes should evolve as the sum of the two evolutions, but I wanted to see that with two different frequencies in the same vector, since that is the situation a real error is in.

Mode A: \(T_i = e^{ji\pi/2}\) with amplitude \(3\), giving \(g(\pi/2) = -2\). Mode B: \(U_i = e^{ji\pi} = (-1)^i\) with amplitude \(2\), giving \(g(\pi) = -1 - 2 - 1 = -4\), where \(T_i\) and \(U_i\) denote the two shapes' values at grid point \(i\) and \(g\) denotes the stencil's rescaling factor at each of the two wavenumbers.

The combined vector is \(W_i = 3T_i + 2U_i\), where \(W_i\) denotes the mixed value at grid point \(i\). Written out for \(i = 0 \ldots 4\) it is \(W = [5,\; -2+3j,\; -1,\; -2-3j,\; 5]\).

Apply the stencil directly at \(i = 1\): \(S(W)_1 = (-1) - 2(-2+3j) + 5 = 8 - 6j\).

Now compute the two modes separately and add: \(3\,g(\pi/2)\,T_1 + 2\,g(\pi)\,U_1 = 3(-2)(j) + 2(-4)(-1) = -6j + 8\). Identical.

This is the conclusion that licenses everything that follows: each wavenumber evolves completely independently, forever. A messy real error is a sum of many modes, each one silently growing or shrinking by its own factor, and they never talk to each other or leak into one another.

Now put the time dimension back in. The assumption — the ansatz — is that the error of a given mode/wavenumber at any time is one fixed spatial shape scaled by a single number that changes from step to step:

\(e_i^n = \hat{e}^{\,n} \cdot T_i, \qquad T_i = e^{j\kappa i}\)(19)

where \(e_i^n\) denotes the error at grid point \(i\) at time step \(n\), \(T_i\) denotes the fixed spatial shape of wavenumber \(\kappa\) — identical at every time step, it never changes — and \(\hat{e}^{\,n}\) denotes the amplitude at time step \(n\): a single number, the same across every grid point for a given \(n\), different for different \(n\).

This is the centre of the whole method. I have a rule for how the error evolves, and I have a claimed form for what the error looks like. Put the second into the first.

The recurrence contains four different error terms, and each one gets replaced:

\(e_i^{n+1} = \hat{e}^{\,n+1} T_i\), \(e_i^n = \hat{e}^{\,n} T_i\), \(e_{i+1}^n = \hat{e}^{\,n} T_i e^{j\kappa}\), and \(e_{i-1}^n = \hat{e}^{\,n} T_i e^{-j\kappa}\), where \(\hat{e}^{\,n+1}\) denotes the amplitude at the next time step, \(\hat{e}^{\,n}\) denotes the amplitude at the current one, \(T_i\) denotes the fixed shape at grid point \(i\), and \(e^{j\kappa}\) and \(e^{-j\kappa}\) denote the right- and left-shift multipliers derived earlier.

Dropping all four into the error recurrence:

\(\hat{e}^{\,n+1} T_i = \hat{e}^{\,n} T_i + \frac{\alpha \Delta t}{\Delta x^2}\left(\hat{e}^{\,n} T_i e^{j\kappa} - 2\hat{e}^{\,n} T_i + \hat{e}^{\,n} T_i e^{-j\kappa}\right)\)(20)

where \(\hat{e}^{\,n+1}\) and \(\hat{e}^{\,n}\) denote the amplitudes at the next and current time steps, \(T_i\) denotes the fixed spatial shape at grid point \(i\), \(\alpha\) denotes the diffusivity, \(\Delta t\) the time step, \(\Delta x\) the grid spacing, and the two exponentials the shift multipliers.

The goal is to separate the part that depends on position and time from the part that depends only on the wavenumber. Every single term on the right contains \(\hat{e}^{\,n} T_i\), so pull it out:

\(\hat{e}^{\,n+1} T_i = \hat{e}^{\,n} T_i \left[1 + \frac{\alpha \Delta t}{\Delta x^2}\left(e^{j\kappa} - 2 + e^{-j\kappa}\right)\right]\)(21)

where the bracket denotes everything that survived the factoring, \(1\) comes from the standalone \(\hat{e}^{\,n} T_i\) term, and the remaining symbols are as above. The bracket contains no \(i\) and no \(n\).

The goal is to get rid of the complex numbers, since the bracket is exactly \(g(\kappa)\) and we already simplified that to a real expression:

\(\hat{e}^{\,n+1} T_i = \hat{e}^{\,n} T_i \left[1 + \frac{\alpha \Delta t}{\Delta x^2}\left(2\cos(\kappa) - 2\right)\right]\)(22)

where \(2\cos(\kappa) - 2\) denotes the simplified stencil factor \(g(\kappa)\), \(\kappa\) denotes the wavenumber, and all other symbols are unchanged.

The goal here is to remove space from the problem entirely. \(T_i\) appears once on the left and once on the right, in exactly the same form, so it can be divided out. It is never zero: \(|e^{j\kappa i}| = 1\) for every \(i\), where the left-hand side denotes the magnitude of the complex shape value, so there is no division by zero anywhere on the grid.

\(\hat{e}^{\,n+1} = \hat{e}^{\,n} \left[1 + \frac{\alpha \Delta t}{\Delta x^2}\left(2\cos(\kappa) - 2\right)\right]\)(23)

where \(\hat{e}^{\,n+1}\) denotes the amplitude at the next time step, \(\hat{e}^{\,n}\) the amplitude at the current step, and the bracket the same wavenumber-dependent factor as before. The grid index is gone. What is left is a pure scalar recurrence for the amplitude of a single mode, based only on the wavenumber, the diffusivity, the time step and the grid spacing.

\(G(\kappa) := 1 + \frac{\alpha \Delta t}{\Delta x^2}\left(2\cos(\kappa) - 2\right)\)(24)

where \(G(\kappa)\) denotes the growth factor for the mode of wavenumber \(\kappa\), \(\alpha\) denotes the diffusivity, \(\Delta t\) the time step, \(\Delta x\) the grid spacing, and \(2\cos(\kappa) - 2\) the stencil factor. With that name, the recurrence reads:

\(\hat{e}^{\,n+1} = G(\kappa) \cdot \hat{e}^{\,n}\)(25)

where \(\hat{e}^{\,n+1}\) denotes the amplitude one step later, \(\hat{e}^{\,n}\) denotes the amplitude now, and \(G(\kappa)\) denotes the growth factor.

The thing to hold onto: \(G(\kappa)\) is built from \(\kappa\), \(\alpha\), \(\Delta t\) and \(\Delta x\) — all of which are fixed before the simulation starts. It never depends on \(n\). For a given mode it is the same number at step 1 and at step 10,000.

Take \(G(\kappa) = 0.8\) as given — the point here isn't where that number comes from, it's what happens when the scheme applies it over and over. Each step multiplies the current error amplitude by \(G\) again:

\(\hat{e}^{\,n} = G(\kappa)^n \cdot \hat{e}^{\,0}\)(26)

where \(\hat{e}^{\,n}\) denotes the error amplitude after \(n\) steps and \(\hat{e}^{\,0}\) denotes the starting amplitude.

One might have assumed error accumulation means each step adds a fixed chunk of new error to a running total. It doesn't. The same multiplier gets reapplied every step. With \(G = 1.05\), after 100 steps the error is multiplied by roughly \(131\); after 200 steps, by roughly \(17{,}000\). A 5% overshoot per step is not a 5% problem.

Given \(\hat{e}^{\,n} = G(\kappa)^n \hat{e}^{\,0}\), the condition is immediate. If \(|G(\kappa)| > 1\), repeated multiplication drives the amplitude to infinity. If \(|G(\kappa)| < 1\), it drives it to zero. The absolute value is there because \(G\) can be negative, and a factor of \(-1.4\) flips the sign every step while still growing in size.

Since modes never mix, no mode can be rescued by any other. So:

\(|G(\kappa)| \leq 1 \quad \text{for every } \kappa\)(27)

where \(|G(\kappa)|\) denotes the magnitude of the growth factor for the mode of wavenumber \(\kappa\). If even a single \(\kappa\) violates this, that one mode grows geometrically without bound and eventually dominates the entire field, no matter how tiny it started and no matter how beautifully every other mode behaves. Rounding error alone will seed it.

"Every \(\kappa\)" sounds like infinitely much work. It collapses in three moves.

On a discrete grid, \(\kappa\) and \(\kappa + 2\pi\) produce identical values at every grid point — adding a full turn around the unit circle changes nothing an integer grid index can see. So one window of width \(2\pi\) already covers every distinct mode, and the conventional choice is \(\kappa \in (-\pi,\, \pi]\).

\(G(\kappa) := 1 + \frac{\alpha \Delta t}{\Delta x^2}\left(2\cos(\kappa) - 2\right)\) touches \(\kappa\) only through \(\cos(\kappa)\), and \(\cos(-\kappa) = \cos(\kappa)\), where the two sides denote the cosine of a wavenumber and of its negative. For example \(\cos(-\pi/3) = \cos(\pi/3) = 0.5\). Positive and negative wavenumbers therefore give exactly the same growth factor, and the window shrinks to \(\kappa \in [0,\, \pi]\).

Sample the cosine across that interval at \(\kappa = 0,\, \pi/4,\, \pi/2,\, 3\pi/4,\, \pi\): the values are \(1,\; 0.71,\; 0,\; -0.71,\; -1\). Steadily decreasing.

Since the coefficient \(\frac{\alpha \Delta t}{\Delta x^2}\) is positive — all three of \(\alpha\), \(\Delta t\) and \(\Delta x^2\) are positive quantities — multiplying by it preserves the ordering. So \(G(\kappa)\) moves in lockstep with \(\cos(\kappa)\) and is monotonically decreasing across \([0, \pi]\) too.

A monotonic function on an interval takes its largest and smallest values at the two endpoints. Nothing can hide in the middle, because to be extreme in the middle it would have to turn around, and a monotonic function does not turn around. So checking all of \([0, \pi]\) reduces to checking \(\kappa = 0\) and \(\kappa = \pi\). Two numbers.

At \(\kappa = 0\): \(\cos(0) = 1\), so the stencil factor is \(2(1) - 2 = 0\), so the entire second term of \(G\) vanishes regardless of what \(\alpha\), \(\Delta t\) or \(\Delta x\) are:

\(G(0) = 1\)(28)

where \(G(0)\) denotes the growth factor of the zero-wavenumber mode — the constant mode, same value at every grid point. A stencil measures how a point differs from its neighbours, and a constant has no difference to measure, so it returns zero. A uniform offset in the error neither grows nor shrinks; it just sits there.

At \(\kappa = \pi\): \(\cos(\pi) = -1\), so the stencil factor is \(2(-1) - 2 = -4\):

\(G(\pi) = 1 - \frac{4\alpha \Delta t}{\Delta x^2}\)(29)

With \(\alpha = 1\), \(\Delta t = 0.1\), \(\Delta x = 1\) this gives \(G(\pi) = 0.6\).

It is worth saying what \(\kappa = \pi\) physically is. The shape is \(e^{ji\pi} = (-1)^i\), which alternates \(+1, -1, +1, -1\) from one grid point to the next. That is the sharpest pattern the grid is capable of representing at all — you cannot oscillate faster than flipping sign every single point. And it is this mode, the most oscillatory one, that fails first. That connects directly to the concern I raised in the wave-equation posts about sharp features: the roughest content in the field is always the first thing to go unstable.

\(G(0) = 1\) sits exactly at the boundary and never exceeds it, so it can never be the thing that breaks. The only way to fail is for \(G(\pi)\) to fall below \(-1\), since increasing \(\Delta t\) only pushes it downward. So require:

\(1 - \frac{4\alpha \Delta t}{\Delta x^2} \geq -1\)(30)

where the left-hand side denotes \(G(\pi)\) written out and \(-1\) denotes the lowest value a growth factor may take without its magnitude exceeding one.

Subtract \(1\) from both sides, which isolates the term containing the parameters:

\(-\frac{4\alpha \Delta t}{\Delta x^2} \geq -2\)(31)

where the left-hand side denotes the negative parameter term and \(-2\) denotes what remains on the right after subtracting.

Now multiply both sides by \(-1\) to clear the minus signs. This flips the inequality, because multiplying an inequality by a negative number reverses the direction of the comparison — \(3 > 2\) but \(-3 < -2\):

\(\frac{4\alpha \Delta t}{\Delta x^2} \leq 2\)(32)

Divide both sides by \(4\alpha\), which is positive so the direction stays:

\(\frac{\Delta t}{\Delta x^2} \leq \frac{1}{2\alpha}\)(33)

where \(\Delta t\) denotes the time step, \(\Delta x^2\) the squared grid spacing, and \(\alpha\) the diffusivity.

Finally multiply both sides by \(\Delta x^2\) to get the time step alone on the left, which is the form you actually use when setting up a simulation:

\(\Delta t \leq \frac{\Delta x^2}{2\alpha}\)(34)

where \(\Delta t\) denotes the largest time step the explicit scheme may take without blowing up, \(\Delta x\) denotes the grid spacing, and \(\alpha\) denotes the thermal diffusivity of the material.

Look at that fraction again. The time step is bounded by the grid spacing squared. Halve \(\Delta x\) and the allowed \(\Delta t\) drops to a quarter. So doubling the spatial resolution costs you twice as many grid points and four times as many time steps, which in one dimension is eight times the work, and worse in two or three dimensions.

That scaling law is the direct reason high-fidelity simulations are expensive, and it is the practical motivation behind this entire series. Generating training data for a neural PDE surrogate means running exactly these solvers, at exactly these resolutions, under exactly these constraints. When people say the data is the bottleneck for operator learning, this inequality is a large part of what they mean.

Earlier I assumed the true solution satisfies the discrete formula exactly. It does not — the stencil is a truncated Taylor expansion, so plugging the true solution in leaves a small leftover:

\(u_i^{n+1} = u_i^n + \frac{\alpha \Delta t}{\Delta x^2}\left(u_{i+1}^n - 2u_i^n + u_{i-1}^n\right) + \tau_i^n\)(38)

where \(u_i^{n+1}\) denotes the true value at the next time step and \(\tau_i^n\) denotes the truncation error — the gap left by cutting off the Taylor series, freshly generated at every point and every step, regardless of whether the scheme is stable.

Carrying this through the same mode decomposition as before, each step's amplitude picks up one extra term:

\(\hat{e}^{\,n+1} = G(\kappa) \cdot \hat{e}^{\,n} - \hat{\tau}^{\,n}\)(40)

where \(\hat{\tau}^{\,n}\) denotes the truncation error injected at step \(n\). Once injected, it's indistinguishable from ordinary error — it gets hit by \(G\) at every step from then on, same as anything already present. So \(\tau\) is just a constant drip of new error the scheme adds on its own, and whether that drip fades or snowballs is decided by the same \(G\) as everything else: \(|G| \leq 1\) and each dose decays away; \(|G| > 1\) and every dose ever injected keeps growing forever.

Everything in this post rested on one structural fact I did not draw attention to at the time: there was a single field. One error field, one amplitude \(\hat{e}^{\,n}\) per mode, one growth factor \(G(\kappa)\) per mode. "Does it blow up" reduced to "is \(|G| > 1\)", a comparison between two ordinary numbers.

The wave equation does not have one field. As established some post before, stepping it forward requires carrying displacement and velocity, two coupled grids, because a displacement snapshot alone cannot predict its own future. So the error is a pair of fields, and each mode has a pair of amplitudes — one for displacement error, one for velocity error. Both coupled update equations feed into both new amplitudes. That is where the Courant number appears, and with it the actual CFL condition I went looking for in the first place. Next post.