088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 1The sine transform solved the five-point system in $O(N^2\log N)$, but it needed a rectangle, a uniform grid and constant coefficients. Without those, and with $N^2$ unknowns too many to eliminate, the remaining option is to iterate: start from a guess and improve it, using only the ability to apply the operator.
3This chapter is about how much difference the choice of iteration makes, and about a question that only appears once one iterates: when to stop.
5### Five iterations
7**Jacobi.** Replace every value by the average of its four old neighbours. It is the simplest thing that could work, and it converges: the eigenvalues of the iteration matrix are $\frac{1}{2}(\cos\pi k h + \cos \pi l h)$, whose largest magnitude is $1 - \pi^2h^2/2 + O(h^4)$. Reducing the error by a fixed factor therefore takes $O(N^2)$ sweeps, and since a sweep costs $O(N^2)$, the total is $O(N^4)$: worse than the dense solve is not, but not by much.
9**Gauss-Seidel.** Use each new value as soon as it is computed. No extra storage, no extra work, and the asymptotic rate is exactly the square of Jacobi's, which halves the number of sweeps. Still $O(N^2)$ sweeps.
11**SOR.** Overshoot the Gauss-Seidel update by a factor $\omega$. For this matrix the optimal choice is known,
13$$\omega_{\text{opt}} = \frac{2}{1 + \sin \pi h}, \qquad \rho = \omega_{\text{opt}} - 1 = 1 - 2\pi h + O(h^2),$$
15which turns $O(N^2)$ sweeps into $O(N)$. That is a genuine change of order, and it is also the method's weakness: $\omega_{\text{opt}}$ is known here only because the spectrum of this particular matrix is known in closed form, and the rate degrades quickly if $\omega$ is wrong in the wrong direction.
17**Conjugate gradients.** The matrix is symmetric positive definite, so CG applies. It needs one operator application and a few vector operations per step, and its convergence depends on $\sqrt{\kappa} = O(N)$ rather than on $\kappa$, giving $O(N)$ iterations like optimal SOR, but with no parameter to choose and no knowledge of the spectrum required. With a good preconditioner it does better; the labs here use none.
19**Multigrid.** Gauss-Seidel is a bad solver but a good *smoother*: two sweeps remove the oscillatory part of the error and leave something smooth. Smooth functions are well represented on a grid twice as coarse, where the remaining problem costs a quarter as much to work on, and the same argument applies again there. The resulting V-cycle costs $O(N^2)$, proportional to the number of unknowns, and reduces the error by a constant factor, independent of $N$, each time. It is the only method here whose iteration count does not grow with the mesh.
21<!--lab:history-->
23### Reading the three figures
25The first shows the residual. On $N = 64$ the classical relaxations have barely started after a hundred and twenty sweeps, SOR is doing what its rate predicts, CG is somewhat better than SOR, and multigrid has reached $10^{-12}$ in ten cycles, at which point it stops because it has hit the tolerance. Raise $N$ to 128 and everything except multigrid gets worse, exactly in proportion to the rates quoted above; multigrid does not notice.
27The second shows the true error, which is the quantity anybody actually cares about, together with the dashed line marking the error of the *exact* solution of the discrete system. Every curve levels off at that line and cannot go below it. That is the practical point of the chapter:
29> Iterating past the discretisation error is wasted work. It buys a more accurate solution of the wrong problem.
31A useful stopping rule follows. The discretisation error here is $O(h^2)$, so a relative residual tolerance of about $h^2$ is enough, and anything smaller is spent for nothing. The labs stop at a fixed relative residual so that the methods can be compared, and the second figure shows how much of that work each method wasted.
33The third figure puts the same errors against measured milliseconds, and adds the sine-transform solve of Chapter 6 as a single point. On the square, with a uniform grid, the direct solver wins outright and no iteration is worth considering. That is the honest comparison: the iterative methods in this chapter earn their place on the problems where the direct solver is unavailable, which is where the next two chapters go.
35### What carries over, and what does not
37Of these five, only three would survive a change of domain or of coefficients unchanged. Jacobi, Gauss-Seidel and CG need nothing but the ability to apply the operator, so they run on an L-shaped grid, on a variable-coefficient problem, on an unstructured mesh. SOR needs a spectrum it will not have in closed form. Multigrid needs a hierarchy of grids and transfer operators between them, which is straightforward on a rectangle, real work on a complicated domain, and the subject of a large literature (algebraic multigrid) when there is no geometry to hierarchise.
39The convergence *rates* also carry over, with the caveat that they were derived for this operator on this domain. A re-entrant corner, for instance, changes the constants without changing the orders, which the next chapter makes visible in a different way.