e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 1Many of the steady states of classical physics satisfy one equation. The temperature in a conducting plate that has stopped changing, the electrostatic potential in a region free of charge, the velocity potential of an incompressible irrotational flow, and the displacement of a stretched membrane at rest are all solutions of Laplace's equation,
3$$\Delta u \;=\; \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \;=\; 0 .$$
5A function satisfying it is called *harmonic*. If the interior carries a source, a heat source, a charge density, a load on the membrane, the equation gains a right-hand side and becomes Poisson's equation, $-\Delta u = f$. The two differ only in that right-hand side, and every method in this tutorial applies to both, though not always at the same order of accuracy.
7The equation by itself does not determine $u$. Adding any harmonic function to a solution gives another solution, so a condition on the boundary $\partial\Omega$ of the region $\Omega$ is needed. Three kinds are standard. A **Dirichlet** condition prescribes $u = g$ on the boundary: the plate is held at a known temperature along its edge. A **Neumann** condition prescribes the outward normal derivative $\partial u/\partial n = g$: the flux through the edge is known instead. A **Robin**, or mixed, condition prescribes a combination of the two. We say the *Dirichlet problem* is the first of these:
9$$\Delta u = 0 \ \ \text{in } \Omega, \qquad u = g \ \ \text{on } \partial\Omega,$$
6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 11or the same with $-\Delta u = f$ inside. It is well posed: for a reasonable domain and continuous data a solution exists, it is unique, and it depends continuously on $g$. Neumann and mixed conditions are not treated anywhere in these chapters, and neither are variable coefficients, nonlinearity, or time dependence. What is left is the simplest elliptic boundary value problem there is, which is exactly why it is worth ten chapters and a closing benchmark: nearly every technique used on harder elliptic problems can be introduced here, where the answer is known and nothing else is in the way.
13<!--fig:dirichlet-problem-->
15### What is being measured
17A tutorial on numerical methods has to decide what counts as success. Showing that a method converges is not difficult and not very informative, since almost anything consistent and stable converges. The question a user actually faces is narrower: for the accuracy I need, on the problem I have, which method gets there in the least time, and under what conditions does that answer change.
19To ask that question honestly one needs an error, and an error needs a truth to measure against. Reference solutions computed at higher resolution by the same family of methods are the usual substitute, and they are least reliable exactly where the interesting behaviour is. Here we avoid them entirely. Every problem in this tutorial has an exact solution known in closed form, either manufactured by choosing $u$ first and defining the data from it, or taken from the supply of harmonic functions that can be written down. The error of every computation is therefore known at every point, to full precision, at every resolution.
6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 21What the chapters then compare is accuracy against cost, with the cost measured as time for one solve inside the numerical layer. The spread turns out to be large. Across the eleven chapters, methods that are reasonable choices for the same problem differ by roughly fourteen orders of magnitude in the accuracy they deliver for comparable work.
23<!--fig:accuracy-cost-->
25Two things drive that spread, and they recur in every chapter. The first is the order of the method, which sets the slope of the curve rather than its position, and therefore decides who wins only once the accuracy demand is serious. The second is structure: a method that exploits a special feature of the problem, a rectangle, a separable operator, an analytic solution, a smooth boundary, can be faster than a general method by factors that no amount of tuning recovers. The conditions attached to each such method are as much the subject of this tutorial as the speedups are.
27### The domains
29The geometry is what removes those special features, one at a time, so the chapters are organised by domain rather than by method.
31<!--fig:domains-->
33**Part I, one dimension.** Poisson's equation on $[0,1]$. In one dimension Laplace's equation is $u'' = 0$, whose solution is the straight line through the two boundary values, so there is nothing to compute and the model problem has to be Poisson's. It has the same conditioning and the same convergence theory as the two-dimensional problem at a hundredth of the cost, which makes it the right place to develop the tools: manufactured solutions and the three-point scheme (Chapter 1), the roundoff floor and Richardson extrapolation (Chapter 2), a compact fourth-order scheme (Chapter 3), and Chebyshev collocation (Chapter 4).
35**Part II, the unit square.** Laplace's equation, with harmonic ground truth. The five-point Laplacian and the cost of a dense solve (Chapter 5), the same system diagonalised by the sine transform (Chapter 6), a sixth-order compact stencil for the same price (Chapter 7), and what to do when no direct solver is available (Chapter 8).
6292f12A benchmark chapter, and offline benchmarks to compare it againstJeremy Magland 37**Part III, beyond the square.** An L-shaped domain, where a re-entrant corner puts a singularity into the solution and second order is lost (Chapter 9), and a smooth closed curve, where the unknown moves onto the boundary and the answer comes back at full double precision (Chapter 10). A closing benchmark then runs one representative of every family on a domain that has both difficulties at once, and its work-precision plot can also take in timings recorded outside the browser, in real MATLAB, Python or Julia (Chapter 11).
39### How the pages work
41Each chapter is prose with runnable sections dropped into it. Every method is a short MATLAB-syntax script that you can read, edit, and re-run, with sliders for the resolution and for which exact solution is being asked for. The code runs in your browser through [numbl](https://numbl.org), a MATLAB-to-JavaScript compiler: there is no server, nothing to install, and nothing is sent anywhere. Heavier sections wait for a Run button rather than starting on sight.
43The timings you see are your own machine's, measured around the solve and nothing else. Any absolute number quoted in the prose was measured at the command line, where a native linear algebra library is available; in the browser the same code is a few times slower. Nothing in the tutorial depends on the absolute numbers, only on the ratios.