1From here on the equation is Laplace's,
3$$\Delta u = u_{xx} + u_{yy} = 0 \ \ \text{in } \Omega, \qquad u = g \ \ \text{on } \partial\Omega,$$
5and $\Omega$ is the unit square. Poisson's equation is no harder to discretise, and everything below applies to it unchanged, but with $f = 0$ the only error in the computation is the error of the method, which is what we are trying to measure.
7### Ground truth in two dimensions
9The manufactured-solution trick needs an adjustment. Choosing $u$ freely and setting $f := -\Delta u$ would give a Poisson problem, not a Laplace problem. Instead we choose $u$ from the supply of functions that are *already* harmonic, and take its boundary trace as the data. Four are available in these labs:
11| choice | $u$ | why it is here |
12|---|---|---|
13| `expsin` | $e^{x}\sin y$ | entire, gentle, the standard test |
14| `poly` | $\operatorname{Re}(x+iy)^k$ | a harmonic polynomial; exact for a scheme of order $k$ or higher |
15| `wave` | $e^{ky}\sin(kx)$ | oscillates along one edge and grows by $e^{k}$ across the square |
16| `log` | $\log|z - z_0|$, $z_0$ just outside | harmonic inside, but with a singularity a distance $d$ from the boundary |
18The real and imaginary parts of any analytic function are harmonic, which is where all four come from, and it is worth noticing how much freedom that gives: any solution one can write down in closed form is available as ground truth, including ones that are nearly singular.
20### The five-point stencil
22Put a uniform grid of spacing $h = 1/N$ on the square and replace both second derivatives by their three-point differences. The result is the five-point formula
24$$\frac{u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1} - 4u_{i,j}}{h^2} \;=\; 0,$$
26one equation for each of the $(N-1)^2$ interior points, with the boundary values moved to the right-hand side. Its truncation error is $\frac{h^2}{12}(u_{xxxx} + u_{yyyy}) + O(h^4)$, so the scheme is second-order accurate, and the argument of Chapter 1 carries over unchanged: the error satisfies the same discrete equation with the truncation error as its source, and therefore inherits its size.
28Written in terms of the array of interior unknowns, the matrix has the Kronecker form $A = I \otimes T + T \otimes I$ with $T = \mathrm{tridiag}(-1,2,-1)$. That structure will matter a great deal in the next chapter. In this one we ignore it and hand the assembled matrix to Gaussian elimination.
30<!--lab:solve-->
32The error field is worth looking at rather than summarising. It is smooth, largest in the interior, and zero on the boundary, because the boundary values are exact by construction. Its shape is again the solution of a Poisson problem with the truncation error as source. Try `poly` with $k = 3$: the fourth derivatives of a cubic vanish, the leading truncation error is identically zero, and the scheme is exact to roundoff. Then try `log` with the singularity close to the edge, where $u_{xxxx}$ is large near that edge and the error concentrates there.
34### The wall
36<!--lab:cost-->
38The convergence is second order, cleanly. The cost is the problem. Dense elimination on $(N-1)^2$ unknowns costs $O(N^6)$ operations and $O(N^4)$ storage, and the measured times follow that faithfully: each doubling of $N$ multiplies the time by about 64. The last row of the table is a 2209 by 2209 dense solve, which buys five and a half correct digits and takes a fraction of a second at the command line, several seconds in a browser. $N = 128$ would need about 30 GB and is not a question of patience.
40Two remarks on what this does and does not show.
42The matrix is extremely sparse, with at most five entries per row, and treating it as dense is a waste. A real sparse direct solver with a good ordering, nested dissection, would solve this system in $O(N^3)$ operations with $O(N^2 \log N)$ fill, which is a completely different curve. numbl's backslash converts a sparse matrix to dense before solving, so that option is not available in this page, and the naive version is what the lab measures.
44But even nested dissection is not the right answer here. The matrix is not merely sparse: it is the Kronecker sum of two identical tridiagonal matrices, whose eigenvectors are known in closed form. The next chapter uses that to solve the same system in $O(N^2 \log N)$ operations, with the same answer to the last bit.