A second-order method with a perfect solver still needs $10^{12}$ unknowns for twelve digits. The way out is the same as in Chapter 3: keep the structure that makes the solver fast, and change the stencil.

The compact nine-point stencil#

Add the four diagonal neighbours, with weights chosen to cancel as much of the truncation error as possible:

This is the Mehrstellen operator. For Poisson's equation the $h^2$ term is $\frac{h^2}{12}\Delta f$, which is known and can be subtracted, giving a fourth-order scheme with the same nine points. For Laplace's equation, $\Delta^2 u = \Delta(\Delta u) = 0$ identically, so that term is not merely computable, it is zero, and the next one vanishes as well: the scheme is sixth-order accurate.

The stencil is still a sum of tensor products of tridiagonal operators, so the sine vectors are still its eigenvectors. With $a = \pi k/N$ and $b = \pi l/N$ the eigenvalue is

and the solver of the previous chapter applies verbatim, with one line changed. A sixth-order method for the price of a second-order one, on the same grid, in the same time.

Chebyshev in two dimensions#

The one-dimensional collocation of Chapter 4 extends to the square by taking the tensor product of two Chebyshev grids. Written on the array $U$ of interior values, the discrete problem is a Sylvester equation,

where $A$ is the one-dimensional operator on interior nodes and $F$ collects the boundary contributions. Solving it by forming $I \otimes A + A \otimes I$ and calling backslash costs $O(N^6)$ and is the two-dimensional dense solve all over again. Diagonalising $A = V\Lambda V^{-1}$ once instead turns the equation into a division,

at a cost of one eigendecomposition and four matrix products, $O(N^3)$ in total. This is the matrix-decomposition method, and it is the same idea as the sine transform of Chapter 6, with the eigenvectors computed at run time instead of known in advance. Unlike the sine transform it does not care that the operator is a difference operator, only that the problem separates.

What the comparison says#

On expsin, at the default settings:

• the five-point scheme reaches $4 \times 10^{-8}$ with $10^{6}$ unknowns;
• the nine-point scheme reaches $3 \times 10^{-14}$ with $10^{3}$ unknowns, in under a millisecond;
• Chebyshev collocation reaches $10^{-14}$ with $225$ unknowns.

Three methods, the same problem, the same machine, and a factor of a million in unknowns between the first and the last. The nine-point curve is a straight line of slope $-6$ until it hits the roundoff floor near $10^{-14}$; the Chebyshev curve does not have a slope, it falls off a cliff.

The qualifications matter as much as the numbers.

• The sixth order is specific to Laplace's equation. Solve Poisson's equation with the same stencil and the order drops to four unless the source term is averaged too.
• Everything here needs the square. All three methods are tensor-product constructions. Chapter 9 removes one corner from the square and Chapter 10 replaces it with a curve, and none of these three survive that.
• The roundoff floor arrives earlier for higher order. Both high-order curves flatten near $10^{-13}$ and then creep upward as $N$ grows, exactly as in Chapters 2 and 4. There is no reason to run a sixth-order scheme past the resolution where it reaches the floor.
• Chebyshev's cost grows faster. $O(N^3)$ against $O(N^2\log N)$, and the eigendecomposition has a large constant. At the accuracies reachable here it never matters, because the required $N$ is so small, but on a problem needing $N = 500$ it would.