Second order says that halving $h$ divides the error by four. Taken literally it promises any accuracy at all, given enough points. It does not deliver, and the two reasons it does not are the subject of this chapter: arithmetic is not exact, and the matrix gets worse as the mesh gets finer.

The condition number of the difference matrix#

The eigenvalues of the $(N-1) \times (N-1)$ matrix $\mathrm{tridiag}(-1,2,-1)$ are known in closed form,

so after the factor $1/h^2$ the extreme eigenvalues of $A$ are about $\pi^2$ and $4/h^2$, and

This is not a defect of the discretisation. The continuous operator $-d^2/dx^2$ is unbounded, and any consistent discretisation of it must have a condition number that grows as the mesh is refined. The practical consequence is a floor: solving in double precision, with unit roundoff $\varepsilon \approx 2.2 \times 10^{-16}$, cannot be expected to give a relative accuracy better than roughly $\varepsilon \kappa(A) \approx 10^{-16} N^2$. Setting that equal to the discretisation error $C h^2 = C N^{-2}$ gives

About half the digits, and no more, whatever the machine and however long one is willing to wait.

The prediction is quantitative and the lab confirms it. The error falls at exactly second order until $N \approx 10^5$, reaches a little under $10^{-8}$, and then climbs again. The climb is roundoff amplified by the condition number, and its detailed rate depends on how the rounding errors happen to fall, but the shape of the curve is completely reproducible: past the minimum, refining the mesh makes the answer worse and takes longer to do it.

Note also the cost. Each row of the table is a solve of $O(N)$ work, so the wall clock in the last column grows linearly, and the whole 20-fold improvement in accuracy from $N = 10^3$ to $N = 10^5$ costs a factor of 100 in time.

Getting a better rate without a better scheme#

The error expansion of the previous chapter says more than "the error is $O(h^2)$". For smooth $u$ the nodal error has an asymptotic expansion in even powers of $h$,

with coefficients independent of $h$. Two solves, at $h$ and $h/2$, then give two equations at the points the meshes share, and the combination

cancels $c_2$ exactly and leaves $O(h^4)$. This is Richardson extrapolation. It costs one extra solve, which for a tridiagonal system is a factor of three in total work, and it needs no new discretisation, no new stencil, and no new code.

On the smooth solution the observed order goes from 2 to 4, and at $N = 512$ the extrapolated answer is about five orders of magnitude better than the one it was built from. Switch to lowreg and the gain in order disappears entirely: both curves converge at about $1.5$, because $u''''$ does not exist at $x = 1/2$ and the expansion in powers of $h^2$ that justifies the combination is not valid there. The extrapolation still improves the constant, by a factor of about seven, which is a fair return for one extra solve but is not what it was asked to do.

That is the honest summary of extrapolation. It is nearly free when it works, it is easy to check whether it worked, since the observed order tells you, and when the assumed expansion is absent it degenerates into a mildly better constant. The next chapter takes the other route: change the discretisation itself.