Every chapter so far has posed a problem shaped to show one method at its best. This closing chapter does the opposite. It fixes one problem that is unfair to everybody and runs the tutorial's four families of methods on it, so that the accuracy-against-cost question the introduction promised can finally be answered on level ground: same domain, same data, same error points, same clock.

The domain is the Pac-Man,

the unit disk with the quadrant $x > 0$, $y < 0$ removed. Its boundary is a smooth arc and two straight edges; the edges meet the arc in two right-angle corners, and each other in a re-entrant corner of interior angle $3\pi/2$ at the origin. Both difficulties of Part "Beyond the square" are therefore present at once: the boundary cuts through grid cells, as in Chapter 10, and the re-entrant corner carries the $r^{2/3}$ singularity of Chapter 9. The corner has the same angle as the L-shape's, at the same origin, so the exact solutions of Chapter 9 (lp_lexact) work here unchanged: the singular function $r^{2/3}\sin(2\theta/3)$, the smooth $\log|z - z_0|$, and their sum, which is the generic case. As always, the data is the exact solution's boundary trace, and the error is known everywhere.

Four methods, four handicaps#

Each method is the one you have already met, pushed off its home ground.

• Shortley-Weller + SOR. The five-point scheme, with the stencil legs shortened wherever a neighbour falls outside the arc, so the boundary condition is imposed on the arc itself rather than on a staircase. The price is symmetry: the shortened stencil is not a symmetric matrix, conjugate gradients no longer applies, and the solver here is SOR. On smooth data the scheme is cleanly second order, which is Shortley and Weller's achievement; on generic data the corner takes that away, exactly as in Chapter 9.
• The harmonic basis. The method of particular solutions, with the Chapter 9 basis: harmonic polynomials plus the corner functions $r^{2k/3}\sin(2k\theta/3)$, fitted by least squares at boundary points clustered towards the corner. Nothing about a curved arc troubles it, since the basis never sees the boundary's shape, only samples of the data on it.
• The double layer on graded panels. The Nyström method of Chapter 10 loses its trapezoid rule here, because the rule's spectral accuracy needed a smooth periodic integrand and the corners break both the geometry and the density: at each corner the density itself is singular, near the re-entrant one like $r^{-1/3}$. The replacement is composite 12-point Gauss-Legendre on panels refined dyadically towards the three corners, to a depth $d$ that is the method's resolution parameter.
• Chebyshev on the polar rectangle. In polar coordinates the Pac-Man is the rectangle $(r, \theta) \in (0,1) \times (0, 3\pi/2)$, Laplace's equation becomes $r^2 u_{rr} + r u_r + u_{\theta\theta} = 0$, and the matrix decomposition of Chapter 7 applies with two different one-dimensional operators. The whole side $r = 0$ of the rectangle is the corner, a single boundary point with a known value, so the coordinate degeneracy costs nothing to impose. What it costs is convergence: the $r^{2/3}$ singularity now sits on an edge of the expansion domain, and the geometric convergence of Chapters 4 and 7 drops to an algebraic rate.

The benchmark#

The sweep below runs all four methods over their own resolution ladders and measures the maximum error at a fixed lattice of interior points, kept a little away from the boundary (the reason is the close-evaluation problem of Chapter 10, which would otherwise charge the double layer for a difficulty every method of this family has and this chapter is not about) and reaching within $0.125$ of the corner. Timings are for the solve only, warmed up, as everywhere in this tutorial.

On the generic data the four curves tell one story with four morals. The difference scheme converges at the corner-limited order $4/3$ and buys its accuracy with by far the most unknowns: at $N = 128$ it holds about $38\,000$ of them and stands at $4 \times 10^{-4}$. The polar Chebyshev method does somewhat better than second order in this range, but it has lost its exponential rate to the corner; on the smooth data, where nothing is singular, it reaches $10^{-10}$ with $40$ points per direction and is the second-best method on the board, which is a useful reminder that the ranking is a property of the problem, not of the methods. The graded double layer converges geometrically in the depth $d$, but with the modest ratio $2^{2/3} \approx 1.6$ per unit of depth, which is the corner exponent speaking: the leftover error is what the unresolved $r^{-1/3}$ tail of the density costs, about $2^{-2d/3}$. And the harmonic basis, whose every element satisfies the equation and whose two singular families are both in the span, is exact to roundoff with 97 unknowns in a few milliseconds, ten to twelve orders of magnitude ahead of everything else at equal cost. The conditions attached to that verdict are the ones stated at the end of Chapter 9, and they have not weakened: the equation must be homogeneous and linear with known solutions, and the singularities must be known in advance.

Note that the first plot cannot rank the methods, because each refines a different parameter; only the plot against time can. That is the plot this tutorial has been building towards since the introduction.

Benchmarks beyond the browser#

Everything above ran in your browser through numbl, and the browser is a particular place to do linear algebra: the dense solves run on a TypeScript fallback rather than a native LAPACK, and the loops run through numbl's JIT. How much of the picture is the method, and how much the platform? To answer that, the repository carries offline ports of the two best methods, the harmonic basis and the graded double layer, for real MATLAB, Python with NumPy, and Julia:

python3 bench/python/bench_pacman.py
julia bench/julia/bench_pacman.jl
matlab -batch "run('bench/matlab/bench_pacman.m')"

Each writes a small JSON file of resolutions, errors and times; the "load offline results" button under the work-precision plot draws them into it as dashed curves. The MATLAB harness does not port anything: it runs this chapter's driver and script, unchanged, in real MATLAB, which is the most direct numbl comparison available. The Python and Julia scripts are line-by-line ports, and the check that they are honest ports is in the plot itself: their error curves must land on top of the browser's, since error depends only on the mathematics, while their time curves may not.

The reason for the upload, rather than a table of times measured on the authors' machine, is comparability. A timing from someone else's hardware shares no clock with the numbers your browser just produced; a timing from the scripts above, run on the machine you are reading this on, does. Loading a file recorded elsewhere is still allowed, and the legend will say what it was, but the one comparison the plot can defend is browser and native on the same machine. In general, on this problem, expect the native runs to shift the dashed curves left by a factor that is large for the dense-solve-heavy double layer and modest for the small least-squares fits of the harmonic basis; how large is precisely what the plot is for.