The last post ended by naming the one thing it could not do. The heat equation has a single field, so each Fourier mode of the error had a single amplitude, and the object that mapped that amplitude to the next step was a single number \(G(\kappa)\). Stability then reduced to the comparison \(|G(\kappa)| \leq 1\) between two ordinary numbers, and that comparison produced \(\Delta t \leq \frac{\Delta x^2}{2\alpha}\).

The wave equation is not stepped forward as one field. It is stepped forward as two coupled fields, displacement and velocity, for the reasons worked out two posts ago. So each mode carries two amplitudes, one for the displacement error and one for the velocity error, and the object mapping one pair to the next pair cannot be a single number. This post works out what it is instead. The Courant number appears on its own during that work, without being inserted by hand, and the CFL condition that opened the last post finally gets derived rather than asserted.

Everything below leans on results already established, so let me list them rather than rebuild them.

The three quantities: \(u_i^n\) denotes the true solution at grid point \(i\) and time step \(n\), \(U_i^n\) denotes the numerical solution there, and the error is their difference. The superscript \(n\) is an index, not a power.

The error obeys the same update rule as the field itself, because the update rule is linear in the field values: pushing two quantities through the same linear rule and subtracting afterwards is the same as subtracting first and pushing the difference through once.

You cannot test stability by feeding arbitrary error vectors into the scheme, because the stencil mixes three neighbours into every output point and a messy shape in gives an unrelated messy shape out. The way around it is Fourier decomposition: any error vector on the grid is an exact sum of modes \(T_i = e^{j\kappa i}\), where \(T_i\) denotes the value of the mode's spatial shape at grid point \(i\), \(\kappa\) denotes the wavenumber, and \(j\) denotes the imaginary unit (the letter \(i\) is already taken by the grid index throughout this series).

Shifting the grid index by one is multiplication by a fixed constant: \(T_{i+1} = T_i \cdot e^{j\kappa}\) and \(T_{i-1} = T_i \cdot e^{-j\kappa}\), where \(e^{j\kappa}\) and \(e^{-j\kappa}\) denote the right-shift and left-shift multipliers, neither of which depends on \(i\).

Feeding a mode through the bare second-derivative stencil rescales it by

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

It depends only on \(\kappa\) — never on the position \(i\), never on the time step \(n\).

The range of \(g(\kappa)\) is the single fact this post uses most often, so it is worth restating with its numbers. At \(\kappa = 0\) it is \(2(1) - 2 = 0\). At \(\kappa = \pi\) it is \(2(-1) - 2 = -4\). In between it decreases monotonically. So \(g(\kappa)\) is never positive and never smaller than \(-4\). The mode at \(\kappa = \pi\) is the checkerboard pattern \((-1)^i\), and it is always the first mode to go unstable.

Finally, truncation error \(\tau\), which denotes the leftover from cutting off the Taylor expansion, is injected fresh at every point and every step and is thereafter indistinguishable from ordinary error. It is generated identically whether the scheme is stable or not. Nothing in this post changes how much \(\tau\) appears; what it changes is what happens to it afterwards.

This is the two-field Euler rollout built in the post before last. The state is a displacement grid and a velocity grid, and both are stepped forward together:

\(V_i^{n+1} = V_i^n + \Delta t \, c^2 \, \frac{U_{i+1}^n - 2U_i^n + U_{i-1}^n}{\Delta x^2}\)(1) \(U_i^{n+1} = U_i^n + \Delta t \, V_i^n\)(2)

where \(U_i^n\) denotes the numerical displacement at grid point \(i\) at time step \(n\), \(V_i^n\) denotes the numerical velocity at the same point and step, \(c\) denotes the wave speed, \(\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 displacements at the right and left neighbours at the current step. The superscripts \(n\) and \(n+1\) are time indices, not exponents.

The detail that matters most here is the one that looks least important: the displacement update on the second line reads \(V_i^n\), the old velocity, not the value the first line just computed. Both updates read from the same frozen snapshot of the state. Neither one sees the other's fresh result.

The goal is a rule for how the pair of errors at one step becomes the pair of errors at the next — the same bookkeeping as last time, done twice.

Two fields means two error labels:

\(e_i^{\,u,n} := U_i^n - u_i^n, \qquad e_i^{\,v,n} := V_i^n - v_i^n\)(3)

where \(e_i^{\,u,n}\) denotes the displacement error and \(e_i^{\,v,n}\) the velocity error, both at grid point \(i\), time step \(n\).

The superscript now carries two labels separated by a comma: which field, then which time step. Neither is an exponent — \(e_i^{\,u,n}\) is not something raised to the power \(n\).

Subtracting the true velocity equation from the numerical one, with the coefficient surviving untouched and each bracketed difference becoming a displacement error:

\(e_i^{\,v,n+1} = e_i^{\,v,n} + \Delta t \, c^2 \, \frac{e_{i+1}^{\,u,n} - 2e_i^{\,u,n} + e_{i-1}^{\,u,n}}{\Delta x^2}\)(4)

the new velocity error, driven by the stencil applied to the current displacement error.

Same move on the displacement equation, which has no stencil at all:

\(e_i^{\,u,n+1} = e_i^{\,u,n} + \Delta t \, e_i^{\,v,n}\)(5)

the new displacement error, driven by the current velocity error.

That is the coupling, in both directions: velocity's next value depends on displacement's stencil, displacement's next value depends on velocity. Neither field evolves alone — neither equation says anything about growth by itself.

The ansatz is the same as before, applied twice:

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

where \(T_i\) denotes the fixed spatial shape of wavenumber \(\kappa\) evaluated at grid point \(i\), \(\hat{e}^{\,u,n}\) denotes the amplitude of the displacement error at time step n and \(\hat{e}^{\,v,n}\) denotes the amplitude of the velocity error at that step.

Two things about this need saying explicitly.

The shape is shared. The same \(T_i\) appears in both lines, with the same \(\kappa\). So the question is: does a single wavenumber reproduce itself in both equations at once, and by how much?

The amplitudes are separate. Displacement error and velocity error are different physical quantities in different units, and there is no reason for their amplitudes to agree. So each mode carries the pair \((\hat{e}^{\,u,n}, \hat{e}^{\,v,n})\).

Just think about what that means: we're essentially giving both error fields the possibility to include a certain fourier mode, however, the amplitude of that mode in the displacement field can be different from the amplitude of that mode in the velocity field. The shape is the same, but the amplitudes are independent. So the amplitude of one mode in one field might be large while the amplitude of that same mode in the other field might be 0.

Five error terms appear in the velocity recurrence and each gets replaced. The two neighbour terms additionally use the shift rules, so \(e_{i+1}^{\,u,n} = \hat{e}^{\,u,n} T_i e^{j\kappa}\) and \(e_{i-1}^{\,u,n} = \hat{e}^{\,u,n} T_i e^{-j\kappa}\). Written out in full before anything is tidied:

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

where \(\hat{e}^{\,v,n+1}\) denotes the velocity amplitude at the next step and \(e^{j\kappa}\) and \(e^{-j\kappa}\) denote the right- and left-shift multipliers.

The goal now is to remove space from the problem, exactly as in the heat case. Every term inside the bracket contains \(\hat{e}^{\,u,n} T_i\), so pull that out:

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

where the bracket denotes everything that survived the factoring and contains no \(i\) and no \(n\). That bracket is the definition of \(g(\kappa)\).

Now \(T_i\) sits as a factor in every single term, on both sides, so divide it out. The division is legal everywhere on the grid because \(|e^{j\kappa i}| = 1\) for every \(i\):

\(\hat{e}^{\,v,n+1} = \hat{e}^{\,v,n} + \frac{\Delta t \, c^2}{\Delta x^2}\, g(\kappa)\, \hat{e}^{\,u,n}\)(9)

where \(g(\kappa)\) denotes the stencil's rescaling factor for wavenumber \(\kappa\) and all other symbols are as above. Space is gone.

Same procedure, shorter. Substitute the ansatz into \(e_i^{\,u,n+1} = e_i^{\,u,n} + \Delta t \, e_i^{\,v,n}\), giving \(\hat{e}^{\,u,n+1} T_i = \hat{e}^{\,u,n} T_i + \Delta t \, \hat{e}^{\,v,n} T_i\), then divide by \(T_i\) on the same justification as before:

\(\hat{e}^{\,u,n+1} = \hat{e}^{\,u,n} + \Delta t \, \hat{e}^{\,v,n}\)(10)

No \(g(\kappa)\) appears in this one as there is no stencil.

Two numbers go in and two numbers come out, so collect them into a column vector:

\(\hat{e}^{\,n} = \begin{bmatrix} \hat{e}^{\,u,n} \\ \hat{e}^{\,v,n} \end{bmatrix}\)(11)

where \(\hat{e}^{\,n}\) denotes the pair of amplitudes for one wavenumber at time step \(n\).

Before writing down the matrix, it is worth restating what a matrix acting on a vector actually computes, because the four entries are read off by matching against exactly that pattern:

\(\begin{bmatrix} a & b \\ c & d \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} ax + by \\ cx + dy \end{bmatrix}\)(12)

where \(a, b, c, d\) denote the four matrix entries, \(x\) and \(y\) denote the two entries of the input vector, and the output's top entry \(ax + by\) is built from the top row while the bottom entry \(cx + dy\) is built from the bottom row.

The two amplitude equations, ordered displacement first, are \(\hat{e}^{\,u,n+1} = 1 \cdot \hat{e}^{\,u,n} + \Delta t \cdot \hat{e}^{\,v,n}\) and \(\hat{e}^{\,v,n+1} = \frac{\Delta t c^2}{\Delta x^2}g(\kappa) \cdot \hat{e}^{\,u,n} + 1 \cdot \hat{e}^{\,v,n}\). Matching the first against the top row gives \(a = 1\) and \(b = \Delta t\); matching the second against the bottom row gives \(c = \frac{\Delta t c^2}{\Delta x^2}g(\kappa)\) and \(d = 1\). So:

\(\hat{e}^{\,n+1} = M(\kappa)\, \hat{e}^{\,n}, \qquad M(\kappa) = \begin{bmatrix} 1 & \Delta t \\ \frac{\Delta t \, c^2}{\Delta x^2} g(\kappa) & 1 \end{bmatrix}\)(13)

where \(M(\kappa)\) denotes the amplification matrix of the mode with wavenumber \(\kappa\). Entry by entry: the top-left \(1\) is the displacement error carrying itself forward, the bottom-right \(1\) is the velocity error carrying itself forward, the top-right \(\Delta t\) is the displacement's dependence on the velocity, and the bottom-left \(\frac{\Delta t c^2}{\Delta x^2}g(\kappa)\) is the velocity's dependence on the displacement through the stencil.

In the heat case, repeated application gave \(\hat{e}^{\,n} = G^n \hat{e}^{\,0}\), where \(G^n\) denotes the growth factor raised to the power \(n\) and \(\hat{e}^{\,0}\) denotes the starting amplitude. Stability was then just "is \(|G| > 1\)".

The same repetition here gives

\(\hat{e}^{\,n} = M(\kappa)^n\, \hat{e}^{\,0}\)(15)

where \(M(\kappa)^n\) denotes the amplification matrix multiplied by itself \(n\) times and \(\hat{e}^{\,0}\) denotes the starting pair of amplitudes.

Raising a number to a power is trivial. Raising a matrix to a power is not. You do not raise the four entries individually — that is simply not what matrix multiplication does — and multiplying \(M(\kappa)\) by itself ten times produces four expressions with no visible pattern in them. So the goal of this whole section is to turn the question "does \(M(\kappa)^n\) blow up" back into a question about plain numbers raised to the power \(n\).

There is one kind of matrix for which powers are easy. Take

\(D = \begin{bmatrix} 5 & 0 \\ 0 & 3 \end{bmatrix}\)(16)

where \(D\) denotes a diagonal matrix, meaning both off-diagonal entries are zero.

Apply it to \([2, 4]^{\mathsf{T}}\), where the superscript \(\mathsf{T}\) denotes transposition, so this is a column vector written sideways to save space. Using the rule above, the result is \([5(2) + 0(4),\; 0(2) + 3(4)]^{\mathsf{T}} = [10, 12]^{\mathsf{T}}\). The top entry only ever met the \(5\) and the bottom entry only ever met the \(3\), because the zeros mean the rows never pull in the other coordinate.

Since that is true at every application, it is true after \(n\) of them:

\(D^n \begin{bmatrix} 2 \\ 4 \end{bmatrix} = \begin{bmatrix} 5^n \cdot 2 \\ 3^n \cdot 4 \end{bmatrix}\)(17)

where \(D^n\) denotes the diagonal matrix applied \(n\) times, and \(5^n\) and \(3^n\) denote the two diagonal entries raised to that power. Each coordinate independently becomes a plain number raised to \(n\). That is the easy scalar case from the heat post, running twice in parallel.

Whoever has watched the video series on linear algebra by 3Blue1Brown knows, how he visualized matrices by where they send the two basis vectors: the first column is where \(\hat{\imath} = [1, 0]^{\mathsf{T}}\) lands and the second column is where \(\hat{\jmath} = [0, 1]^{\mathsf{T}}\) lands, where \(\hat{\imath}\) denotes the unit vector along the horizontal axis and \(\hat{\jmath}\) denotes the unit vector along the vertical axis.

For \(D\), \(\hat{\imath}\) lands at \([5, 0]^{\mathsf{T}}\): five times longer, still exactly on the horizontal axis. And \(\hat{\jmath}\) lands at \([0, 3]^{\mathsf{T}}\): three times longer, still exactly on the vertical axis. Diagonal means both basis vectors land on their own line, only longer or shorter, never rotated away.

\(M(\kappa)\), however, is not diagonal. Its first column says \(\hat{\imath}\) lands at \([1,\; \frac{\Delta t c^2}{\Delta x^2}g(\kappa)]^{\mathsf{T}}\), which is off the horizontal axis entirely whenever \(g(\kappa) \neq 0\).

Feel free to skip this section if you already know what an eigenvector is. It is a vector that lands on its own line when multiplied by a matrix, only rescaled.

So the question becomes: even for a matrix that mixes the coordinates, is there some other direction — not \(\hat{\imath}\), not \(\hat{\jmath}\), that lands back on its own line, only rescaled?

Work it on a concrete and simple example matrix that will be carried through the rest of this section:

\(A = \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}\)(18)

where \(A\) is chosen such that its numbers are clean.

\(\hat{\imath}\) lands at \([2, 1]^{\mathsf{T}}\), which is rotated off the horizontal axis, so \(\hat{\imath}\) is not one of the special directions. But try \([1, 1]^{\mathsf{T}}\):

\(A\begin{bmatrix} 1 \\ 1 \end{bmatrix} = \begin{bmatrix} 2(1) + 1(1) \\ 1(1) + 2(1) \end{bmatrix} = \begin{bmatrix} 3 \\ 3 \end{bmatrix} = 3\begin{bmatrix} 1 \\ 1 \end{bmatrix}\)(19)

The output points in the same direction as the input and is three times as long. A vector with that property is called an eigenvector of the matrix, and the number it gets multiplied by is called the corresponding eigenvalue. Written generally:

\(Aw = \lambda w\)(20)

where \(A\) denotes the matrix, \(w\) denotes a nonzero vector, and \(\lambda\) denotes a plain number — the eigenvalue — such that applying the matrix to \(w\) has the same effect as multiplying \(w\) by that number.

Apply the matrix twice to an eigenvector:

\(A(Aw) = A(\lambda w) = \lambda (Aw) = \lambda(\lambda w) = \lambda^2 w\)(21)

where the second equality replaces \(Aw\) by \(\lambda w\) using the eigenvector property, and the third pulls the plain number \(\lambda\) out in front, which is allowed precisely because it is a number and not a matrix.

So for an eigenvector, applying the matrix \(n\) times is exactly raising a plain number to the power \(n\), which is exactly what we want.

However, the actual initial error pair \(\hat{e}^{\,0}\) will most likely not happen to be an eigenvector of \(M\).

But \(A\) for example has a second one (and so does \(M\)): \(A[1, -1]^{\mathsf{T}} = [2(1) + 1(-1),\; 1(1) + 2(-1)]^{\mathsf{T}} = [1, -1]^{\mathsf{T}}\), so \([1,-1]^{\mathsf{T}}\) is an eigenvector with eigenvalue \(1\).

Now write an arbitrary vector as a weighted sum of the two. Take \([2, 4]^{\mathsf{T}}\) and solve

\(\begin{bmatrix} 2 \\ 4 \end{bmatrix} = a \begin{bmatrix} 1 \\ 1 \end{bmatrix} + b \begin{bmatrix} 1 \\ -1 \end{bmatrix}\)(23)

This works under the condition, that the two eigenvectors must point in genuinely different directions and are not parallel.

Put the two pieces together. Decompose the starting pair as \(\hat{e}^{\,0} = a w_1 + b w_2\), then apply the matrix \(n\) times:

\(M(\kappa)^n \hat{e}^{\,0} = a\, \lambda_1^n\, w_1 + b\, \lambda_2^n\, w_2\)(24)

where \(w_1\) and \(w_2\) denote the two eigenvectors of \(M(\kappa)\), \(\lambda_1\) and \(\lambda_2\) denote their eigenvalues, and \(a\) and \(b\) denote the weights of the initial decomposition.

The weights and the eigenvectors are fixed once, at the start, and never change again. The only things on that right-hand side that depend on \(n\) are \(\lambda_1^n\) and \(\lambda_2^n\), two plain numbers raised to a power. So:

\(\text{stability} \iff |\lambda| \leq 1 \ \text{ for both eigenvalues, for every } \kappa\)(25)

Start from \(Aw = \lambda w\) and move everything to one side: \(Aw - \lambda w = 0\). To factor \(w\), the Eigenvector of \(A\) out, both terms have to be of the form "a matrix times \(w\)", and \(\lambda w\) is currently a number times \(w\). The fix is the identity matrix:

\(I = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}, \qquad I\begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} 1(x) + 0(y) \\ 0(x) + 1(y) \end{bmatrix} = \begin{bmatrix} x \\ y \end{bmatrix}\)(26)

Factoring is now allowed:

\((A - \lambda I)\, w = 0\)(27)

where \(A - \lambda I\) denotes the matrix \(A\) with \(\lambda\) subtracted from each of its two diagonal entries, \(w\) denotes the eigenvector and \(0\) denotes the zero vector.

Here's the intuition for why this works, before the general argument. Subtracting \(\lambda I\) never touches \(A\)'s own unique direction of the eigenvectors. It simply scales the grid, such that its eigenvectors stay exactly where they were:

\((A - \lambda I) w = (\lambda_i - \lambda)\, w\)(27.2)

for whichever \(w\) is an eigenvector of \(A\) with eigenvalue \(\lambda_i\). So pick \(\lambda\) equal to one of those actual eigenvalues, and that eigenvector's stretch factor hits exactly zero — it gets sent straight to the origin, not shrunk toward it. The other eigenvector gets turned down too, but by a different starting amount, so it lands somewhere other than zero and survives.

Subtracting \(\lambda\) off the diagonal doesn't touch \(A\)'s two eigen-directions at all — they stay exactly where they are. It only turns down how much stretching happens along each of them, and it turns both down by the same amount as \(\lambda\) climbs.

For the eigenvector itself, that stretch amount counts straight down as \(\lambda\) increases, and hits exactly zero the moment \(\lambda\) equals its own eigenvalue. At that exact point nothing is left riding along that direction, so the eigenvector doesn't just shrink toward the origin — it lands exactly on it.

The other eigen-direction gets turned down by the same \(\lambda\) too, but it started from a different eigenvalue, so it doesn't hit zero at the same time — it just ends up stretched by whatever's left, and everything the matrix touches ends up somewhere along that surviving direction instead.

That's the whole reason it collapses onto a line, and why that line is \(A\)'s other eigenvector. But this only tells us what happens once we already know which eigenvector and eigenvalue we're dealing with — which is exactly what we don't know yet. So we still need a way to test for this using only \(A\) itself.

What is being asked for is a nonzero \(w\) that this matrix sends to the zero vector. The only way for something nonzero to reach zero is for the whole plane to be squashed onto a line or onto a single point.

Squashing the plane onto a line crushes every area to zero, and the quantity that measures how a matrix scales areas is the determinant. So the condition is:

\(\det(A - \lambda I) = 0\)(28)

where \(\det\) denotes the determinant, that is the factor by which the matrix scales areas, and the equation says that factor is zero.

For a two-by-two matrix the determinant is:

\(\det \begin{bmatrix} p & q \\ r & s \end{bmatrix} = ps - qr\)(29)

Subtract \(\lambda\) from the diagonal of \(M(\kappa)\) and take the determinant. Abbreviate the bottom-left entry as \(x := \frac{\Delta t c^2}{\Delta x^2}g(\kappa)\), where \(x\) denotes that single entry, purely to keep the line readable:

\(\det \begin{bmatrix} 1 - \lambda & \Delta t \\ x & 1 - \lambda \end{bmatrix} = (1-\lambda)^2 - \Delta t \, x = 0\)(31)

where \((1-\lambda)^2\) denotes the product of the two diagonal entries and \(\Delta t \, x\) denotes the product of the two off-diagonal entries.

Now expand that product of off-diagonal entries, because it is where the interesting structure is hiding:

\(\Delta t \, x = \Delta t \cdot \frac{\Delta t \, c^2}{\Delta x^2} g(\kappa) = \frac{\Delta t^2 c^2}{\Delta x^2} g(\kappa)\)(32)

where \(\Delta t^2\) denotes the time step squared, produced by the two factors of \(\Delta t\) multiplying, \(c\) denotes the wave speed and \(\Delta x\) the grid spacing.

That fraction is a square. Define

\(C := \frac{c \, \Delta t}{\Delta x} \quad \Longrightarrow \quad C^2 = \frac{c^2 \Delta t^2}{\Delta x^2}\)(33)

where \(C\) denotes the Courant number, which brings us one step closer to the CFL condition. It is a dimensionless number that measures how far a wave travels in one time step, relative to the grid spacing.

So \(\Delta t \, x = C^2 g(\kappa)\) and the characteristic equation becomes

\((1 - \lambda)^2 = C^2 g(\kappa) \quad \Longrightarrow \quad 1 - \lambda = \pm\sqrt{C^2 g(\kappa)} \quad \Longrightarrow \quad \lambda = 1 \mp \sqrt{C^2 g(\kappa)}\)(34)

where \(\lambda\) denotes an eigenvalue of \(M(\kappa)\), \(C\) denotes the Courant number and \(g(\kappa)\) denotes the stencil factor. Since both signs occur, this is the same as \(\lambda = 1 \pm \sqrt{C^2 g(\kappa)}\).

Looking at what's under the root. \(C^2 \geq 0\) because it is a square, and \(g(\kappa) \leq 0\) always, as established at the top of this post. A non-negative number times a non-positive number is non-positive, so \(C^2 g(\kappa) \leq 0\) for every choice of parameters and every wavenumber.

No real number squares to a negative number, since a positive times a positive and a negative times a negative are both positive. That is exactly the gap the imaginary unit fills: \(j := \sqrt{-1}\), where \(j\) denotes the number whose square is \(-1\).

Concretely: \(\sqrt{-4} = \sqrt{(-1)(4)} = \sqrt{-1}\sqrt{4} = 2j\).

The same move symbolically. Write \(g(\kappa) = -|g(\kappa)|\), where \(|g(\kappa)|\) denotes the size of the stencil factor with its sign stripped off, a non-negative number between \(0\) and \(4\). Then:

\(\sqrt{C^2 g(\kappa)} = \sqrt{-C^2 |g(\kappa)|} = j\sqrt{C^2 |g(\kappa)|}\)(35)

where everything under the final root is non-negative, so that root is an ordinary real number. The two eigenvalues are therefore

\(\lambda_1 = 1 + j\sqrt{C^2|g(\kappa)|}, \qquad \lambda_2 = 1 - j\sqrt{C^2|g(\kappa)|}\)(36)

where \(\lambda_1\) and \(\lambda_2\) denote the two eigenvalues, \(1\) denotes their common real part and the root term denotes the size of their imaginary parts.

Both have real part exactly \(1\). With \(C = 1\) and \(\kappa = \pi\), where \(|g(\pi)| = 4\), the root is \(\sqrt{1 \cdot 4} = 2\), so \(\lambda = 1 \pm 2j\).

In Formula 24 I have established, that error gets multiplied by an eigenvalue at every step, so after \(n\) steps it has been multiplied by that eigenvalue raised to the \(n\)-th power. A number larger than 1, raised to a large power, grows without bound. A number of size 1 or smaller, raised to a large power, stays bounded. This is why the stability criterion \(|\lambda| \leq 1\) is stated in terms of size.

Here \(\lambda\) is complex rather than real, so "size" requires a definition that still applies. A complex number \(a+bj\) can be treated as a point in a plane, with \(a\) as the coordinate along one axis and \(b\) as the coordinate along the other. Its size is defined as the distance of that point from the origin, given by the Pythagorean formula:

\(|a + bj| = \sqrt{a^2 + b^2}\)(37)

This distance is the quantity to raise to the \(n\)-th power.

Both eigenvalues have real part \(1\) and imaginary part \(\sqrt{C^2|g(\kappa)|}\), established above. Squaring a square root returns what was under it, so the squared magnitude is:

\(|\lambda|^2 = 1^2 + \left(\sqrt{C^2|g(\kappa)|}\right)^2 = 1 + C^2|g(\kappa)|\)(38)

where \(|\lambda|^2\) is the squared magnitude of either eigenvalue, \(C\) is the Courant number and \(|g(\kappa)|\) is the size of the stencil factor.

\(C^2\) is a square, so it cannot be negative. \(|g(\kappa)|\) is an absolute value, so it cannot be negative either. Their product is therefore never negative, which means \(|\lambda|^2 \geq 1\) for every parameter choice and every wavenumber. Equality holds only when \(C = 0\) (no time-stepping at all) or when \(g(\kappa) = 0\), the flat mode \(\kappa = 0\) that the stencil cannot detect.

This scheme is unconditionally unstable. \(|\lambda|\) is never below \(1\), so error never shrinks — at best it holds steady, and generally it grows. Shrinking \(\Delta t\) shrinks \(C\), which pushes \(|\lambda|\) back toward \(1\) and slows the blow-up per step — but never below \(1\). And a smaller \(\Delta t\) means more steps are needed to reach the same simulated time, so the slower per-step growth is offset by a larger number of steps. The instability never actually disappears, at any time step.

The mechanism is the update ordering. Both lines of the scheme read from the same frozen snapshot: the velocity update uses the old displacement, and the displacement update uses the old velocity. Neither ever sees the other's freshly computed value.

For an oscillating system, each update uses a value that is already slightly out of date by the time it is applied, and because the state is going round in a loop rather than settling toward anything, being consistently one step behind pushes the state slightly outward each time round.

The fix is to update the velocity first and then feed the freshly computed velocity into the displacement update:

\(V_i^{n+1} = V_i^n + \Delta t \, c^2 \, \frac{U_{i+1}^n - 2U_i^n + U_{i-1}^n}{\Delta x^2}\)(39) \(U_i^{n+1} = U_i^n + \Delta t \, V_i^{n+1}\)(40)

where every symbol means what it did before, and the only change in the entire scheme is that the second line now reads \(V_i^{n+1}\), the velocity at the next time step, instead of \(V_i^n\). This variant is called symplectic Euler.

The velocity line is untouched, so its error recurrence carries over word for word. Redoing the subtraction on the displacement line gives:

\(e_i^{\,u,n+1} = e_i^{\,u,n} + \Delta t \, e_i^{\,v,n+1}\)(41)

where \(e_i^{\,v,n+1}\) denotes the velocity error at the next step and the other symbols are as before.

That right-hand side refers to a quantity at step \(n+1\), which is the thing being solved for, so nothing can be put into matrix form yet. It is not a real obstacle, because the velocity recurrence already expresses \(e_i^{\,v,n+1}\) entirely in terms of step-\(n\) quantities. Substitute it in:

\(e_i^{\,u,n+1} = e_i^{\,u,n} + \Delta t \left( e_i^{\,v,n} + \Delta t \, c^2 \frac{e_{i+1}^{\,u,n} - 2e_i^{\,u,n} + e_{i-1}^{\,u,n}}{\Delta x^2} \right)\)(42)

Distribute the outer \(\Delta t\) across both terms in the bracket, and combine it with the \(\Delta t\) already inside the second one:

\(e_i^{\,u,n+1} = e_i^{\,u,n} + \Delta t \, e_i^{\,v,n} + \Delta t^2 c^2 \frac{e_{i+1}^{\,u,n} - 2e_i^{\,u,n} + e_{i-1}^{\,u,n}}{\Delta x^2}\)(43)

Something appeared that was not there before: a spatial stencil now sits inside the displacement error equation, which previously contained no spatial term at all. That extra term is the essence of using the updated velocity, and it is the entire difference between the two schemes.

Apply the ansatz, factor out \(T_i\) and divide it out on the same justification as before. The stencil turns into \(g(\kappa)\) and the coefficient \(\frac{\Delta t^2 c^2}{\Delta x^2}\) is recognizable as \(C^2\). Grouping the two terms that multiply \(\hat{e}^{\,u,n}\):

\(\hat{e}^{\,u,n+1} = \left(1 + C^2 g(\kappa)\right)\hat{e}^{\,u,n} + \Delta t \, \hat{e}^{\,v,n}\)(44) \(\hat{e}^{\,v,n+1} = \frac{\Delta t \, c^2}{\Delta x^2} g(\kappa)\, \hat{e}^{\,u,n} + \hat{e}^{\,v,n}\)(45)

where \(1 + C^2 g(\kappa)\) denotes the new coefficient of the displacement amplitude, with the \(1\) coming from the displacement carrying itself forward and the \(C^2 g(\kappa)\) coming from the new stencil term, and all other symbols are as before.

Reading off the four entries the same way as before:

\(M_{\text{sympl}}(\kappa) = \begin{bmatrix} 1 + C^2 g(\kappa) & \Delta t \\ \frac{\Delta t \, c^2}{\Delta x^2} g(\kappa) & 1 \end{bmatrix}\)(46)

where \(M_{\text{sympl}}(\kappa)\) denotes the amplification matrix of the symplectic scheme for wavenumber \(\kappa\). Exactly one entry differs from the plain scheme: the top-left, which went from \(1\) to \(1 + C^2 g(\kappa)\).

Subtract \(\lambda\) from the diagonal of the symplectic amplification matrix and set the determinant to zero, writing \(g\) for \(g(\kappa)\) to keep things short:

\(\left(1 + C^2 g - \lambda\right)\left(1 - \lambda\right) - C^2 g = 0\)(47)

Multiplying out the two brackets and subtracting the trailing term, the \(C^2 g\) piece cancels completely and the two leftover \(-\lambda\) terms combine into one. What's left is:

\(\lambda^2 - \left(2 + C^2 g(\kappa)\right)\lambda + 1 = 0\)(49)

Much cleaner than the setup suggested: every parameter is buried inside a single middle coefficient, and the constant term is just \(1\). Name that coefficient:

\(B := 2 + C^2 g(\kappa)\)(50)

and the whole thing is now an ordinary quadratic, \(\lambda^2 - B\lambda + 1 = 0\), so the quadratic formula hands over the eigenvalues directly:

\(\lambda = \frac{B \pm \sqrt{B^2 - 4}}{2}\)(51)

Whether those come out as two real numbers or a complex pair depends entirely on the sign of \(B^2 - 4\).

Read \(B\) as an instruction: start at \(2\), then add \(C^2 g(\kappa)\). Now, \(C^2\) is a square and therefore never negative, and \(g(\kappa)\) is never positive. A non-negative number times a non-positive number is non-positive. So the thing being added to \(2\) can only ever pull \(B\) downward. It can never push it up.

So \(B \leq 2\), always, with equality only in the two degenerate cases \(C = 0\) and \(\kappa = 0\). How far down it actually gets pulled depends on two things: how large \(C\) is, and how negative \(g(\kappa)\) is for the mode in question.

If the eigenvalues come out complex, their magnitude will turn out (next section) to be exactly \(1\). If they come out real, one of them will turn out to be larger than \(1\) in size, which means growing instability. So checking whether \(B^2 - 4\) is negative or positive is really checking whether the scheme is stable at all, for this \(\lambda\) and \(\kappa\).

For the eigenvalues to be real, the discriminant has to be non-negative, meaning \(B^2 \geq 4\). Squaring anything strictly between \(-2\) and \(2\) gives less than \(4\), so \(B^2 \geq 4\) requires \(B\) to be at least \(2\) or at most \(-2\).

We have just established that \(B\) cannot exceed \(2\). So the only live route to real eigenvalues is \(B\) being dragged below \(-2\). The stability question has now become one concrete question: how large does \(C\) have to be to drag \(B\) past \(-2\)?

Take the mode that gets pulled down the hardest, which is the checkerboard \(\kappa = \pi\) with \(g(\pi) = -4\). If this mode stays safe, all the others do too, since every other \(\kappa\) has a smaller \(|g(\kappa)|\) and is therefore pulled down less.

For that mode the downward pull is \(C^2 \cdot 4\). And \(B\) starts at \(2\) and has to reach \(-2\), which is a drop of \(4\). Setting the pull equal to the available drop:

\(4C^2 = 4 \quad \Longrightarrow \quad C^2 = 1 \quad \Longrightarrow \quad C = 1\)(52)

where the left-hand side denotes the downward pull on \(B\) at the worst wavenumber, the right-hand side denotes the distance from \(2\) down to \(-2\), and the negative root is discarded because \(C\) is built from positive quantities.

So the three regimes are:

• \(C < 1\): the pull is less than \(4\), \(B\) stays inside \([-2, 2]\), the discriminant is negative, the eigenvalues are a complex pair.
• \(C = 1\): \(B\) lands exactly on \(-2\), the discriminant is exactly zero.
• \(C > 1\): \(B\) is dragged past \(-2\), the discriminant is positive, the eigenvalues are two real numbers.

For \(C \leq 1\), the discriminant is negative, so the two eigenvalues are complex:

\(\lambda = \frac{B}{2} \pm j\frac{\sqrt{4 - B^2}}{2}\)(53)

with real part \(B/2\) and imaginary part \(\sqrt{4-B^2}/2\).

Square both parts and add them, the way magnitude works for any complex number:

\(|\lambda|^2 = \left(\frac{B}{2}\right)^2 + \left(\frac{\sqrt{4 - B^2}}{2}\right)^2 = \frac{B^2 + (4-B^2)}{4} = 1\)(54)

The \(B^2\) in the real part and the \(B^2\) hiding inside the imaginary part cancel each other out exactly, no matter what \(B\) is. So \(|\lambda| = 1\) for every single \(B\) in this range, not just at the edges: every Courant number at or below \(1\), every wavenumber, both eigenvalues land exactly on the unit circle. Changing \(C\) or \(\kappa\) only slides them around that circle, never off it.

Here is why that matters physically. For the heat equation, a stable scheme had \(|G(\kappa)| < 1\) strictly — error actively shrank at every step, because diffusion itself destroys structure and a good scheme mimics that.

Here, multiplying by a complex number of magnitude exactly \(1\) is a pure rotation: nothing grows, nothing shrinks, only the phase turns. That matches the physics too — a wave conserves its energy rather than dissipating it, so a scheme that damped the error would also be damping the actual solution, which is just a different way of getting the wrong answer.

Sitting exactly on that boundary, rather than safely inside it, is called marginal stability. For a wave equation, that's not a compromise — it's the target.

Once \(C > 1\), the eigenvalues are real — and the instability can be read off without even solving for them.

For any quadratic \(\lambda^2 - B\lambda + 1 = 0\), the product of its two roots always equals the constant term — a standard shortcut, and here that constant is \(1\):

\(\lambda_1 \lambda_2 = 1\)(55)

This holds no matter whether the roots are real or complex.

Now, if both roots are real and multiply to \(1\), the only way for both to also have magnitude \(1\) is if they're exactly \(\pm 1\) — which happens at exactly one point, \(B = -2\). Move past that point and they can't both stay at magnitude \(1\) anymore, but the product still has to equal \(1\). The only way to keep multiplying to \(1\) is for one root to grow past \(1\) while the other shrinks below it — they become reciprocals of each other. That larger root is the unstable mode, and it grows geometrically from the moment it exists.

This transition isn't gradual. Below \(C = 1\), both eigenvalues slide along the unit circle as \(C\) and \(\kappa\) change. At \(C = 1\) exactly, they meet at the point \(-1\). Push \(C\) past \(1\) and they split apart along the real axis — one moving inward, one moving outward — and the outward one is already outside the circle. There is no "barely unstable" in between: the instant an eigenvalue leaves the circle, it is gone.

The first mode to fail is the checkerboard, \(\kappa = \pi\) — same as in the heat equation — because it has the largest \(|g(\kappa)|\) and so is the first one pulled past the danger threshold. The sharpest content the grid can represent is always the first casualty.

Unpack the result. \(C \leq 1\) is \(\frac{c \Delta t}{\Delta x} \leq 1\), which rearranges to \(c\Delta t \leq \Delta x\): the wave may not travel more than one grid cell per time step. That's the same domain-of-dependence topic I started with in the last post, which cleanly closes the loop.

That intuitive argument and the algebra above land on the same number, but they aren't doing the same job. Domain-of-dependence only says the scheme must fail past \(C=1\) — it doesn't say everything below \(C=1\) is safe. Plain forward Euler has the exact same domain of dependence as the symplectic scheme and is unstable at every Courant number, which proves the intuition alone isn't enough. Von Neumann analysis is what actually nails the boundary down, and it's what told us which mode fails first, that the failure is sudden rather than gradual, and that update ordering matters as much as step size.

Side by side, the two stability bounds from this series:

\(\text{heat:} \quad \Delta t \leq \frac{\Delta x^2}{2\alpha} \qquad\qquad \text{wave:} \quad \Delta t \leq \frac{\Delta x}{c}\)(56)

where \(\alpha\) is the thermal diffusivity, \(c\) the wave speed, and \(\Delta x\), \(\Delta t\) the grid spacing and largest allowed time step.

The heat bound is quadratic in \(\Delta x\), the wave bound is linear — and that difference is exactly why these solvers get expensive as you refine the grid. Halve \(\Delta x\) to double your spatial resolution, and you don't just double the number of grid points: CFL also forces you to shrink \(\Delta t\), so you need more time steps to reach the same final time. For the wave equation that's a further \(2\times\), for a combined \(4\times\) more work. That's the real reason CFL-bound solvers are expensive: refining resolution taxes you twice — once for the extra grid points, once for the smaller time step the CFL condition then forces on you.