When FTCS Trembles, Upwind Goes Numb — First-Order Advection, Artificial Viscosity, and the CFL Trade

In the early 1960s, people next to von Neumann dropped centered differences into one-dimensional linear advection. The square pulse was supposed to drift. Instead the trailing edge grew oscillations and left the grid in a few steps.

The order is not the cause. Pulling the downstream value with equal weight breaks causality. A one-cell upwind stops the runaway and kills the corners. The boundary of that trade is the CFL number.

Move the CFL dial from 0.05 to 1.5 and, on the same square pulse, FTCS keeps trembling while upwind lies on the dashed line only near 1.

At 0.5 the upwind pulse spreads the most. Past 1.05, upwind and Lax–Wendroff both leave the frame. The stencil has become shorter than the characteristic.

The most natural difference is the first to fall#

Linear advection in one dimension is one line.

$q$ is the transported quantity and $u>0$ a constant speed. The solution is the initial shape translated at speed $u$.

Forward Euler in time and a centered second-order difference in space give this update.

$c \equiv u\,\Delta t / \Delta x$ is the Courant number, the number of cells the wave crosses in one step. The scheme is second-order in space and first-order in time, so the LTE (local truncation error) is $\mathcal{O}(\Delta t,\Delta x^2)$. Shrink $c$ as far as you like. This scheme still runs away.

Substitute $q^n_i = g^n e^{i k x_i}$ and the amplification factor is

so $|g|^2 = 1 + c^2 \sin^2(k\Delta x) > 1$. Every mode grows every step. One mode with $|g|$ above 1 is enough to wreck the field.

At the code's $c=0.8$, the mode with $\sin=1$ has size $\sqrt{1+c^2}=\sqrt{1.64}\approx 1.28$. That is 1.28 per step, and 400 steps is that factor to the 400th power. The rough $1.78^{400}$ quoted below is not this mode. The formula gives 1.28.

Information flows from one side only#

If $u>0$, then $q_i$ at $t_{n+1}$ may look only upstream, at $x<x_i$. Centered differencing pulls the downstream value $q_{i+1}^n$ with equal weight. It reads a side the flow has not reached.

Upwind fixes the direction.

The amplification factor is

and $|g|^2 = 1 - c(1-c)\cdot 4\sin^2(k\Delta x/2)$. For $0 \le c \le 1$, $|g|\le 1$.

The coefficient $c(1-c)$ peaks at 0.25 when $c=0.5$. At $c=0.8$ it is 0.16. On the shortest wave, $|g|^2 = 1 - 4\cdot 0.16 = 0.36$, so $|g|=0.6$. At $c=0.5$ that mode's $|g|$ falls to 0. At $c=1$ the subtraction vanishes and every mode has $|g|=1$. The smearing peaks at 0.5, not at the edge of the stable range.

Upwind = centered + a hidden viscosity#

Split the first-order upwind slope and the identity is

The front term is centered differencing. The back term is a discrete second derivative. So upwind is actually integrating

Advection was the target. The algorithm inserted a viscosity $\nu_{\text{num}} = u\Delta x/2$. Artificial viscosity, also called numerical diffusion.

That viscosity kills the runaway. The corners of the square pulse spread to a thickness $\sqrt{4\nu_{\text{num}} t}$. A smaller $\Delta x$ cuts the viscosity, but reaching the same time takes more steps, and the smear stays.

The later coefficient $u\Delta x(1-c)/2$ disagrees with this display because of $c$. The $u\Delta x/2$ above is the $c\to 0$ piece. At $c=0.8$, $(1-c)=0.2$, so the coefficient that remains is one fifth of that, $u\Delta x/10$. Only the form with $(1-c)$ hits zero at $c=1$.

CFL — when the grid can't keep up#

The value of $q_i$ at $t_{n+1}$ is whatever sat at $x_i - u\Delta t$ at time $t_n$. If that foot lies between $x_{i-1}$ and $x_i$, two grid values can interpolate it. Once $u\Delta t > \Delta x$, the foot is past $x_{i-2}$. The stencil never looks there.

That geometry is the Courant–Friedrichs–Lewy condition.

The numerical domain of dependence (the stencil) has to contain the physical one (the characteristic). CFL is necessary for every explicit scheme and never sufficient. FTCS can satisfy CFL and still detonate. The blast is broken causality, not a short stencil.

Three schemes in fifty lines of NumPy#

Same grid, same square pulse, same CFL.

import numpy as np
 
N = 200
dx = 1.0 / N
u = 1.0
cfl = 0.8
dt = cfl * dx / u
n_steps = 400
 
x = (np.arange(N) + 0.5) * dx
q0 = ((x > 0.2) & (x < 0.4)).astype(float)  # square pulse
 
def ftcs_advect(q, c):
    return q - 0.5 * c * (np.roll(q, -1) - np.roll(q, 1))
 
def upwind_advect(q, c):
    return q - c * (q - np.roll(q, 1))  # u > 0
 
def lax_wendroff_advect(q, c):
    flux_diff = 0.5 * c * (np.roll(q, -1) - np.roll(q, 1))
    visc = 0.5 * c * c * (np.roll(q, -1) - 2 * q + np.roll(q, 1))
    return q - flux_diff + visc
 
q_ftcs, q_up, q_lw = q0.copy(), q0.copy(), q0.copy()
for _ in range(n_steps):
    q_ftcs = ftcs_advect(q_ftcs, cfl)
    q_up = upwind_advect(q_up, cfl)
    q_lw = lax_wendroff_advect(q_lw, cfl)
 
# Exact solution: one full revolution back to start
q_exact = q0
print("FTCS  L1 err:", np.abs(q_ftcs - q_exact).sum() * dx)
print("Upwind L1 err:", np.abs(q_up - q_exact).sum() * dx)
print("LW    L1 err:", np.abs(q_lw - q_exact).sum() * dx)

At cfl=0.8 and 400 steps, the rough figures written with this post are:

• FTCS: the amplitude grows by about $1.78^{400}$ and becomes NaN or inf.
• Upwind: the square fattens into a Gaussian-like blob and the L1 error sits near 0.04.
• Lax–Wendroff: the amplitude holds, but ripples (Gibbs) trail the corners.

c=0.5 spreads upwind the most#

The peak of $c(1-c)$ at 0.5 is the smear on the screen. At $c=1$, $(1-c)$ is zero and first-order upwind shifts the square by one cell, intact. $c=0.8$ sits in between, at 0.16. Stable does not mean less numb. 0.8 is only less viscous than 0.5.

First order is too blunt, but second order has new headaches#

Upwind viscosity is welcome next to a shock. It flattens unphysical wiggles. The same term stays on in smooth regions, so the error is stuck at the size of a second derivative. Production codes pick one of two exits.

• Use a higher-order scheme (Lax–Wendroff, MUSCL, WENO) and add a flux limiter or a slope limiter that drops the order only near shocks. The calculation that pins down Lax–Wendroff ringing on a step with a limiter is that choice.
• Call a Riemann solver (Godunov, HLLC, Roe) on each face and integrate the causality directly.

Both are machines for deciding where viscosity is accepted and where it is refused. First-order upwind accepts it in every cell. Simple, safe, and the most diffuse.

Calling 400 steps one lap is a failure#

With $N=200$, $u=1$, and $c=0.8$, $\Delta t = c/N = 0.004$. One lap of a domain of length 1 takes $1/\Delta t = 250$ steps, which is $N/c$. Four hundred steps travel $400\times 0.004 = 1.6$. A pulse of width 0.2, shifted by 0.6, moves from $[0.2,0.4]$ to $[0.8,1.0]$. Those two squares do not overlap. The L1 of the sharp pair is 0.4.

The upwind L1 of 0.04 quoted above cannot be that comparison. Reading 0.04 as the smear after one clean lap requires 250 steps. First-order upwind in one dimension lives up to $c=1$. Once the thing being carried is an interface, why does the limit drop so far? Where interface advection is left with a CFL of 0.05 shows where stencil length and the interface constraint split.