Software 13.5.1 (Stability of FTCS). FTCS for \(u_t = u_{xx}\) with \(u(x, 0) = \sin \pi x\) plus a sawtooth of size \(10^{-12}\), on \(50\) intervals, up to \(t = 0.1\): the largest error for \(r = 0.45\), \(0.5\), \(0.55\); then the error for \(r = 0.4\) on four grids:

>>> import numpy as np
>>> def ftcs(N, r, T):              # u_t = u_xx on (0, 1), u = 0 at the ends
...     x = np.linspace(0, 1, N + 1); dx = 1 / N; dt = r * dx**2
...     u = np.sin(np.pi * x) + 1e-12 * (-1)**np.arange(N + 1)  # tiny noise
...     u[0] = u[-1] = 0
...     steps = int(round(T / dt))
...     for _ in range(steps):
...         u[1:-1] = u[1:-1] + r * (u[2:] - 2*u[1:-1] + u[:-2])
...     exact = np.exp(-np.pi**2 * steps * dt) * np.sin(np.pi * x)
...     return np.abs(u - exact).max()
...
>>> for r in (0.45, 0.50, 0.55):           # N = 50, up to t = 0.1
...     print(r, f"{ftcs(50, r, 0.1):.3e}")
...
0.45 2.059e-04
0.5 2.422e-04
0.55 5.947e+23
>>> [f"{ftcs(N, 0.4, 0.1):.3e}" for N in (10, 20, 40, 80)]   # second order
['4.294e-03', '1.063e-03', '2.649e-04', '6.620e-05']

Reading the output. For \(r \leq \frac 12\) the error is \(2 \times 10^{-4}\). At \(r = 0.55\) the sawtooth is multiplied by \(|1 - 4r| = 1.2\) at each of \(455\) steps, a factor of \(10^{36}\), and the \(10^{-12}\) has grown to \(6 \times 10^{23}\). With \(r\) fixed at \(0.4\), each halving of \(\Delta x\) divides the error by \(4\): second order in \(\Delta x\), as Theorem 13.2.1 predicts, since \(\Delta t \propto \Delta x^2\).

MATLAB/Octave: u(2:end-1) = u(2:end-1) + r*diff(u,2) in a loop.

Software 13.5.2 (Amplification factors). The largest \(|G(\theta )|\) over \(0 \leq \theta \leq \pi \), and \(G\) at the sawtooth mode \(\theta = \pi \), for the schemes of this section:

>>> th = np.linspace(0, np.pi, 2001); s = np.sin(th / 2)**2
>>> G = {                                  # amplification factors G(theta)
...   'FTCS heat':    lambda r: 1 - 4*r*s,
...   'backward':     lambda r: 1 / (1 + 4*r*s),
...   'Crank-Nic.':   lambda r: (1 - 2*r*s) / (1 + 2*r*s),
...   'upwind':       lambda v: 1 - v * (1 - np.exp(-1j*th)),
...   'FTCS advect':  lambda v: 1 - 1j * v * np.sin(th),
...   'Lax-Wendroff': lambda v: 1 - 1j*v*np.sin(th) - 2*v**2*s}
>>> for name, p in (('FTCS heat', 0.4), ('FTCS heat', 0.6), ('backward', 5),
...                 ('Crank-Nic.', 5), ('upwind', 0.8), ('upwind', 1.2),
...                 ('FTCS advect', 0.5), ('Lax-Wendroff', 0.8)):
...     g = G[name](p)                     # worst case, and the highest mode
...     print(f"{name:13s} {p:4}  max|G| = {np.abs(g).max():.4f}"
...           f"  G(pi) = {complex(g[-1]).real:+.4f}")
...
FTCS heat      0.4  max|G| = 1.0000  G(pi) = -0.6000
FTCS heat      0.6  max|G| = 1.4000  G(pi) = -1.4000
backward         5  max|G| = 1.0000  G(pi) = +0.0476
Crank-Nic.       5  max|G| = 1.0000  G(pi) = -0.8182
upwind         0.8  max|G| = 1.0000  G(pi) = -0.6000
upwind         1.2  max|G| = 1.4000  G(pi) = -1.4000
FTCS advect    0.5  max|G| = 1.1180  G(pi) = +1.0000
Lax-Wendroff   0.8  max|G| = 1.0000  G(pi) = -0.2800

Reading the output. Every scheme has \(|G(0)| = 1\), so the maximum is at least \(1\); stability means it is not more. FTCS for heat fails at \(r = 0.6\) and upwind at \(\nu = 1.2\), both with \(G(\pi ) = -1.4\); centred FTCS for advection fails at \(\nu = 0.5\) (\(|G| = 1.118\) at \(\theta = \pi /2\)). Backward Euler and Crank–Nicolson are stable at \(r = 5\), but Crank–Nicolson leaves the sawtooth at \(-0.82\) per step.

MATLAB/Octave: max(abs(1 - 4*r*sin(th/2).^2)).

Software 13.5.3 (Upwind against Lax–Wendroff). A square pulse carried once round a periodic interval by \(u_t + u_x = 0\), on \(200\) cells with \(\nu = 0.8\):

>>> N, nu = 200, 0.8                       # u_t + u_x = 0, periodic on [0, 1)
>>> x = np.arange(N) / N
>>> u0 = ((x > 0.25) & (x < 0.5)).astype(float)          # a square pulse
>>> up, lw = u0.copy(), u0.copy()
>>> for _ in range(int(round(N / nu))):    # one full period
...     up = up - nu * (up - np.roll(up, 1))
...     lw = (lw - nu/2 * (np.roll(lw, -1) - np.roll(lw, 1))
...           + nu**2/2 * (np.roll(lw, -1) - 2*lw + np.roll(lw, 1)))
...
>>> for name, v in (('upwind', up), ('Lax-Wendroff', lw)):
...     print(f"{name:12s} max {v.max():.3f}  min {v.min():+.3f}"
...           f"  L1 error {np.abs(v - u0).mean():.4f}")
...
upwind       max 1.000  min +0.000  L1 error 0.0504
Lax-Wendroff max 1.195  min -0.195  L1 error 0.0347

Reading the output. Upwind keeps the solution between \(0\) and \(1\), but its numerical diffusion smears the edges, with mean error \(0.050\). Lax–Wendroff is more accurate overall (\(0.035\)) but overshoots to \(1.195\) and undershoots to \(-0.195\) next to the jumps: Godunov’s theorem in action.

MATLAB/Octave: circshift in place of np.roll.