Crank-Nicolson is second order in time, backward Euler first

The rod is $0\le x\le 1$ with both ends held at $u=0$. The PDE's solution is the Fourier sine series

since $x(1-x)$ has sine coefficients $8/({\pi }^3{k}^3)$ for odd $k$ and zero for even $k$, and each mode decays at its own rate $k^2{\pi }^2$.

The second derivative is replaced by the centered difference ${\delta }^2{u}_j=(u_{j-1}-2u_j+u_{j+1})/h^2$ on 199 interior points, $h=1/200$. The time steps are then the only error left to measure. The Python reference, posts/refs/crank-nicolson-is-second-order-in-time.py, steps the same equations with scipy's banded solver and computes the exact space-discretized solution with a matrix exponential (output beside it).

Two schemes, one family

Both schemes average the space difference between the old and the new time level, with weight $\theta$ on the new one:

$\theta =1$ is backward Euler; $\theta =1/2$ is Crank-Nicolson, which is the trapezoidal rule in time. The left side is a tridiagonal system, with $1+2\theta r$ on the diagonal and $-\theta r$ beside it, $r=\Delta t/h^2$. The Thomas algorithm solves it in one sweep down and one back up.

Why the order differs: backward Euler evaluates the right side at the end of the step, which makes a local error of order $\Delta t^2$ per step and $\Delta t$ after the $T/\Delta t$ steps. The trapezoidal rule evaluates it at both ends, the first-order error terms of the two ends cancel, and what remains is $\Delta t^3$ per step and $\Delta t^2$ in total.

The exact answers

The error of a time-stepping scheme is measured against the exact solution of the equations it steps, $d{u}_j/dt={\delta }^2{u}_j$, which differs from the PDE's solution by the space error. That solution is also a sine series. The grid vectors $\sin (k\pi x_j)$ are eigenvectors of ${\delta }^2$, with eigenvalues $-{\lambda }_k$, ${\lambda }_k=(4/h^2){\sin }^2(k\pi h/2)$, which tend to $k^2{\pi }^2$ as $h\to 0$. The second difference of a quadratic is exact, so $-{\delta }^2{u}_0=2$ at every grid point. Expanding the constant 2 in the eigenvectors then gives the coefficients of $u_0$ in closed form:

which tends to $8/({\pi }^3{k}^3)$ as $h\to 0$. At $t=0$ the sum returns $u_0$ to rounding (tested below). The Python reference computes the same vector as $e^{-At}{u}_0$ with scipy's expm and agrees to $7\times 10^{-14}$.

Halving the time step

The cell runs both schemes at $\Delta t=0.02,0.01,\dots ,0.000625$ (5 to 160 steps) and takes the largest error over the grid at $t=0.1$. The measured order between two runs is ${\log }_2$ of the ratio of their errors.

Backward Euler's order climbs from 0.94 to 0.996, toward 1. Crank-Nicolson's is 2.14, 2.09, 2.11, 2.23 and then 2.000. The excess over 2 comes from the fast modes. For a mode with $\Delta t{\lambda }_k>2$, Crank-Nicolson's factor per step, $(1-\Delta t{\lambda }_k/2)/(1+\Delta t{\lambda }_k/2)$, is negative and close to $-1$, so the scheme barely damps a mode the equation kills at once. Backward Euler's factor $1/(1+\Delta t{\lambda }_k)$ goes to 0 instead. At $\Delta t=0.02$ that is every mode from $k=5$ up, and it adds 47% (fast_share) to the error of the slowest mode alone. As $\Delta t$ shrinks those modes have decayed before they matter; at $\Delta t=0.000625$ the error is the slowest mode's to 1%, and its order is 2.000.

Against the PDE's own Fourier series the Crank-Nicolson error falls with the rest until it meets the space error, $1.95\times 10^{-6}$ (the horizontal line). The two errors have opposite signs, so at $\Delta t=0.00125$ they partly cancel and the total dips below both. After that the total is the space error. Refining $\Delta t$ further buys nothing unless $h$ is refined with it.