The five-point matrix is not a general sparse matrix that happens to have five diagonals. It is

and everything about $T$ is known in closed form. Its eigenvectors are the discrete sine vectors

for $k = 1, \dots, N-1$. A Kronecker sum has the products of the eigenvectors as its eigenvectors and the sums of the eigenvalues as its eigenvalues, so

The system is therefore diagonal in a basis we can write down. Expand the right-hand side in that basis, divide each coefficient by $\lambda_k + \lambda_l$, and transform back. The expansion is a two-dimensional discrete sine transform, and each one-dimensional sine transform is one FFT of the odd extension, so the whole solve costs $O(N^2 \log N)$ for $N^2$ unknowns, with no matrix formed and no storage beyond the grid.

This is the fast Poisson solver, and the important thing about it is that it is not an approximation. It returns the solution of the same linear system the dense solve returns, to within roundoff, in the same sense that a better route to the same place is not a different place.

Where the two methods overlap, the errors agree to every digit shown, which is the check that matters: the fast solver has not changed the discretisation, only the arithmetic used to invert it. Where they do not overlap, the dense solve is simply unavailable, and the sine transform is solving a million unknowns in about half a second.

What it buys, and what it costs#

The accuracy at fixed $N$ is unchanged, because the discretisation is unchanged. What changes is which $N$ is affordable, and with second-order accuracy that is worth about four digits: $N = 32$ gives $10^{-5}$, $N = 1024$ gives $10^{-8}$.

It also shows the limit of the whole approach. Second order means the error is $C N^{-2}$, so twelve digits would need $N \approx 10^6$, which is $10^{12}$ unknowns. No solver, however fast, rescues a second-order discretisation at that accuracy. If we want twelve digits we have to stop refining and change the stencil, which is Chapter 7.

The conditions#

The derivation used every one of the following, and each is a real restriction.

• A rectangle, so that the grid is a tensor product and the matrix is a Kronecker sum.
• A uniform grid, so that $T$ is the same in every line and its eigenvectors are sine vectors.
• Constant coefficients, so that the operator separates. A variable coefficient $\nabla\cdot(a\nabla u)$ destroys the structure at once.
• Dirichlet conditions on all four sides. Neumann or periodic conditions are equally fine, with the cosine or Fourier transform in place of the sine transform, but a mixture along one side is not.

Break any of them and this chapter's method disappears, while the previous chapter's method merely gets slower. That trade, between methods that exploit structure and methods that tolerate its absence, is the theme of the rest of the tutorial. Chapter 8 asks what to do when the structure is gone and the matrix is too large for elimination.