088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 1Second 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.
3### The condition number of the difference matrix
5The eigenvalues of the $(N-1) \times (N-1)$ matrix $\mathrm{tridiag}(-1,2,-1)$ are known in closed form,
7$$\lambda_k \;=\; 4\sin^2\!\Big(\frac{k\pi h}{2}\Big), \qquad k = 1, \dots, N-1,$$
9so after the factor $1/h^2$ the extreme eigenvalues of $A$ are about $\pi^2$ and $4/h^2$, and
11$$\kappa(A) \;=\; \frac{\lambda_{N-1}}{\lambda_1} \;\approx\; \frac{4}{\pi^2 h^2} \;\approx\; 0.4\,N^2 .$$
13This 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
15$$N_{\text{best}} \;\sim\; \varepsilon^{-1/4} \approx 10^4, \qquad \text{error}_{\text{best}} \;\sim\; \sqrt{\varepsilon} \approx 10^{-8}.$$
17About half the digits, and no more, whatever the machine and however long one is willing to wait.
19<!--lab:floor-->
21The 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.
23Note 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.
25### Getting a better rate without a better scheme
27The 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$,
29$$u_h(x) = u(x) + c_2(x)\,h^2 + c_4(x)\,h^4 + \cdots,$$
31with coefficients independent of $h$. Two solves, at $h$ and $h/2$, then give two equations at the points the meshes share, and the combination
33$$u_R \;=\; \frac{4\,u_{h/2} - u_h}{3}$$
35cancels $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.
37<!--lab:richardson-->
39On 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.
41That 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.