A benchmark chapter, and offline benchmarks to compare it against
Chapter 11 fixes one problem that is unfair to every method and runs the
tutorial's four families on it: the Pac-Man domain, the unit disk with a
quadrant removed, whose boundary carries a smooth arc, two right-angle
corners, and the re-entrant corner of Chapter 9. The exact solutions of the
L-shaped chapter apply unchanged, since the corner has the same angle at the
same origin.
The four methods are the earlier ones pushed off their home ground:
Shortley-Weller five-point differences with SOR, since the shortened stencils
are not symmetric and CG no longer applies; the harmonic basis of Chapter 9;
the double layer of Chapter 10 with the trapezoid rule replaced by
Gauss-Legendre panels graded towards the three corners; and Chebyshev
collocation on the polar rectangle, where the corner singularity sits on an
edge of the expansion domain and the geometric convergence of Chapters 4 and 7
becomes algebraic.
bench/ holds offline ports of the two best of those methods for real MATLAB,
Python and Julia. The MATLAB harness runs this chapter's own driver and script
unchanged; the other two are line-by-line ports, and the check that they are
faithful is that their error curves land on the browser's, since the error
depends only on the mathematics. Each writes a small JSON file that the
chapter's work-precision plot loads as dashed curves. The upload exists rather
than a table of the author's timings because a reader's browser numbers and
their own offline numbers share a clock, and someone else's hardware does not.
With several series then sharing a colour and differing only in dash, the
legend stops being enough to tell them apart, so the work-precision panel
gained two ways to ask: pointing at a marker names its series and reads out
its resolution, error and time, and pointing at a legend entry picks that
curve out of the rest. Plot reports the pointer in frame pixels for it and the
panel converts back through the frame it is already handed.
test:matlab gains nine assertions, one for every number the new chapter's
prose claims.
26 changed files+1938−66
.gitignoremodified+1−0View file
@@ -3,3 +3,4 @@ dist
33 .DS_Store
44 *.local
55 tmp
6+bench-pacman-*.json
CLAUDE.mdmodified+5−0View file
@@ -23,6 +23,11 @@ src/
2323 onedim/driver.m the driver for chapters 1-4, plus their shared panels
2424 twodim/driver.m the driver for chapters 5-7, plus their shared panels
2525 chNN/ per-chapter text.md, editable .m script, chapter.ts
26+ ch11/panels.tsx also the offline-results upload for the benchmark plot
27+bench/ offline ports of the two best methods for real MATLAB,
28+ Python and Julia; the JSON they write (schema
29+ "laplace-dirichlet-bench-1", documented in bench/README.md)
30+ loads into chapter 11's work-precision plot
2631 ```
2732
2833 A **chapter** is prose in one `text.md` with `<!--lab:id-->` markers where its labs
README.mdmodified+10−1View file
@@ -14,7 +14,7 @@ answer changes.
1414
1515 An introduction, before Chapter 1, defines Laplace's and Poisson's equations and the
1616 Dirichlet, Neumann and Robin boundary conditions, states what the tutorial measures,
17-and lays out the four domains the chapters work on. Every chapter also carries one or
17+and lays out the five domains the chapters work on. Every chapter also carries one or
1818 two static diagrams, drawn as inline SVG, so that a page says something before its
1919 labs are run.
2020
@@ -70,6 +70,11 @@ domains the chapters work on.
7070 10. **Boundary integral equations** — a smooth curved domain, the double-layer
7171 potential, a second-kind Fredholm equation, the spectral accuracy of the
7272 periodic trapezoid rule, and the close-evaluation problem that comes with it.
73+11. **Benchmarks: one problem, every method** — the Pac-Man domain, whose curved
74+ boundary and re-entrant corner are unfair to everybody; the four families of
75+ methods on one work-precision plot; and offline runs of the two best methods
76+ in real MATLAB, Python and Julia (see `bench/`), which the plot can load so
77+ that browser and native timings from the same machine share an axis.
7378
7479 A **shared MATLAB library** page lists the small function files every lab boots
7580 with.
@@ -101,6 +106,10 @@ npm run test:browser # the built app in headless Chrome
101106 npm run test:site # the deployed site, in headless Chrome
102107 ```
103108
109+The offline benchmarks for Chapter 11 live in `bench/`, with their own README;
110+they need only whichever of MATLAB, Python + NumPy, or Julia is at hand, and each
111+writes a JSON file that the benchmark chapter's plot loads.
112+
104113 `test:matlab` needs a numbl clone; set `NUMBL_DIR` if it is not at `~/src/numbl`,
105114 and build its CLI (`npm run build:cli` in that checkout) after updating it — the
106115 tests run `dist-cli/cli.js`, not the TypeScript sources.
bench/README.mdadded+60−0View file
@@ -0,0 +1,60 @@
1+# Offline benchmarks
2+
3+The benchmark chapter of the tutorial runs its four methods in the browser
4+through numbl. The scripts here run the two best of those methods, the
5+harmonic basis and the graded-panel double layer, outside the browser, on the
6+same Pac-Man problem, the same error lattice, and the same timing protocol
7+(solve only, warmed up, minimum over repeats when a solve is under 5 ms).
8+
9+Run whichever of these your machine has, from the repository root:
10+
11+```
12+python3 bench/python/bench_pacman.py # needs NumPy
13+julia bench/julia/bench_pacman.jl # standard library only
14+matlab -batch "run('bench/matlab/bench_pacman.m')"
15+```
16+
17+Each writes a `bench-pacman-*.json` file (the MATLAB one lands in the
18+repository root, the others in the current directory). Load those files from
19+the "load offline results" button under the benchmark chapter's
20+work-precision plot, and they appear as dashed curves beside the browser's.
21+
22+The point of running the scripts yourself, rather than reading times measured
23+on someone else's machine, is that your browser results and your offline
24+results then share a clock and can honestly share an axis.
25+
26+The MATLAB harness runs the chapter's own driver and editable script,
27+unchanged, in real MATLAB. The Python and Julia scripts are ports; the check
28+that they are faithful is that their error columns land on top of the
29+browser's curves, since the error depends only on the mathematics. All three
30+were verified to reproduce the browser's errors to plotting accuracy.
31+
32+## The file format
33+
34+```json
35+{
36+ "schema": "laplace-dirichlet-bench-1",
37+ "label": "what the legend should call this run",
38+ "kind": "both",
39+ "date": "2026-08-20",
40+ "methods": [
41+ {
42+ "name": "harmonic basis, least squares",
43+ "slot": "mps",
44+ "ns": [4, 8, 16],
45+ "E": [8.1e-4, 6.3e-6, 3.6e-10],
46+ "T": [0.2, 0.4, 0.9]
47+ }
48+ ]
49+}
50+```
51+
52+`slot` picks the plot colour and marker and is one of `fd2`, `fd4`,
53+`spectral`, `fast`, `iterative`, `mps`, `bie`; anything else falls back to
54+`fd2`. `E` is the maximum absolute error over the chapter's fixed interior
55+lattice, `T` the solve time in milliseconds, and entries may be `null` for a
56+resolution a method skipped. `kind` records which ground truth was measured
57+(`both`, `sing`, or `smooth`); the page warns if it differs from the one the
58+lab is currently showing. Nothing stops you from writing a file in this
59+format from any other language or machine; the legend will say whatever
60+`label` says.
bench/julia/bench_pacman.jladded+245−0View file
@@ -0,0 +1,245 @@
1+#!/usr/bin/env julia
2+# Offline benchmark of the Pac-Man problem in Julia.
3+#
4+# A port of the two best methods from the in-browser benchmark, the harmonic
5+# basis and the graded-panel double layer, on the same problem, the same error
6+# lattice, and the same timing protocol as the chapter driver. Run
7+#
8+# julia bench/julia/bench_pacman.jl
9+#
10+# from anywhere; it writes bench-pacman-julia.json into the current directory.
11+# Load that file from the benchmark chapter's "load offline results" button.
12+# Because it ran on the same machine as your browser, the timings share a
13+# clock and can share an axis. Uses only the standard library.
14+
15+using LinearAlgebra
16+using Dates
17+using Printf
18+
19+const KIND = "both"
20+const Z0 = 1.6 + 1.3im
21+
22+function exact(x, y)
23+ th = atan(y, x)
24+ th < -1e-13 && (th += 2pi)
25+ r = hypot(x, y)
26+ sing = r^(2 / 3) * sin(2 * th / 3)
27+ smooth = 0.5 * log((x - real(Z0))^2 + (y - imag(Z0))^2)
28+ KIND == "sing" && return sing
29+ KIND == "smooth" && return smooth
30+ return sing + smooth
31+end
32+exact(x::AbstractVector, y::AbstractVector) = [exact(x[k], y[k]) for k in eachindex(x)]
33+
34+"""The driver's error lattice: spacing 1/8, margins from the boundary."""
35+function lattice()
36+ px = Float64[]
37+ py = Float64[]
38+ for x in -1.0:0.125:1.0, y in -1.0:0.125:1.0
39+ r2 = x^2 + y^2
40+ keep = r2 <= 0.7225 + 1e-9 && r2 > 1e-12 &&
41+ !(x > 1e-12 && y < -1e-12) &&
42+ !(x > 1e-12 && y < 0.25 - 1e-9) &&
43+ !(y < -1e-12 && x > -0.25 + 1e-9)
44+ keep || continue
45+ push!(px, x)
46+ push!(py, y)
47+ end
48+ return px, py
49+end
50+
51+# ------------------------------------------- the method of particular sols --
52+
53+function pac_basis(x, y, M, K)
54+ np = length(x)
55+ A = zeros(np, 2M + 1 + K)
56+ z = complex.(x, y)
57+ col = 1
58+ p = ones(ComplexF64, np)
59+ for mdeg in 0:M
60+ A[:, col] = real.(p)
61+ col += 1
62+ if mdeg > 0
63+ A[:, col] = imag.(p)
64+ col += 1
65+ end
66+ p .*= z
67+ end
68+ th = [t < -1e-13 ? t + 2pi : t for t in atan.(y, x)]
69+ r = abs.(z)
70+ k, q = 0, 0
71+ while q < K
72+ k += 1
73+ k % 3 == 0 && continue
74+ q += 1
75+ A[:, col] = r .^ (2k / 3) .* sin.(2k .* th ./ 3)
76+ col += 1
77+ end
78+ return A
79+end
80+
81+function pac_bpoints(ntot)
82+ ne = max(6, round(Int, ntot / 6.7))
83+ na = max(12, ntot - 2ne)
84+ t = ((1:ne) .- 0.5) ./ ne
85+ c = t .^ 2
86+ s = 1.5pi .* ((1:na) .- 0.5) ./ na
87+ bx = vcat(c, cos.(s), zeros(ne))
88+ by = vcat(zeros(ne), sin.(s), -c)
89+ return bx, by
90+end
91+
92+function mps_solve(n, px, py)
93+ nb = 3n + 1
94+ bx, by = pac_bpoints(3nb)
95+ A = pac_basis(bx, by, n, n)
96+ c = A \ exact(bx, by) # pivoted QR least squares
97+ return pac_basis(px, py, n, n) * c
98+end
99+
100+# --------------------------------------------- boundary integral, panels ----
101+
102+"""Gauss-Legendre nodes and weights on [-1, 1], by Newton on P_p."""
103+function gauss(p)
104+ x = [cos(pi * (k - 0.25) / (p + 0.5)) for k in 1:p]
105+ dp = ones(p)
106+ for _ in 1:50
107+ p0 = ones(p)
108+ p1 = copy(x)
109+ for k in 2:p
110+ p2 = ((2k - 1) .* x .* p1 .- (k - 1) .* p0) ./ k
111+ p0 = p1
112+ p1 = p2
113+ end
114+ dp = p .* (x .* p1 .- p0) ./ (x .^ 2 .- 1)
115+ dx = p1 ./ dp
116+ x .-= dx
117+ maximum(abs, dx) < 1e-15 && break
118+ end
119+ w = 2 ./ ((1 .- x .^ 2) .* dp .^ 2)
120+ return x, w
121+end
122+
123+function graded_breaks(T, d, base)
124+ half = T / 2
125+ left = vcat(0.0, half .* 2.0 .^ (-d:0))
126+ bp0 = vcat(left, T .- reverse(left[1:end-1]))
127+ bp = [bp0[1]]
128+ for k in 2:length(bp0)
129+ a, b = bp0[k-1], bp0[k]
130+ ns = max(1, ceil(Int, (b - a) / base))
131+ for q in 1:ns
132+ push!(bp, a + (b - a) * q / ns)
133+ end
134+ end
135+ return bp
136+end
137+
138+function panel_nodes(bp, gx, gw)
139+ t = Float64[]
140+ w = Float64[]
141+ for k in 1:length(bp)-1
142+ a, b = bp[k], bp[k+1]
143+ append!(t, (a + b) / 2 .+ (b - a) / 2 .* gx)
144+ append!(w, (b - a) / 2 .* gw)
145+ end
146+ return t, w
147+end
148+
149+function pac_panels(d)
150+ gx, gw = gauss(12)
151+ be = graded_breaks(1.0, d, 0.35)
152+ ba = graded_breaks(1.5pi, d, 0.35)
153+ t1, w1 = panel_nodes(be, gx, gw)
154+ t2, w2 = panel_nodes(ba, gx, gw)
155+ t3, w3 = panel_nodes(be, gx, gw)
156+ n1, n3 = length(t1), length(t3)
157+ qx = vcat(t1, cos.(t2), zeros(n3))
158+ qy = vcat(zeros(n1), sin.(t2), t3 .- 1)
159+ qnx = vcat(zeros(n1), cos.(t2), ones(n3))
160+ qny = vcat(-ones(n1), sin.(t2), zeros(n3))
161+ qw = vcat(w1, w2, w3)
162+ qkap = vcat(zeros(n1), ones(length(t2)), zeros(n3))
163+ return qx, qy, qnx, qny, qw, qkap
164+end
165+
166+function bie_solve(d, px, py)
167+ qx, qy, qnx, qny, qw, qkap = pac_panels(d)
168+ n = length(qx)
169+ A = zeros(n, n)
170+ for j in 1:n, i in 1:n
171+ if i == j
172+ A[i, j] = -0.5 - qkap[j] / (4pi) * qw[j]
173+ else
174+ dx = qx[i] - qx[j]
175+ dy = qy[i] - qy[j]
176+ A[i, j] = (qnx[j] * dx + qny[j] * dy) / (2pi * (dx^2 + dy^2)) * qw[j]
177+ end
178+ end
179+ sigma = A \ exact(qx, qy)
180+ v = zeros(length(px))
181+ for k in eachindex(px)
182+ s = 0.0
183+ for j in 1:n
184+ dx = px[k] - qx[j]
185+ dy = py[k] - qy[j]
186+ s += (qnx[j] * dx + qny[j] * dy) / (2pi * (dx^2 + dy^2)) * qw[j] * sigma[j]
187+ end
188+ v[k] = s
189+ end
190+ return v
191+end
192+
193+# ------------------------------------------------------------- the sweep ----
194+
195+function sweep(solve, ns, px, py, ve)
196+ E = Float64[]
197+ T = Float64[]
198+ solve(ns[1], px[1:1], py[1:1]) # warm-up, which in Julia is compilation
199+ for n in ns
200+ el = @elapsed v = solve(n, px, py)
201+ # a second timed run always: stray compilation and GC pauses
202+ # otherwise land in single measurements
203+ el = min(el, @elapsed v = solve(n, px, py))
204+ rep = 0
205+ while el < 0.005 && rep < 4
206+ rep += 1
207+ el = min(el, @elapsed v = solve(n, px, py))
208+ end
209+ push!(E, maximum(abs.(v .- ve)))
210+ push!(T, el * 1000)
211+ end
212+ return E, T
213+end
214+
215+function json_num(v)
216+ isfinite(v) || return "null"
217+ return string(v)
218+end
219+json_arr(v) = "[" * join(json_num.(float.(v)), ",") * "]"
220+
221+function main()
222+ px, py = lattice()
223+ ve = exact(px, py)
224+ entries = String[]
225+ for (name, slot, solve, ns) in [
226+ ("harmonic basis, least squares", "mps", mps_solve, [4, 8, 16, 32, 48]),
227+ ("double layer, graded panels", "bie", bie_solve, [4, 8, 12, 16, 20, 24]),
228+ ]
229+ E, T = sweep(solve, ns, px, py, ve)
230+ println(name)
231+ for k in eachindex(ns)
232+ @printf(" n=%3d err=%.3e %8.2f ms\n", ns[k], E[k], T[k])
233+ end
234+ push!(entries,
235+ """{"name":"$(name)","slot":"$(slot)","ns":$(json_arr(ns)),"E":$(json_arr(E)),"T":$(json_arr(T))}""")
236+ end
237+ label = "Julia $(VERSION.major).$(VERSION.minor)"
238+ open("bench-pacman-julia.json", "w") do f
239+ write(f,
240+ """{"schema":"laplace-dirichlet-bench-1","label":"$(label)","kind":"$(KIND)","date":"$(Dates.today())","methods":[$(join(entries, ","))]}""")
241+ end
242+ println("wrote bench-pacman-julia.json")
243+end
244+
245+main()
bench/matlab/bench_pacman.madded+97−0View file
@@ -0,0 +1,97 @@
1+% Offline benchmark of the Pac-Man problem in real MATLAB.
2+%
3+% Runs the tutorial's own MATLAB layer, unchanged: the chapter driver and the
4+% editable script are concatenated into main.m exactly as the page does, and
5+% executed here instead of in numbl. Only the two best methods from the
6+% in-browser benchmark are swept (the harmonic basis and the graded-panel
7+% double layer), with a deeper panel ladder than the page uses, since a native
8+% LAPACK has no trouble with the larger dense systems.
9+%
10+% matlab -batch "run('bench/matlab/bench_pacman.m')"
11+%
12+% from the repository root writes bench-pacman-matlab.json into the
13+% repository root (run() changes the working directory while it executes, so
14+% "the current directory" would be a moving target). Load that file from the
15+% benchmark chapter's "load offline results" button; because it ran on the
16+% same machine as your browser, the timings share a clock and can share an
17+% axis.
18+
19+here = fileparts(mfilename('fullpath'));
20+root = fileparts(fileparts(here));
21+outfile = fullfile(root, 'bench-pacman-matlab.json');
22+
23+% the deep basis ladders are legitimately rank deficient at the 1e-12 level;
24+% Chapter 9 discusses why, and the warning would repeat once per solve
25+warning('off', 'MATLAB:rankDeficientMatrix');
26+
27+kind = 'both';
28+% ladders: zeros skip a method (rows: SW+SOR, basis, panels, Chebyshev)
29+nsm = [0, 0, 0, 0, 0, 0;
30+ 4, 8, 16, 32, 48, 0;
31+ 4, 8, 12, 16, 20, 24;
32+ 0, 0, 0, 0, 0, 0];
33+
34+work = tempname;
35+mkdir(work);
36+libdir = fullfile(root, 'src', 'matlab');
37+libfiles = dir(fullfile(libdir, '*.m'));
38+for k = 1:numel(libfiles)
39+ copyfile(fullfile(libdir, libfiles(k).name), work);
40+end
41+driver = fileread(fullfile(root, 'src', 'chapters', 'ch11', 'driver.m'));
42+script = fileread(fullfile(root, 'src', 'chapters', 'ch11', 'pacman.m'));
43+fid = fopen(fullfile(work, 'main.m'), 'w');
44+fwrite(fid, sprintf('%s\n%s', driver, script));
45+fclose(fid);
46+
47+params = struct('kind', kind, 'mode', 'sweep', 'nsm', nsm);
48+fid = fopen(fullfile(work, 'params.json'), 'w');
49+fwrite(fid, jsonencode(params));
50+fclose(fid);
51+
52+old = cd(work);
53+cleanup = onCleanup(@() cd(old));
54+run('main.m');
55+res = jsondecode(fileread('out.json'));
56+cd(old);
57+
58+% out.json holds one row per method; keep the two that actually ran
59+slots = {'fd2', 'mps', 'bie', 'spectral'};
60+keep = [2, 3];
61+methods = cell(1, numel(keep));
62+for q = 1:numel(keep)
63+ k = keep(q);
64+ ns = pick_row(res.NS, k);
65+ ok = isfinite(ns) & ns > 0;
66+ m = struct();
67+ m.name = res.names{k};
68+ m.slot = slots{k};
69+ m.ns = ns(ok);
70+ m.E = pick_row(res.E, k);
71+ m.E = m.E(ok);
72+ m.T = pick_row(res.T, k);
73+ m.T = m.T(ok);
74+ methods{q} = m;
75+end
76+
77+out = struct();
78+out.schema = 'laplace-dirichlet-bench-1';
79+out.label = ['MATLAB ' version('-release')];
80+out.kind = kind;
81+out.date = char(datetime('now', 'Format', 'yyyy-MM-dd'));
82+out.methods = methods;
83+
84+fid = fopen(outfile, 'w');
85+fwrite(fid, jsonencode(out));
86+fclose(fid);
87+fprintf('wrote %s\n', outfile);
88+
89+function r = pick_row(M, k)
90+% jsondecode gives a matrix when the rows are rectangular and a cell of
91+% columns when they are not; accept both.
92+if iscell(M)
93+ r = M{k}(:)';
94+else
95+ r = M(k, :);
96+end
97+end
bench/python/bench_pacman.pyadded+207−0View file
@@ -0,0 +1,207 @@
1+#!/usr/bin/env python3
2+"""Offline benchmark of the Pac-Man problem in Python (NumPy).
3+
4+A port of the two best methods from the in-browser benchmark, the harmonic
5+basis and the graded-panel double layer, on the same problem, the same error
6+lattice, and the same timing protocol as the chapter driver. Run
7+
8+ python3 bench/python/bench_pacman.py
9+
10+from anywhere; it writes bench-pacman-python.json into the current directory.
11+Load that file from the benchmark chapter's "load offline results" button.
12+Because it ran on the same machine as your browser, the timings share a clock
13+and can share an axis.
14+"""
15+
16+import json
17+import platform
18+import sys
19+import time
20+from datetime import date
21+
22+import numpy as np
23+
24+KIND = "both"
25+Z0 = 1.6 + 1.3j
26+
27+
28+def exact(x, y, kind=KIND):
29+ """lp_lexact: the exact solution, singular and/or smooth part."""
30+ th = np.arctan2(y, x)
31+ th = np.where(th < -1e-13, th + 2 * np.pi, th)
32+ r = np.hypot(x, y)
33+ sing = r ** (2.0 / 3.0) * np.sin(2.0 * th / 3.0)
34+ smooth = 0.5 * np.log((x - Z0.real) ** 2 + (y - Z0.imag) ** 2)
35+ if kind == "sing":
36+ return sing
37+ if kind == "smooth":
38+ return smooth
39+ return sing + smooth
40+
41+
42+def lattice():
43+ """The driver's error lattice: spacing 1/8, margins from the boundary."""
44+ s = np.arange(-8, 9) / 8.0
45+ X, Y = np.meshgrid(s, s)
46+ R2 = X**2 + Y**2
47+ keep = (
48+ (R2 <= 0.7225 + 1e-9)
49+ & (R2 > 1e-12)
50+ & ~((X > 1e-12) & (Y < -1e-12))
51+ & ~((X > 1e-12) & (Y < 0.25 - 1e-9))
52+ & ~((Y < -1e-12) & (X > -0.25 + 1e-9))
53+ )
54+ return X[keep], Y[keep]
55+
56+
57+# ------------------------------------------- the method of particular sols --
58+
59+
60+def pac_basis(x, y, M, K):
61+ z = x + 1j * y
62+ cols = [np.ones_like(x)]
63+ p = z.copy()
64+ for _ in range(M):
65+ cols.append(p.real)
66+ cols.append(p.imag)
67+ p = p * z
68+ th = np.arctan2(y, x)
69+ th = np.where(th < -1e-13, th + 2 * np.pi, th)
70+ r = np.abs(z)
71+ k, q = 0, 0
72+ while q < K:
73+ k += 1
74+ if k % 3 != 0:
75+ q += 1
76+ cols.append(r ** (2.0 * k / 3.0) * np.sin(2.0 * k * th / 3.0))
77+ return np.column_stack(cols)
78+
79+
80+def pac_bpoints(ntot):
81+ ne = max(6, round(ntot / 6.7))
82+ na = max(12, ntot - 2 * ne)
83+ t = (np.arange(1, ne + 1) - 0.5) / ne
84+ c = t**2
85+ s = 1.5 * np.pi * (np.arange(1, na + 1) - 0.5) / na
86+ bx = np.concatenate([c, np.cos(s), np.zeros(ne)])
87+ by = np.concatenate([np.zeros(ne), np.sin(s), -c])
88+ return bx, by
89+
90+
91+def mps_solve(n, px, py):
92+ nb = 3 * n + 1
93+ bx, by = pac_bpoints(3 * nb)
94+ A = pac_basis(bx, by, n, n)
95+ # QR-based least squares, the closest match to MATLAB's backslash
96+ c, *_ = np.linalg.lstsq(A, exact(bx, by), rcond=None)
97+ return pac_basis(px, py, n, n) @ c
98+
99+
100+# --------------------------------------------- boundary integral, panels ----
101+
102+
103+def graded_breaks(T, d, base):
104+ half = T / 2.0
105+ left = np.concatenate([[0.0], half * 2.0 ** np.arange(-d, 1)])
106+ bp0 = np.concatenate([left, T - left[-2::-1]])
107+ bp = [bp0[0]]
108+ for a, b in zip(bp0[:-1], bp0[1:]):
109+ ns = max(1, int(np.ceil((b - a) / base)))
110+ for q in range(1, ns + 1):
111+ bp.append(a + (b - a) * q / ns)
112+ return np.array(bp)
113+
114+
115+def panel_nodes(bp, gx, gw):
116+ a = bp[:-1]
117+ b = bp[1:]
118+ t = ((a + b) / 2)[:, None] + ((b - a) / 2)[:, None] * gx[None, :]
119+ w = ((b - a) / 2)[:, None] * gw[None, :]
120+ return t.ravel(), w.ravel()
121+
122+
123+def pac_panels(d):
124+ gx, gw = np.polynomial.legendre.leggauss(12)
125+ be = graded_breaks(1.0, d, 0.35)
126+ ba = graded_breaks(1.5 * np.pi, d, 0.35)
127+ t1, w1 = panel_nodes(be, gx, gw)
128+ t2, w2 = panel_nodes(ba, gx, gw)
129+ t3, w3 = panel_nodes(be, gx, gw)
130+ z1 = np.zeros_like(t1)
131+ z3 = np.zeros_like(t3)
132+ qx = np.concatenate([t1, np.cos(t2), z3])
133+ qy = np.concatenate([z1, np.sin(t2), t3 - 1.0])
134+ qnx = np.concatenate([z1, np.cos(t2), np.ones_like(t3)])
135+ qny = np.concatenate([-np.ones_like(t1), np.sin(t2), z3])
136+ qw = np.concatenate([w1, w2, w3])
137+ qkap = np.concatenate([z1, np.ones_like(t2), z3])
138+ return qx, qy, qnx, qny, qw, qkap
139+
140+
141+def bie_solve(d, px, py):
142+ qx, qy, qnx, qny, qw, qkap = pac_panels(d)
143+ n = qx.size
144+ DX = qx[:, None] - qx[None, :]
145+ DY = qy[:, None] - qy[None, :]
146+ R2 = DX**2 + DY**2
147+ np.fill_diagonal(R2, 1.0)
148+ K = (qnx[None, :] * DX + qny[None, :] * DY) / (2 * np.pi * R2)
149+ np.fill_diagonal(K, -qkap / (4 * np.pi))
150+ A = -0.5 * np.eye(n) + K * qw[None, :]
151+ sigma = np.linalg.solve(A, exact(qx, qy))
152+ EX = px[:, None] - qx[None, :]
153+ EY = py[:, None] - qy[None, :]
154+ KE = (qnx[None, :] * EX + qny[None, :] * EY) / (2 * np.pi * (EX**2 + EY**2))
155+ return (KE * qw[None, :]) @ sigma
156+
157+
158+# ------------------------------------------------------------- the sweep ----
159+
160+
161+def sweep(solve, ns, px, py, ve):
162+ E, T = [], []
163+ solve(ns[0], px[:1], py[:1]) # warm-up
164+ for n in ns:
165+ t0 = time.perf_counter()
166+ v = solve(n, px, py)
167+ el = time.perf_counter() - t0
168+ rep = 0
169+ while el < 0.005 and rep < 4:
170+ rep += 1
171+ t0 = time.perf_counter()
172+ v = solve(n, px, py)
173+ el = min(el, time.perf_counter() - t0)
174+ E.append(float(np.max(np.abs(v - ve))))
175+ T.append(el * 1000.0)
176+ return E, T
177+
178+
179+def main():
180+ px, py = lattice()
181+ ve = exact(px, py)
182+ methods = []
183+ for name, slot, solve, ns in [
184+ ("harmonic basis, least squares", "mps", mps_solve, [4, 8, 16, 32, 48]),
185+ ("double layer, graded panels", "bie", bie_solve, [4, 8, 12, 16, 20, 24]),
186+ ]:
187+ E, T = sweep(solve, ns, px, py, ve)
188+ methods.append({"name": name, "slot": slot, "ns": ns, "E": E, "T": T})
189+ print(name)
190+ for n, e, t in zip(ns, E, T):
191+ print(f" n={n:3d} err={e:.3e} {t:8.2f} ms")
192+
193+ out = {
194+ "schema": "laplace-dirichlet-bench-1",
195+ "label": f"Python {sys.version_info.major}.{sys.version_info.minor} + NumPy {np.__version__}",
196+ "kind": KIND,
197+ "date": date.today().isoformat(),
198+ "machine": platform.machine(),
199+ "methods": methods,
200+ }
201+ with open("bench-pacman-python.json", "w") as f:
202+ json.dump(out, f)
203+ print("wrote bench-pacman-python.json")
204+
205+
206+if __name__ == "__main__":
207+ main()
scripts/matlab-test.mjsmodified+53−0View file
@@ -231,5 +231,58 @@ console.log('\nChapter 10: boundary integral equations')
231231 )
232232 }
233233
234+console.log('\nChapter 11: the Pac-Man benchmark')
235+{
236+ // zeros in a ladder skip that resolution (or the whole method)
237+ const { out } = runLab(src('chapters/ch11/driver.m'), src('chapters/ch11/pacman.m'), {
238+ kind: 'both',
239+ mode: 'sweep',
240+ nsm: [
241+ [16, 32, 64, 0],
242+ [4, 8, 16, 32],
243+ [4, 8, 0, 0],
244+ [8, 16, 32, 0],
245+ ],
246+ })
247+ const qfd = last(out.orders[0])
248+ check('Shortley-Weller is corner-limited near order 4/3', near(qfd, 1.37, 0.12), `q = ${qfd.toFixed(3)}`)
249+ check('the harmonic basis reaches 1e-13 with 97 unknowns', out.E[1][3] < 1e-13, out.E[1][3].toExponential(2))
250+ check('and converges geometrically', out.E[1][1] / out.E[1][2] > 100)
251+ const ratio = out.E[2][0] / out.E[2][1]
252+ check(
253+ 'graded panels gain about 2^(8/3) per four levels of depth',
254+ ratio > 4 && ratio < 9,
255+ `ratio = ${ratio.toFixed(2)}`,
256+ )
257+ check('and reach 3e-5 by depth 8', out.E[2][1] < 3e-5, out.E[2][1].toExponential(2))
258+ check('polar Chebyshev is algebraic on the singular data', out.E[3][2] > 1e-5 && out.E[3][2] < 3e-4, out.E[3][2].toExponential(2))
259+
260+ const { out: sm } = runLab(src('chapters/ch11/driver.m'), src('chapters/ch11/pacman.m'), {
261+ kind: 'smooth',
262+ mode: 'sweep',
263+ nsm: [
264+ [16, 32, 64, 0],
265+ [0, 0, 0, 0],
266+ [0, 0, 0, 0],
267+ [6, 16, 40, 0],
268+ ],
269+ })
270+ const qsw = last(sm.orders[0])
271+ check('with smooth data Shortley-Weller is second order', near(qsw, 2, 0.15), `q = ${qsw.toFixed(3)}`)
272+ check('and polar Chebyshev is geometric again, 1e-9 by N = 40', sm.E[3][2] < 1e-9, sm.E[3][2].toExponential(2))
273+
274+ const { out: sg } = runLab(src('chapters/ch11/driver.m'), src('chapters/ch11/pacman.m'), {
275+ kind: 'sing',
276+ mode: 'sweep',
277+ nsm: [
278+ [0, 0, 0, 0],
279+ [4, 0, 0, 0],
280+ [0, 0, 0, 0],
281+ [0, 0, 0, 0],
282+ ],
283+ })
284+ check('the singular function is in the basis exactly', sg.E[1][0] < 1e-13, sg.E[1][0].toExponential(2))
285+}
286+
234287 console.log(failures === 0 ? '\nall checks passed' : `\n${failures} check(s) failed`)
235288 process.exit(failures === 0 ? 0 : 1)
src/chapters/ch10/text.mdmodified+3−1View file
@@ -56,9 +56,11 @@ The first row is the point of the whole tutorial in one line: full double precis
5656
5757 ### Where this leaves the tutorial
5858
59-The methods in these ten chapters differ by fourteen orders of magnitude in the accuracy they reach for a comparable amount of work, on problems whose exact solutions were known throughout. The pattern behind that spread is worth stating plainly.
59+The methods in these chapters differ by fourteen orders of magnitude in the accuracy they reach for a comparable amount of work, on problems whose exact solutions were known throughout. The pattern behind that spread is worth stating plainly.
6060
6161 - **Accuracy is bounded by the smoothness of the solution, and only then by the method.** Chapters 4 and 9 are the same lesson twice: a spectral method on a solution with a $|x-\tfrac12|^{5/2}$ kink is no better than a three-point difference, and second-order differences on a re-entrant corner are no better than order $4/3$.
6262 - **Where the solution is smooth, high order is nearly free and enormously profitable.** The compact nine-point stencil costs the same per grid point as the five-point stencil and is ten orders of magnitude more accurate at $N = 32$.
6363 - **Structure, when it exists, is worth more than any amount of brute force.** The sine transform, the Chebyshev matrix decomposition and the second-kind integral equation are all the same move: identify what makes the problem special and let the solver use it.
6464 - **Every one of these statements has conditions, and the conditions are the useful part.** A separable operator on a rectangle, an analytic solution, a smooth boundary, a target not too close to it. Knowing which condition a method needs is what tells you when to reach for it, and the labs above are set up so that each condition can be switched off and the consequences watched.
65+
66+The next chapter is the test of those statements: one domain that grants nobody their favourite condition, all four families of methods on it, and one plot of accuracy against time.
src/chapters/ch11/chapter.tsadded+132−0View file
@@ -0,0 +1,132 @@
1+import { createElement } from 'react'
2+import type { Chapter, Control } from '../../ui/types.ts'
3+import driver from './driver.m?raw'
4+import script from './pacman.m?raw'
5+import text from './text.md?raw'
6+import { Fields2D, type Single2DOut } from '../twodim/panels.tsx'
7+import { fmtErr } from '../onedim/panels.tsx'
8+import { BenchSweep, type BenchOut } from './panels.tsx'
9+import type { Num } from '../../engine/types.ts'
10+
11+interface FieldOut extends Single2DOut {
12+ cx: number[]
13+ cy: number[]
14+}
15+
16+const KIND: Control = {
17+ kind: 'select',
18+ key: 'kind',
19+ label: 'ground truth',
20+ options: [
21+ { value: 'both', label: 'singular + smooth (generic)' },
22+ { value: 'sing', label: 'singular only: r^(2/3) sin(2θ/3)' },
23+ { value: 'smooth', label: 'smooth only: log|z − z₀|' },
24+ ],
25+}
26+
27+/** Each method's resolution ladder behind the shared coarse-to-finest slider. */
28+const LEVELS: number[][] = [
29+ [32, 64, 96, 128], // grid intervals per unit
30+ [8, 16, 32, 48], // basis degree
31+ [4, 6, 8, 10], // refinement depth
32+ [10, 20, 30, 40], // Chebyshev points per direction
33+]
34+
35+const NSM = [
36+ [8, 16, 32, 64, 128],
37+ [4, 8, 16, 32, 48],
38+ [2, 4, 6, 8, 10],
39+ [6, 10, 16, 26, 40],
40+]
41+
42+function bestOf(d: BenchOut): { v: number; name: string } {
43+ let v = Infinity
44+ let name = ''
45+ d.E.forEach((row, k) => {
46+ for (const e of row) {
47+ if (e != null && isFinite(e) && e > 0 && e < v) {
48+ v = e
49+ name = d.names[k]
50+ }
51+ }
52+ })
53+ return { v, name }
54+}
55+
56+export const ch11: Chapter = {
57+ id: 'benchmarks',
58+ nav: 'One problem, every method',
59+ part: 'Beyond the square',
60+ title: 'Benchmarks: one problem, every method',
61+ blurb:
62+ 'A curved boundary with three corners, the four families of methods side by side, and timings you can extend beyond the browser.',
63+ text,
64+ labs: [
65+ {
66+ id: 'field',
67+ title: 'The solution, and where each method hurts',
68+ caption:
69+ 'The error picture is the method’s signature: the difference scheme concentrates its error at the corner, the double layer at the boundary (that ring is the close-evaluation problem of Chapter 10), and the collocation methods spread theirs.',
70+ driver,
71+ script,
72+ manual: true,
73+ controls: [
74+ KIND,
75+ {
76+ kind: 'select',
77+ key: 'mid',
78+ label: 'method',
79+ options: [
80+ { value: 1, label: 'Shortley-Weller + SOR' },
81+ { value: 2, label: 'harmonic basis' },
82+ { value: 3, label: 'double layer, graded panels' },
83+ { value: 4, label: 'Chebyshev, polar rectangle' },
84+ ],
85+ },
86+ {
87+ kind: 'slider',
88+ key: 'level',
89+ label: 'resolution',
90+ values: [1, 2, 3, 4],
91+ min: 1,
92+ max: 4,
93+ format: (v) => ['coarse', 'medium', 'fine', 'finest'][v - 1],
94+ },
95+ ],
96+ defaults: { kind: 'both', mode: 'field', mid: 1, level: 2 },
97+ derive: (p) => ({ ...p, n: LEVELS[Number(p.mid) - 1][Number(p.level) - 1] }),
98+ render: (d: FieldOut) =>
99+ createElement(Fields2D, {
100+ data: d,
101+ xDomain: [-1.0625, 1.0625],
102+ yDomain: [-1.0625, 1.0625],
103+ curves: [{ x: d.cx, y: d.cy }],
104+ }),
105+ summary: (d: FieldOut) => [
106+ { label: 'unknowns', value: String(d.unknowns) },
107+ { label: 'max error', value: fmtErr(d.maxerr) },
108+ ],
109+ },
110+ {
111+ id: 'bench',
112+ title: 'The benchmark',
113+ caption:
114+ 'The first plot is only meaningful curve by curve, because each method refines its own parameter; the second, against wall clock, is the comparison. Results from the offline scripts land in the second plot as dashed curves.',
115+ driver,
116+ script,
117+ manual: true,
118+ editorOpen: true,
119+ controls: [KIND],
120+ defaults: { kind: 'both', mode: 'sweep', nsm: NSM },
121+ render: (d: BenchOut, p) => createElement(BenchSweep, { data: d, kind: String(p.kind) }),
122+ summary: (d: BenchOut) => {
123+ const b = bestOf(d)
124+ const t = d.T.flat().reduce((s: number, v: Num) => s + (v ?? 0), 0)
125+ return [
126+ { label: 'best error', value: `${fmtErr(b.v)} (${b.name})` },
127+ { label: 'time in solves', value: t < 1000 ? `${t.toFixed(0)} ms` : `${(t / 1000).toFixed(1)} s` },
128+ ]
129+ },
130+ },
131+ ],
132+}
src/chapters/ch11/driver.madded+143−0View file
@@ -0,0 +1,143 @@
1+% Driver for the benchmark chapter: one hard problem, every method.
2+%
3+% The editable script supplies
4+%
5+% names = method_names()
6+% v = solve_pac(mid, n, g, px, py)
7+% d = dof_pac(mid, n)
8+%
9+% and must return the computed solution at the points (px, py), whatever grid,
10+% basis, or boundary discretisation it used internally. n is the method's own
11+% resolution parameter: intervals per unit length for the difference scheme,
12+% basis terms for the expansion, refinement depth for the panels, points per
13+% direction for the collocation. dof_pac reports the number of unknowns behind
14+% that parameter, which the page shows so that the parameters stay comparable.
15+%
16+% g(x, y) is the Dirichlet data: the exact solution restricted to the
17+% boundary. The domain is the Pac-Man
18+%
19+% 0 < r < 1, 0 < theta < 3*pi/2,
20+%
21+% and the exact solutions are those of the L-shaped chapter (lp_lexact), whose
22+% corner sits at the same origin with the same interior angle.
23+p = lp_params();
24+kind = p.kind;
25+names = method_names();
26+nm = numel(names);
27+g = @(xx, yy) lp_lexact(kind, xx, yy);
28+
29+out = struct();
30+out.names = names;
31+
32+% the boundary, drawn over the field pictures: the edge theta = 0, the arc,
33+% and the edge theta = 3*pi/2, walked as one closed polyline
34+ta = (0:120)' * (3*pi/2) / 120;
35+e1 = linspace(0, 1, 15)';
36+e2 = linspace(-1, 0, 15)';
37+out.cx = [e1; cos(ta); zeros(15, 1)]';
38+out.cy = [zeros(15, 1); sin(ta); e2]';
39+
40+if strcmp(p.mode, 'field')
41+ % a picture on a fixed grid of spacing 1/32, masked to the open domain;
42+ % pixels on the two straight edges are left out so that the boundary
43+ % methods are never evaluated exactly on the curve they integrate over
44+ m = 34;
45+ xs = (-m:m)' / 32;
46+ [X, Y] = meshgrid(xs, xs);
47+ R2 = X.^2 + Y.^2;
48+ onedge = (abs(Y) < 1e-12 & X > -1e-12) | (abs(X) < 1e-12 & Y < 1e-12);
49+ inside = R2 < 1 - 1e-9 & ~(X > 1e-12 & Y < -1e-12) & ~onedge;
50+ px = X(inside);
51+ py = Y(inside);
52+ t0 = tic;
53+ v = solve_pac(p.mid, p.n, g, px, py);
54+ out.ms = toc(t0) * 1000;
55+ ve = g(px, py);
56+ out.maxerr = max(abs(v - ve));
57+ out.unknowns = dof_pac(p.mid, p.n);
58+
59+ U = NaN(2*m+1, 2*m+1);
60+ E = NaN(2*m+1, 2*m+1);
61+ U(inside) = v;
62+ E(inside) = abs(v - ve);
63+ out.x = xs';
64+ out.y = xs';
65+ out.U = lp_rows(U);
66+ out.Err = lp_rows(E);
67+else
68+ % errors at the nodes of a fixed lattice of spacing 1/8, kept at least
69+ % 0.15 from the arc and 0.25 from the straight edges, and reaching within
70+ % 0.125 of the re-entrant corner (but excluding the corner itself, which
71+ % is a boundary point)
72+ [X, Y] = meshgrid((-8:8)' / 8);
73+ R2 = X.^2 + Y.^2;
74+ keep = R2 <= 0.7225 + 1e-9 & R2 > 1e-12 & ~(X > 1e-12 & Y < -1e-12) ...
75+ & ~(X > 1e-12 & Y < 0.25 - 1e-9) & ~(Y < -1e-12 & X > -0.25 + 1e-9);
76+ px = X(keep);
77+ py = Y(keep);
78+ ve = g(px, py);
79+
80+ % Each method gets its own list of resolutions. numbl's jsondecode turns
81+ % a JSON array of arrays into a cell of column vectors where MATLAB gives
82+ % a matrix; the offline harness runs this same file in real MATLAB, so
83+ % both shapes are accepted.
84+ NSc = p.nsm;
85+ if iscell(NSc)
86+ nn = numel(NSc{1});
87+ NS = zeros(nm, nn);
88+ for k = 1:nm
89+ row = NSc{k};
90+ NS(k, :) = row(:)';
91+ end
92+ else
93+ NS = NSc;
94+ nn = size(NS, 2);
95+ end
96+ E = NaN(nm, nn);
97+ T = NaN(nm, nn);
98+ Q = NaN(nm, nn);
99+ DOF = NaN(nm, nn);
100+ for k = 1:nm
101+ % a ladder of zeros skips the method entirely, which is how the
102+ % offline harness runs a subset of the methods through this driver
103+ nw = 0;
104+ for j = 1:nn
105+ if isfinite(NS(k, j)) && NS(k, j) > 0
106+ nw = NS(k, j);
107+ break;
108+ end
109+ end
110+ if nw == 0
111+ continue;
112+ end
113+ vw = solve_pac(k, nw, g, px(1), py(1)); % warm-up
114+ for j = 1:nn
115+ n = NS(k, j);
116+ if ~isfinite(n) || n <= 0
117+ continue;
118+ end
119+ t0 = tic;
120+ v = solve_pac(k, n, g, px, py);
121+ el = toc(t0);
122+ rep = 0;
123+ while el < 0.005 && rep < 4
124+ rep = rep + 1;
125+ t0 = tic;
126+ v = solve_pac(k, n, g, px, py);
127+ el = min(el, toc(t0));
128+ end
129+ T(k, j) = el * 1000;
130+ E(k, j) = max(abs(v - ve));
131+ DOF(k, j) = dof_pac(k, n);
132+ end
133+ Q(k, :) = lp_orders(E(k, :), NS(k, :));
134+ end
135+ out.ns = NS(1, :);
136+ out.NS = lp_rows(NS);
137+ out.DOF = lp_rows(DOF);
138+ out.E = lp_rows(E);
139+ out.T = lp_rows(T);
140+ out.orders = lp_rows(Q);
141+end
142+
143+lp_write(out);
src/chapters/ch11/pacman.madded+407−0View file
@@ -0,0 +1,407 @@
1+% Four methods for the Dirichlet problem on the Pac-Man domain: the unit disk
2+% with the quadrant x > 0, y < 0 removed. The boundary is a curved arc and two
3+% straight edges; the edges meet the arc in two right-angle corners, and each
4+% other in a re-entrant corner of interior angle 3*pi/2 at the origin. Every
5+% difficulty this tutorial has met is present at once: a boundary that cuts
6+% through grid cells (Chapter 10's curve), and a corner singularity r^(2/3)
7+% (Chapter 9's corner).
8+%
9+% 1. SHORTLEY-WELLER + SOR. The five-point scheme, with the legs of the
10+% stencil shortened wherever a neighbour falls outside the arc, so the
11+% boundary condition is imposed on the arc itself and not on a staircase
12+% approximation of it. The shortened stencil is no longer symmetric, so
13+% conjugate gradients is not available; the solver is SOR with the
14+% overrelaxation parameter near its optimum.
15+%
16+% 2. THE HARMONIC BASIS. The method of particular solutions from Chapter 9,
17+% unchanged in spirit: harmonic polynomials Re z^m, Im z^m for the smooth
18+% part, plus the corner's own functions r^(2k/3) sin(2k theta/3), fitted
19+% to the boundary data by least squares.
20+%
21+% 3. DOUBLE LAYER, GRADED PANELS. The Nystrom method of Chapter 10 cannot
22+% use the trapezoid rule here, because the boundary is not smooth: the
23+% density has its own singularity at each corner. The quadrature is
24+% composite Gauss-Legendre on panels refined dyadically towards the three
25+% corners, which resolves the density algebraically panel by panel and
26+% geometrically in the refinement depth.
27+%
28+% 4. CHEBYSHEV IN POLAR COORDINATES. In polar coordinates the Pac-Man is a
29+% rectangle, (r, theta) in (0,1) x (0, 3*pi/2), and Laplace's equation
30+% becomes r^2 u_rr + r u_r + u_tt = 0, which is separable: the matrix
31+% decomposition of Chapter 7 applies verbatim. The whole side r = 0 of
32+% the rectangle maps to the corner, where u is known, so the coordinate
33+% degeneracy costs nothing to impose; what it costs is accuracy, because
34+% the r^(2/3) singularity now sits on the edge of the expansion domain.
35+
36+function names = method_names()
37+names = {'Shortley-Weller + SOR', 'harmonic basis, least squares', ...
38+ 'double layer, graded panels', 'Chebyshev, polar rectangle'};
39+end
40+
41+function v = solve_pac(mid, n, g, px, py)
42+if mid == 1
43+ v = fd_pac(n, g, px, py);
44+elseif mid == 2
45+ v = mps_pac(n, g, px, py);
46+elseif mid == 3
47+ v = bie_pac(n, g, px, py);
48+else
49+ v = spec_pac(n, g, px, py);
50+end
51+end
52+
53+function d = dof_pac(mid, n)
54+if mid == 1
55+ d = fd_count(n);
56+elseif mid == 2
57+ d = 3*n + 1;
58+elseif mid == 3
59+ [qx, qy, qnx, qny, qw, qkap] = pac_panels(n);
60+ d = numel(qx);
61+else
62+ d = (n - 1)^2;
63+end
64+end
65+
66+% --------------------------------------------------- finite differences -----
67+
68+function v = fd_pac(N, g, px, py)
69+% N intervals per unit length on (-1:1/N:1)^2. Nodes on the two straight
70+% edges are grid points and carry the data directly; nodes adjacent to the
71+% arc get Shortley-Weller legs, shortened to where the arc actually is.
72+h = 1 / N;
73+m = 2*N + 1;
74+xs = (-N:N)' * h;
75+[X, Y] = meshgrid(xs, xs);
76+R2 = X.^2 + Y.^2;
77+
78+% start from the data everywhere: boundary nodes keep it, and for the
79+% unknowns it is a continuous extension, which is a decent initial guess
80+U = g(X, Y);
81+inquad = X > 1e-12 & Y < -1e-12;
82+U(inquad) = 0;
83+onedge = (abs(Y) < 1e-12 & X > -1e-12) | (abs(X) < 1e-12 & Y < 1e-12);
84+unknown = R2 < 1 - 1e-12 & ~inquad & ~onedge;
85+mask = double(unknown);
86+
87+% Shortley-Weller coefficients. Every leg is h except where the neighbour
88+% lies beyond the arc; there the leg is shortened to the crossing and the
89+% boundary value at the crossing moves to the right-hand side.
90+CE = zeros(m, m);
91+CW = zeros(m, m);
92+CN = zeros(m, m);
93+CS = zeros(m, m);
94+CP = ones(m, m);
95+B = zeros(m, m);
96+for j = 2:m-1
97+ for i = 2:m-1
98+ if mask(i, j) == 0
99+ continue;
100+ end
101+ x = xs(j);
102+ y = xs(i);
103+ hE = h; hW = h; hN = h; hS = h;
104+ gE = 0; gW = 0; gN = 0; gS = 0;
105+ if (x + h)^2 + y^2 >= 1 - 1e-12
106+ hE = max(sqrt(1 - y^2) - x, 1e-8 * h);
107+ gE = 1;
108+ end
109+ if (x - h)^2 + y^2 >= 1 - 1e-12
110+ hW = max(x + sqrt(1 - y^2), 1e-8 * h);
111+ gW = 1;
112+ end
113+ if x^2 + (y + h)^2 >= 1 - 1e-12
114+ hN = max(sqrt(1 - x^2) - y, 1e-8 * h);
115+ gN = 1;
116+ end
117+ if x^2 + (y - h)^2 >= 1 - 1e-12
118+ hS = max(y + sqrt(1 - x^2), 1e-8 * h);
119+ gS = 1;
120+ end
121+ cE = 2 / (hE * (hE + hW));
122+ cW = 2 / (hW * (hE + hW));
123+ cN = 2 / (hN * (hN + hS));
124+ cS = 2 / (hS * (hN + hS));
125+ CP(i, j) = cE + cW + cN + cS;
126+ % a shortened leg reaches a point of the arc, not a grid node: its
127+ % (known) value multiplies the coefficient and moves to the rhs
128+ if gE == 1
129+ B(i, j) = B(i, j) + cE * g(x + hE, y);
130+ else
131+ CE(i, j) = cE;
132+ end
133+ if gW == 1
134+ B(i, j) = B(i, j) + cW * g(x - hW, y);
135+ else
136+ CW(i, j) = cW;
137+ end
138+ if gN == 1
139+ B(i, j) = B(i, j) + cN * g(x, y + hN);
140+ else
141+ CN(i, j) = cN;
142+ end
143+ if gS == 1
144+ B(i, j) = B(i, j) + cS * g(x, y - hS);
145+ else
146+ CS(i, j) = cS;
147+ end
148+ end
149+end
150+
151+% SOR, swept until the largest Gauss-Seidel correction is far below the
152+% discretisation error. The shortened stencils make the matrix
153+% nonsymmetric, which is what rules out the CG solver of Chapter 9.
154+om = 2 / (1 + sin(pi * h / 2));
155+maxit = 60 * N;
156+for it = 1:maxit
157+ du = 0;
158+ for j = 2:m-1
159+ for i = 2:m-1
160+ if mask(i, j) > 0
161+ unew = (CE(i,j) * U(i,j+1) + CW(i,j) * U(i,j-1) ...
162+ + CN(i,j) * U(i+1,j) + CS(i,j) * U(i-1,j) + B(i,j)) / CP(i,j);
163+ d = unew - U(i, j);
164+ U(i, j) = U(i, j) + om * d;
165+ if abs(d) > du
166+ du = abs(d);
167+ end
168+ end
169+ end
170+ end
171+ if du < 1e-9
172+ break;
173+ end
174+end
175+
176+% the requested points are grid points, so sampling is a lookup
177+v = zeros(numel(px), 1);
178+for k = 1:numel(px)
179+ i = round(py(k) / h) + N + 1;
180+ j = round(px(k) / h) + N + 1;
181+ v(k) = U(i, j);
182+end
183+end
184+
185+function d = fd_count(N)
186+h = 1 / N;
187+xs = (-N:N)' * h;
188+[X, Y] = meshgrid(xs, xs);
189+R2 = X.^2 + Y.^2;
190+inquad = X > 1e-12 & Y < -1e-12;
191+onedge = (abs(Y) < 1e-12 & X > -1e-12) | (abs(X) < 1e-12 & Y < 1e-12);
192+d = sum(sum(double(R2 < 1 - 1e-12 & ~inquad & ~onedge)));
193+end
194+
195+% ------------------------------------------- the method of particular sols --
196+
197+function v = mps_pac(n, g, px, py)
198+% Harmonic polynomials up to degree n plus n corner functions, fitted by
199+% least squares at about three boundary points per unknown.
200+nb = 3*n + 1;
201+[bx, by] = pac_bpoints(3 * nb);
202+A = pac_basis(bx, by, n, n);
203+c = A \ g(bx, by);
204+v = pac_basis(px, py, n, n) * c;
205+end
206+
207+function A = pac_basis(x, y, M, K)
208+% Re z^m, Im z^m for m = 0..M, then r^(2k/3) sin(2k theta/3) for K values of
209+% k that are not multiples of 3 (those already are harmonic polynomials, and
210+% listing them twice would make the matrix exactly rank deficient).
211+np = numel(x);
212+A = zeros(np, 2*M + 1 + K);
213+z = x(:) + 1i * y(:);
214+col = 1;
215+p = ones(np, 1);
216+for mdeg = 0:M
217+ A(:, col) = real(p);
218+ col = col + 1;
219+ if mdeg > 0
220+ A(:, col) = imag(p);
221+ col = col + 1;
222+ end
223+ p = p .* z;
224+end
225+th = atan2(y(:), x(:));
226+th(th < -1e-13) = th(th < -1e-13) + 2*pi;
227+r = abs(z);
228+k = 0;
229+q = 0;
230+while q < K
231+ k = k + 1;
232+ if mod(k, 3) ~= 0
233+ q = q + 1;
234+ A(:, col) = r.^(2*k/3) .* sin(2*k*th/3);
235+ col = col + 1;
236+ end
237+end
238+end
239+
240+function [bx, by] = pac_bpoints(ntot)
241+% Boundary points: on each straight edge, clustered quadratically towards
242+% the re-entrant corner, where the data varies fastest; on the arc, spaced
243+% uniformly in angle. Points per piece in proportion to its length.
244+ne = max(6, round(ntot / 6.7));
245+na = max(12, ntot - 2*ne);
246+t = ((1:ne)' - 0.5) / ne;
247+c = t.^2;
248+s = 3*pi/2 * ((1:na)' - 0.5) / na;
249+bx = [c; cos(s); zeros(ne, 1)];
250+by = [zeros(ne, 1); sin(s); -c];
251+end
252+
253+% --------------------------------------------- boundary integral, panels ----
254+
255+function v = bie_pac(d, g, px, py)
256+% The double-layer representation and jump relation of Chapter 10, with the
257+% trapezoid rule replaced by 12-point Gauss-Legendre panels, refined
258+% dyadically to depth d towards each of the three corners. On the straight
259+% edges the kernel between two points of the same edge is identically zero,
260+% and on the arc it is the constant -1/(4*pi); everything nontrivial in this
261+% matrix happens where two pieces meet at a corner.
262+[qx, qy, qnx, qny, qw, qkap] = pac_panels(d);
263+n = numel(qx);
264+
265+DX = qx * ones(1, n) - ones(n, 1) * qx';
266+DY = qy * ones(1, n) - ones(n, 1) * qy';
267+R2 = DX.^2 + DY.^2;
268+R2(1:n+1:end) = 1; % keep 0/0 off the diagonal
269+K = ((ones(n, 1) * qnx') .* DX + (ones(n, 1) * qny') .* DY) ./ (2*pi*R2);
270+K(1:n+1:end) = -qkap / (4*pi); % the smooth-limit diagonal
271+A = -0.5 * eye(n) + K .* (ones(n, 1) * qw');
272+
273+sigma = A \ g(qx, qy);
274+
275+np = numel(px);
276+EX = px(:) * ones(1, n) - ones(np, 1) * qx';
277+EY = py(:) * ones(1, n) - ones(np, 1) * qy';
278+KE = ((ones(np, 1) * qnx') .* EX + (ones(np, 1) * qny') .* EY) ...
279+ ./ (2*pi*(EX.^2 + EY.^2));
280+v = (KE .* (ones(np, 1) * qw')) * sigma;
281+end
282+
283+function [qx, qy, qnx, qny, qw, qkap] = pac_panels(d)
284+% Panel endpoints on each piece: dyadic towards both ends, with no panel
285+% longer than 0.35, then 12-point Gauss-Legendre on every panel. All three
286+% pieces are parametrised by arclength, so the speed is 1 throughout.
287+[gx, gw] = lp_gauss(12);
288+be = graded_breaks(1, d, 0.35);
289+ba = graded_breaks(3*pi/2, d, 0.35);
290+
291+[t1, w1] = panel_nodes(be, gx, gw);
292+[t2, w2] = panel_nodes(ba, gx, gw);
293+[t3, w3] = panel_nodes(be, gx, gw);
294+
295+% edge theta = 0 from the corner to (1,0); the arc; the edge theta = 3*pi/2
296+% from (0,-1) back towards the corner. Outward normals for a domain kept on
297+% the left of this traversal.
298+qx = [t1; cos(t2); zeros(numel(t3), 1)];
299+qy = [zeros(numel(t1), 1); sin(t2); t3 - 1];
300+qnx = [zeros(numel(t1), 1); cos(t2); ones(numel(t3), 1)];
301+qny = [-ones(numel(t1), 1); sin(t2); zeros(numel(t3), 1)];
302+qw = [w1; w2; w3];
303+qkap = [zeros(numel(t1), 1); ones(numel(t2), 1); zeros(numel(t3), 1)];
304+end
305+
306+function bp = graded_breaks(T, d, base)
307+% Breakpoints on [0, T]: halve towards each end d times, then split anything
308+% still longer than base into equal pieces.
309+half = T / 2;
310+left = [0, half * 2.^(-d:0)];
311+bp0 = [left, T - left(end-1:-1:1)];
312+bp = bp0(1);
313+for k = 2:numel(bp0)
314+ len = bp0(k) - bp0(k-1);
315+ ns = max(1, ceil(len / base));
316+ for q = 1:ns
317+ bp = [bp, bp0(k-1) + len * q / ns];
318+ end
319+end
320+end
321+
322+function [t, w] = panel_nodes(bp, gx, gw)
323+np = numel(bp) - 1;
324+p = numel(gx);
325+t = zeros(np * p, 1);
326+w = zeros(np * p, 1);
327+for k = 1:np
328+ a = bp(k);
329+ b = bp(k+1);
330+ t((k-1)*p+1 : k*p) = (a + b)/2 + (b - a)/2 * gx;
331+ w((k-1)*p+1 : k*p) = (b - a)/2 * gw;
332+end
333+end
334+
335+% --------------------------------------- Chebyshev on the polar rectangle ---
336+
337+function v = spec_pac(N, g, px, py)
338+% Collocation on the tensor Chebyshev grid of the rectangle
339+% (r, theta) in [0,1] x [0, 3*pi/2], where Laplace's equation reads
340+% r^2 u_rr + r u_r + u_tt = 0. The side r = 1 carries the arc data, the
341+% sides theta = 0 and theta = 3*pi/2 the edge data, and the side r = 0 is
342+% the corner itself, where u equals the corner value of the data.
343+[D, s] = lp_cheb(N);
344+r = (s + 1) / 2; % from 1 down to 0
345+th = (s + 1) * (3*pi/4); % from 3*pi/2 down to 0
346+Dr = 2 * D;
347+Dt = (4 / (3*pi)) * D;
348+Ar = diag(r.^2) * (Dr * Dr) + diag(r) * Dr;
349+Bt = Dt * Dt;
350+
351+G = zeros(N+1, N+1);
352+G(1, :) = g(cos(th'), sin(th')); % r = 1: the arc
353+G(N+1, :) = g(0, 0) * ones(1, N+1); % r = 0: the corner
354+G(:, 1) = g(zeros(N+1, 1), -r); % theta = 3*pi/2
355+G(:, N+1) = g(r, zeros(N+1, 1)); % theta = 0
356+
357+ii = 2:N;
358+nn = N - 1;
359+A = Ar(ii, ii);
360+Bi = Bt(ii, ii);
361+F = -(Ar(ii, [1, N+1]) * G([1, N+1], ii)) - (G(ii, [1, N+1]) * Bt(ii, [1, N+1])');
362+
363+% the matrix decomposition of Chapter 7, with two different one-dimensional
364+% operators instead of one
365+[V, LA] = eig(A);
366+[W, LB] = eig(Bi);
367+la = diag(LA);
368+lb = diag(LB);
369+Ft = V \ (F / W.');
370+Ut = Ft ./ (la * ones(1, nn) + ones(nn, 1) * lb.');
371+U = G;
372+U(ii, ii) = real(V * Ut * W.');
373+
374+% evaluate off the grid by barycentric interpolation in each direction
375+rr = sqrt(px.^2 + py.^2);
376+tt = atan2(py, px);
377+tt(tt < -1e-13) = tt(tt < -1e-13) + 2*pi;
378+v = bary2(U, s, 2*rr - 1, tt / (3*pi/4) - 1);
379+end
380+
381+function v = bary2(U, s, sr, st)
382+% Tensor barycentric interpolation on the Chebyshev grid s (both
383+% directions), rows of U indexed by the first coordinate.
384+n1 = numel(s);
385+wb = [0.5; ones(n1-2, 1); 0.5] .* (-1).^((0:n1-1)');
386+np = numel(sr);
387+v = zeros(np, 1);
388+for k = 1:np
389+ dr = sr(k) - s;
390+ dt = st(k) - s;
391+ hr = abs(dr) < 1e-14;
392+ ht = abs(dt) < 1e-14;
393+ if any(hr)
394+ cr = double(hr);
395+ else
396+ cr = wb ./ dr;
397+ cr = cr / sum(cr);
398+ end
399+ if any(ht)
400+ ct = double(ht);
401+ else
402+ ct = wb ./ dt;
403+ ct = ct / sum(ct);
404+ end
405+ v(k) = cr' * U * ct;
406+end
407+end
src/chapters/ch11/panels.tsxadded+142−0View file
@@ -0,0 +1,142 @@
1+// Figures for the benchmark chapter: the usual convergence and
2+// work-precision panels for the in-browser sweep, plus a slot for results
3+// recorded outside the browser and loaded from a JSON file. The point of the
4+// upload is that the reader runs the offline scripts on the same machine as
5+// the page, so the two sets of timings share a clock and may honestly share
6+// an axis.
7+import { useRef, useState } from 'react'
8+import { Convergence, WorkPrecision, type SweepOut, type SweepStyle } from '../onedim/panels.tsx'
9+import { method as methodColor, type MethodKey } from '../../plot/palette.ts'
10+import type { Num } from '../../engine/types.ts'
11+
12+export interface BenchOut extends SweepOut {
13+ DOF: Num[][]
14+}
15+
16+interface OfflineMethod {
17+ name: string
18+ slot?: string
19+ ns: number[]
20+ E: Num[]
21+ T: Num[]
22+}
23+
24+interface OfflineFile {
25+ label: string
26+ kind?: string
27+ methods: OfflineMethod[]
28+}
29+
30+export const BENCH_SLOTS: MethodKey[] = ['fd2', 'mps', 'bie', 'spectral']
31+
32+/** Offline series keep the colour of their method and take a dash per file. */
33+const FILE_DASHES = ['6 4', '2 3', '9 3 2 3']
34+
35+function isSlot(s: string | undefined): s is MethodKey {
36+ return s != null && s in methodColor
37+}
38+
39+/** Accepts the documented schema and rejects everything else with a reason. */
40+function parseOffline(name: string, text: string): OfflineFile | string {
41+ let raw: unknown
42+ try {
43+ raw = JSON.parse(text)
44+ } catch {
45+ return `${name}: not JSON`
46+ }
47+ const o = raw as Partial<OfflineFile> & { schema?: string }
48+ if (o.schema !== 'laplace-dirichlet-bench-1') return `${name}: missing schema "laplace-dirichlet-bench-1"`
49+ if (!Array.isArray(o.methods) || o.methods.length === 0) return `${name}: no methods`
50+ for (const m of o.methods) {
51+ if (typeof m.name !== 'string' || !Array.isArray(m.ns) || !Array.isArray(m.E) || !Array.isArray(m.T))
52+ return `${name}: method entries need name, ns, E, T`
53+ }
54+ return { label: typeof o.label === 'string' ? o.label : name, kind: o.kind, methods: o.methods }
55+}
56+
57+/** The browser sweep merged with whatever offline files are loaded. */
58+function merged(data: BenchOut, files: OfflineFile[]): { out: SweepOut; style: SweepStyle } {
59+ const names = [...data.names]
60+ const NS: Num[][] = data.NS ? data.NS.map((r) => [...r]) : data.names.map(() => [...data.ns])
61+ const E = data.E.map((r) => [...r])
62+ const T = data.T.map((r) => [...r])
63+ const orders = data.orders.map((r) => [...r])
64+ const slots: MethodKey[] = [...BENCH_SLOTS]
65+ const dashes: (string | undefined)[] = data.names.map(() => undefined)
66+ files.forEach((f, i) => {
67+ for (const m of f.methods) {
68+ names.push(`${m.name} · ${f.label}`)
69+ NS.push(m.ns)
70+ E.push(m.E)
71+ T.push(m.T)
72+ orders.push(m.ns.map(() => null))
73+ slots.push(isSlot(m.slot) ? m.slot : 'fd2')
74+ dashes.push(FILE_DASHES[i % FILE_DASHES.length])
75+ }
76+ })
77+ return { out: { names, ns: data.ns, NS, E, T, orders }, style: { slots, dashes } }
78+}
79+
80+export function BenchSweep({ data, kind }: { data: BenchOut; kind: string }) {
81+ const [files, setFiles] = useState<OfflineFile[]>([])
82+ const [problem, setProblem] = useState<string | null>(null)
83+ const input = useRef<HTMLInputElement>(null)
84+
85+ const onFiles = async (list: FileList | null) => {
86+ if (!list) return
87+ const next: OfflineFile[] = []
88+ let bad: string | null = null
89+ for (const f of Array.from(list)) {
90+ const parsed = parseOffline(f.name, await f.text())
91+ if (typeof parsed === 'string') bad = parsed
92+ else next.push(parsed)
93+ }
94+ setProblem(bad)
95+ if (next.length > 0) setFiles((cur) => [...cur, ...next])
96+ if (input.current) input.current.value = ''
97+ }
98+
99+ const m = merged(data, files)
100+ return (
101+ <>
102+ <Convergence
103+ data={data}
104+ style={{ slots: BENCH_SLOTS, refs: [{ p: -4 / 3, label: 'slope −4/3' }] }}
105+ height={320}
106+ xLabel="the method's own resolution parameter n"
107+ />
108+ <WorkPrecision
109+ data={m.out}
110+ style={m.style}
111+ height={340}
112+ labelCount={data.names.length}
113+ resLabel="resolution"
114+ />
115+ <div className="lab-bar">
116+ <button className="ghost" onClick={() => input.current?.click()}>
117+ load offline results...
118+ </button>
119+ <input
120+ ref={input}
121+ type="file"
122+ accept=".json,application/json"
123+ multiple
124+ style={{ display: 'none' }}
125+ onChange={(e) => void onFiles(e.target.files)}
126+ />
127+ {files.map((f, i) => (
128+ <span className="stat" key={`${f.label}-${i}`}>
129+ <span className="stat-label">
130+ {f.label}
131+ {f.kind && f.kind !== kind ? ` (recorded on '${f.kind}')` : ''}
132+ </span>
133+ <button className="ghost" onClick={() => setFiles((cur) => cur.filter((_, j) => j !== i))}>
134+ remove
135+ </button>
136+ </span>
137+ ))}
138+ {problem && <span className="warn">{problem}</span>}
139+ </div>
140+ </>
141+ )
142+}
src/chapters/ch11/text.mdadded+44−0View file
@@ -0,0 +1,44 @@
1+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.
2+
3+The domain is the *Pac-Man*,
4+
5+$$\Omega = \{(r,\theta) : 0 < r < 1,\ 0 < \theta < 3\pi/2\},$$
6+
7+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.
8+
9+<!--fig:pacman-->
10+
11+### Four methods, four handicaps
12+
13+Each method is the one you have already met, pushed off its home ground.
14+
15+- **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.
16+- **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.
17+- **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.
18+- **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.
19+
20+<!--lab:field-->
21+
22+### The benchmark
23+
24+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.
25+
26+<!--lab:bench-->
27+
28+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.
29+
30+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.
31+
32+### Benchmarks beyond the browser
33+
34+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:
35+
36+```
37+python3 bench/python/bench_pacman.py
38+julia bench/julia/bench_pacman.jl
39+matlab -batch "run('bench/matlab/bench_pacman.m')"
40+```
41+
42+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.
43+
44+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.
src/chapters/index.tsmodified+2−1View file
@@ -10,5 +10,6 @@ import { ch07 } from './ch07/chapter.ts'
1010 import { ch08 } from './ch08/chapter.ts'
1111 import { ch09 } from './ch09/chapter.ts'
1212 import { ch10 } from './ch10/chapter.ts'
13+import { ch11 } from './ch11/chapter.ts'
1314
14-export const CHAPTERS: Chapter[] = [intro, ch01, ch02, ch03, ch04, ch05, ch06, ch07, ch08, ch09, ch10]
15+export const CHAPTERS: Chapter[] = [intro, ch01, ch02, ch03, ch04, ch05, ch06, ch07, ch08, ch09, ch10, ch11]
src/chapters/intro/chapter.tsmodified+1−1View file
@@ -6,7 +6,7 @@ export const intro: Chapter = {
66 nav: 'Introduction',
77 part: 'Introduction',
88 title: 'The Dirichlet problem',
9- blurb: 'What the problem is, what the tutorial measures, and how the ten chapters are arranged.',
9+ blurb: 'What the problem is, what the tutorial measures, and how the eleven chapters are arranged.',
1010 text,
1111 labs: [],
1212 }
src/chapters/intro/text.mdmodified+3−3View file
@@ -8,7 +8,7 @@ The equation by itself does not determine $u$. Adding any harmonic function to a
88
99 $$\Delta u = 0 \ \ \text{in } \Omega, \qquad u = g \ \ \text{on } \partial\Omega,$$
1010
11-or 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: nearly every technique used on harder elliptic problems can be introduced here, where the answer is known and nothing else is in the way.
11+or 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.
1212
1313 <!--fig:dirichlet-problem-->
1414
@@ -18,7 +18,7 @@ A tutorial on numerical methods has to decide what counts as success. Showing th
1818
1919 To 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.
2020
21-What 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 ten 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.
21+What 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.
2222
2323 <!--fig:accuracy-cost-->
2424
@@ -34,7 +34,7 @@ The geometry is what removes those special features, one at a time, so the chapt
3434
3535 **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).
3636
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).
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).
3838
3939 ### How the pages work
4040
src/chapters/onedim/panels.tsxmodified+132−48View file
@@ -1,11 +1,11 @@
11 // Figures for the one-dimensional labs. Every lab in Part I sends back one of
22 // these two shapes, so the same two panels serve all of them.
3-import { Fragment } from 'react'
4-import Plot from '../../plot/Plot.tsx'
3+import { Fragment, useState } from 'react'
4+import Plot, { type PointerPx } from '../../plot/Plot.tsx'
55 import Legend from '../../plot/Legend.tsx'
6-import { EndLabel, Line, Markers, RefSlope } from '../../plot/Curves.tsx'
6+import { EndLabel, Line, Markers, PointTip, RefSlope } from '../../plot/Curves.tsx'
77 import { ink, marker, method, type MethodKey } from '../../plot/palette.ts'
8-import { extent, padDomain } from '../../plot/scales.ts'
8+import { extent, padDomain, type Frame } from '../../plot/scales.ts'
99 import type { Num } from '../../engine/types.ts'
1010
1111 export interface SingleOut {
@@ -44,6 +44,13 @@ export function fmtErr(v: Num): string {
4444 return v.toExponential(2)
4545 }
4646
47+/** A time in milliseconds, at three significant figures either side of a second. */
48+export function fmtMs(v: Num): string {
49+ if (v == null || !isFinite(v)) return '--'
50+ if (v >= 1000) return `${(v / 1000).toPrecision(3)} s`
51+ return `${v.toPrecision(3)} ms`
52+}
53+
4754 /** The computed solution over the exact one, and the error underneath. */
4855 export function Solution1D({ data, slot = 'fd2' }: { data: SingleOut; slot?: MethodKey }) {
4956 const color = method[slot]
@@ -203,23 +210,64 @@ export function OrderTable({ data }: { data: SweepOut }) {
203210 )
204211 }
205212
213+/** The data point nearest the cursor, if one is close enough to have been meant. */
214+function nearestPoint(data: SweepOut, f: Frame, p: PointerPx | null): { k: number; j: number } | null {
215+ if (!p) return null
216+ let best: { k: number; j: number } | null = null
217+ let bestD = 26 * 26
218+ data.E.forEach((row, k) => {
219+ const t = data.T[k]
220+ for (let j = 0; j < row.length; j++) {
221+ const e = row[j]
222+ const tt = t[j]
223+ if (e == null || tt == null || !isFinite(e) || !isFinite(tt) || e <= 0 || tt <= 0) continue
224+ const dx = f.sx(tt) - p.px
225+ const dy = f.sy(e) - p.py
226+ const d = dx * dx + dy * dy
227+ if (d < bestD) {
228+ bestD = d
229+ best = { k, j }
230+ }
231+ }
232+ })
233+ return best
234+}
235+
206236 /** Error against wall clock: the only comparison that decides anything. */
207237 export function WorkPrecision({
208238 data,
209239 style,
210240 height = 320,
211241 label = 'milliseconds for one solve',
242+ labelCount,
243+ resLabel = 'N',
212244 }: {
213245 data: SweepOut
214246 style: SweepStyle
215247 height?: number
216248 label?: string
249+ /** how many series get a direct end-label; the rest rely on the legend */
250+ labelCount?: number
251+ /** what the resolution parameter is called in the hover box */
252+ resLabel?: string
217253 }) {
218254 const yd = seriesExtent(data.E)
219255 const xd = seriesExtent(data.T)
256+ // Two ways to ask which curve is which: point at a marker, which names its
257+ // series and reads out that point's numbers, or point at a legend entry,
258+ // which picks the curve out of the rest. Neither is decoration here: once
259+ // offline results are loaded, several series share a colour and differ only
260+ // in their dash pattern, and the legend alone cannot settle it.
261+ const [ptr, setPtr] = useState<PointerPx | null>(null)
262+ const [picked, setPicked] = useState<number | null>(null)
263+ const emph = (k: number, active: number | null) => (active == null || active === k ? 1 : 0.18)
264+
220265 return (
221266 <>
222- <div className="plot-title">accuracy against cost; down and to the left is better</div>
267+ <div className="plot-title">
268+ accuracy against cost; down and to the left is better. Point at a marker for its numbers, or at a legend
269+ entry for its curve.
270+ </div>
223271 <Plot
224272 height={height}
225273 xDomain={xd}
@@ -229,53 +277,87 @@ export function WorkPrecision({
229277 xLabel={label}
230278 yLabel="max |error|"
231279 margin={{ r: 30 }}
232- over={(f) => (
233- <>
234- {data.E.map((row, k) => {
235- // label the cheapest end of each curve, which is where the
236- // curves are furthest apart
237- const t = data.T[k]
238- let first = -1
239- for (let i = 0; i < row.length; i++) {
240- if (row[i] != null && t[i] != null) {
241- first = i
242- break
280+ onPointer={setPtr}
281+ over={(f) => {
282+ const hit = nearestPoint(data, f, ptr)
283+ const active = hit ? hit.k : picked
284+ return (
285+ <>
286+ {data.E.map((row, k) => {
287+ // label the cheapest end of each curve, which is where the
288+ // curves are furthest apart; alternate the vertical offset so
289+ // that labels of curves starting at similar cost do not
290+ // collide. A picked curve is labelled whether or not it is one
291+ // of the few that carry a label by default.
292+ if (active != null ? active !== k : labelCount != null && k >= labelCount) return null
293+ const t = data.T[k]
294+ let first = -1
295+ for (let i = 0; i < row.length; i++) {
296+ if (row[i] != null && t[i] != null) {
297+ first = i
298+ break
299+ }
243300 }
244- }
245- if (first < 0) return null
246- const px = f.sx(t[first] as number)
247- const right = px > f.iw * 0.65
248- return (
249- <EndLabel
250- key={k}
301+ if (first < 0) return null
302+ const px = f.sx(t[first] as number)
303+ const right = px > f.iw * 0.65
304+ return (
305+ <EndLabel
306+ key={k}
307+ f={f}
308+ x={t[first] as number}
309+ y={row[first] as number}
310+ text={data.names[k]}
311+ color={method[style.slots[k] ?? 'fd2']}
312+ anchor={right ? 'end' : 'start'}
313+ dx={right ? -8 : 8}
314+ dy={12 + (k % 2) * 14}
315+ />
316+ )
317+ })}
318+ {hit && (
319+ <PointTip
251320 f={f}
252- x={t[first] as number}
253- y={row[first] as number}
254- text={data.names[k]}
255- color={method[style.slots[k] ?? 'fd2']}
256- anchor={right ? 'end' : 'start'}
257- dx={right ? -8 : 8}
258- dy={12}
321+ px={f.sx(data.T[hit.k][hit.j] as number)}
322+ py={f.sy(data.E[hit.k][hit.j] as number)}
323+ color={method[style.slots[hit.k] ?? 'fd2']}
324+ lines={[
325+ data.names[hit.k],
326+ `${resLabel} = ${xsOf(data, hit.k)[hit.j] ?? '--'}`,
327+ `error ${fmtErr(data.E[hit.k][hit.j])}`,
328+ `${fmtMs(data.T[hit.k][hit.j])} for one solve`,
329+ ]}
259330 />
260- )
261- })}
262- </>
263- )}
331+ )}
332+ </>
333+ )
334+ }}
264335 >
265- {(f) => (
266- <>
267- {data.E.map((row, k) => {
268- const slot = style.slots[k] ?? 'fd2'
269- const t = data.T[k]
270- return (
271- <g key={k}>
272- <Line x={t.map((v) => v ?? NaN)} y={row} f={f} color={method[slot]} dash={style.dashes?.[k]} />
273- <Markers x={t.map((v) => v ?? NaN)} y={row} f={f} color={method[slot]} shape={marker[slot]} />
274- </g>
275- )
276- })}
277- </>
278- )}
336+ {(f) => {
337+ const active = nearestPoint(data, f, ptr)?.k ?? picked
338+ return (
339+ <>
340+ {data.E.map((row, k) => {
341+ const slot = style.slots[k] ?? 'fd2'
342+ const t = data.T[k]
343+ const o = emph(k, active)
344+ return (
345+ <g key={k} opacity={o}>
346+ <Line
347+ x={t.map((v) => v ?? NaN)}
348+ y={row}
349+ f={f}
350+ color={method[slot]}
351+ dash={style.dashes?.[k]}
352+ width={active === k ? 2.8 : 2}
353+ />
354+ <Markers x={t.map((v) => v ?? NaN)} y={row} f={f} color={method[slot]} shape={marker[slot]} />
355+ </g>
356+ )
357+ })}
358+ </>
359+ )
360+ }}
279361 </Plot>
280362 <Legend
281363 items={data.names.map((n, k) => ({
@@ -283,6 +365,8 @@ export function WorkPrecision({
283365 color: method[style.slots[k] ?? 'fd2'],
284366 dash: style.dashes?.[k],
285367 }))}
368+ onHover={setPicked}
369+ active={picked}
286370 />
287371 </>
288372 )
src/figures/beyond.tsxmodified+88−0View file
@@ -225,3 +225,91 @@ export function CloseEvaluation() {
225225 </Figure>
226226 )
227227 }
228+
229+/** The Pac-Man domain, and its image as a rectangle in polar coordinates. */
230+export function PacmanDomain() {
231+ const CX = 158
232+ const CY = 128
233+ const R = 88
234+ // the rectangle: theta across, r up
235+ const RX = 430
236+ const RY = 48
237+ const RW = 240
238+ const RH = 150
239+ return (
240+ <Figure
241+ width={716}
242+ height={252}
243+ caption="The benchmark domain: the unit disk with the quadrant x > 0, y < 0 removed. The boundary is one smooth arc and two straight edges, with right-angle corners where they meet the arc and a re-entrant corner of angle 3π/2 at the origin. In polar coordinates the same domain is a rectangle, on which the whole side r = 0 is the single corner point; that is the rectangle the fourth method collocates on, and the reason the corner singularity sits on the edge of its expansion domain."
244+ >
245+ {/* the domain */}
246+ <path
247+ d={`M ${CX} ${CY} L ${CX + R} ${CY} A ${R} ${R} 0 1 0 ${CX} ${CY + R} Z`}
248+ fill="#20242c"
249+ stroke="none"
250+ />
251+ <line x1={CX} y1={CY} x2={CX + R} y2={CY} stroke={method.fd2} strokeWidth={2.4} />
252+ <path d={`M ${CX + R} ${CY} A ${R} ${R} 0 1 0 ${CX} ${CY + R}`} fill="none" stroke={method.bie} strokeWidth={2.4} />
253+ <line x1={CX} y1={CY + R} x2={CX} y2={CY} stroke={method.fast} strokeWidth={2.4} />
254+ <Label x={CX - 34} y={CY - 40} size={13} fill={ink.secondary}>
255+ Ω
256+ </Label>
257+ <Label x={CX} y={16} size={12.5} weight={600} fill={ink.primary}>
258+ the Pac-Man domain
259+ </Label>
260+
261+ {/* the interior angle at the corner */}
262+ <path d={`M ${CX + 22} ${CY} A 22 22 0 1 0 ${CX} ${CY + 22}`} fill="none" stroke={method.mps} strokeWidth={1.6} />
263+ <Node x={CX} y={CY} r={3.6} fill={method.mps} />
264+ <Label x={CX - 36} y={CY + 40} size={12} fill={method.mps}>
265+ 3π/2
266+ </Label>
267+
268+ {/* right-angle marks where the edges meet the arc */}
269+ <path
270+ d={`M ${CX + R - 11} ${CY} v 11 h 11`}
271+ fill="none"
272+ stroke={ink.muted}
273+ strokeWidth={1.2}
274+ />
275+ <path
276+ d={`M ${CX} ${CY + R - 11} h 11 v 11`}
277+ fill="none"
278+ stroke={ink.muted}
279+ strokeWidth={1.2}
280+ />
281+ <Label x={CX + R + 14} y={CY + 16} size={11.5} anchor="start" fill={ink.muted}>
282+ π/2
283+ </Label>
284+ <Label x={CX + 16} y={CY + R + 14} size={11.5} anchor="start" fill={ink.muted}>
285+ π/2
286+ </Label>
287+
288+ {/* the polar rectangle */}
289+ <Label x={RX + RW / 2} y={16} size={12.5} weight={600} fill={ink.primary}>
290+ the same domain in (θ, r)
291+ </Label>
292+ <rect x={RX} y={RY} width={RW} height={RH} fill="#20242c" stroke="none" />
293+ <line x1={RX} y1={RY} x2={RX + RW} y2={RY} stroke={method.bie} strokeWidth={2.4} />
294+ <line x1={RX + RW} y1={RY} x2={RX + RW} y2={RY + RH} stroke={method.fd2} strokeWidth={2.4} />
295+ <line x1={RX} y1={RY + RH} x2={RX} y2={RY} stroke={method.fast} strokeWidth={2.4} />
296+ <line x1={RX} y1={RY + RH} x2={RX + RW} y2={RY + RH} stroke={method.mps} strokeWidth={2.4} strokeDasharray="3 4" />
297+ <Node x={RX + RW / 2} y={RY + RH} r={3.6} fill={method.mps} />
298+ <Label x={RX + RW / 2} y={RY - 10} size={11.5} fill={method.bie}>
299+ r = 1: the arc
300+ </Label>
301+ <Label x={RX + RW + 10} y={RY + RH / 2} size={11.5} anchor="start" fill={method.fd2}>
302+ θ = 0
303+ </Label>
304+ <Label x={RX - 10} y={RY + RH / 2} size={11.5} anchor="end" fill={method.fast}>
305+ θ = 3π/2
306+ </Label>
307+ <Label x={RX + RW / 2} y={RY + RH + 16} size={11.5} fill={method.mps}>
308+ r = 0: the whole side is the corner
309+ </Label>
310+ <Label x={RX + RW / 2} y={RY + RH + 34} size={11.5} fill={ink.muted} italic>
311+ u here is one known number, g(0, 0)
312+ </Label>
313+ </Figure>
314+ )
315+}
src/figures/index.tsmodified+2−1View file
@@ -4,7 +4,7 @@ import type { ComponentType } from 'react'
44 import { AccuracyCost, DirichletProblem, Domains } from './intro.tsx'
55 import { ChebyshevPoints, CompactScheme, ManufacturedLoop, RoundoffFloor, ThreePointStencil } from './onedim.tsx'
66 import { FastPoissonSteps, FivePointStencil, NinePointStencil, SmoothingIdea, VCycle } from './twodim.tsx'
7-import { AreaAgainstBoundary, CloseEvaluation, ReentrantCorner } from './beyond.tsx'
7+import { AreaAgainstBoundary, CloseEvaluation, PacmanDomain, ReentrantCorner } from './beyond.tsx'
88
99 export const FIGURES: Record<string, ComponentType> = {
1010 'dirichlet-problem': DirichletProblem,
@@ -23,4 +23,5 @@ export const FIGURES: Record<string, ComponentType> = {
2323 'reentrant-corner': ReentrantCorner,
2424 'area-boundary': AreaAgainstBoundary,
2525 'close-evaluation': CloseEvaluation,
26+ pacman: PacmanDomain,
2627 }
src/figures/intro.tsxmodified+24−3View file
@@ -120,15 +120,15 @@ export function DirichletProblem() {
120120 )
121121 }
122122
123-/** The four domains, in the order the tutorial meets them. */
123+/** The five domains, in the order the tutorial meets them. */
124124 export function Domains() {
125125 const c = ink.node
126- const box = (i: number) => 12 + i * 176
126+ const box = (i: number) => 2 + i * 143
127127 return (
128128 <Figure
129129 width={716}
130130 height={196}
131- caption="The four domains. Each step to the right removes a structural assumption: the second dimension, then the convexity of the corner, then the straight edge itself."
131+ caption="The five domains. Each step to the right removes a structural assumption: the second dimension, then the convexity of the corner, then the straight edge itself; the last, for the closing benchmark, has the curved boundary and the re-entrant corner at once."
132132 >
133133 {/* 1-D interval */}
134134 <g transform={`translate(${box(0)} 0)`}>
@@ -207,6 +207,27 @@ export function Domains() {
207207 Chapter 10
208208 </Label>
209209 </g>
210+
211+ {/* the Pac-Man */}
212+ <g transform={`translate(${box(4)} 0)`}>
213+ <Label x={76} y={16} size={12.5} weight={600} fill={ink.primary}>
214+ Pac-Man
215+ </Label>
216+ <path
217+ d="M 76 92 L 130 92 A 54 54 0 1 0 76 146 Z"
218+ fill="#20242c"
219+ stroke={method.mps}
220+ strokeWidth={2.2}
221+ strokeLinejoin="round"
222+ />
223+ <circle cx={76} cy={92} r={5} fill="none" stroke={method.mps} strokeWidth={2} />
224+ <Label x={95} y={106} size={11} anchor="start" fill={method.mps}>
225+ 3π/2
226+ </Label>
227+ <Label x={76} y={168} size={11.5} fill={ink.muted}>
228+ Chapter 11
229+ </Label>
230+ </g>
210231 </Figure>
211232 )
212233 }
src/index.cssmodified+4−0View file
@@ -460,6 +460,10 @@ button.ghost:hover {
460460 gap: 5px;
461461 }
462462
463+.legend-item.dim {
464+ opacity: 0.4;
465+}
466+
463467 .legend-value {
464468 color: var(--muted);
465469 font-variant-numeric: tabular-nums;
src/matlab/lp_gauss.madded+25−0View file
@@ -0,0 +1,25 @@
1+function [x, w] = lp_gauss(p)
2+% LP_GAUSS Gauss-Legendre nodes and weights on [-1, 1].
3+%
4+% Newton's method on the Legendre polynomial P_p, starting from the Chebyshev
5+% estimate of the roots; three or four iterations reach roundoff for the small
6+% p used here. The weights are 2 / ((1 - x^2) P_p'(x)^2).
7+x = cos(pi * ((1:p)' - 0.25) / (p + 0.5));
8+dp = ones(p, 1);
9+for it = 1:50
10+ p0 = ones(p, 1);
11+ p1 = x;
12+ for k = 2:p
13+ p2 = ((2*k - 1) * x .* p1 - (k - 1) * p0) / k;
14+ p0 = p1;
15+ p1 = p2;
16+ end
17+ dp = p * (x .* p1 - p0) ./ (x.^2 - 1);
18+ dx = p1 ./ dp;
19+ x = x - dx;
20+ if max(abs(dx)) < 1e-15
21+ break;
22+ end
23+end
24+w = 2 ./ ((1 - x.^2) .* dp.^2);
25+end
src/plot/Curves.tsxmodified+61−0View file
@@ -153,6 +153,67 @@ export function RefSlope({
153153 )
154154 }
155155
156+/**
157+ * A ring around the data point under the cursor, and a box of its numbers.
158+ *
159+ * Drawn as SVG rather than an HTML overlay, like everything else in this
160+ * layer, so it needs no positioning relative to the host element. The box
161+ * width is estimated from the character count, which is why the lines are kept
162+ * short; it flips to the other side of the point rather than leave the frame.
163+ */
164+export function PointTip({
165+ f,
166+ px,
167+ py,
168+ color,
169+ lines,
170+}: {
171+ f: Frame
172+ px: number
173+ py: number
174+ color: string
175+ /** first line is the series, in its own colour; the rest are the numbers */
176+ lines: string[]
177+}) {
178+ const size = 11.5
179+ const pad = 6
180+ const lh = 15
181+ const w = Math.max(...lines.map((l) => l.length)) * 6.15 + 2 * pad
182+ const h = lines.length * lh + 2 * pad - 3
183+ const right = px + 14 + w < f.iw
184+ const bx = Math.min(Math.max(right ? px + 14 : px - 14 - w, 2), Math.max(2, f.iw - w - 2))
185+ const by = Math.min(Math.max(py - h / 2, 2), Math.max(2, f.ih - h - 2))
186+ return (
187+ <g pointerEvents="none">
188+ <circle cx={px} cy={py} r={7} fill="none" stroke={color} strokeWidth={2} />
189+ <circle cx={px} cy={py} r={7} fill="none" stroke={surface} strokeWidth={1} opacity={0.7} />
190+ <rect
191+ x={bx}
192+ y={by}
193+ width={w}
194+ height={h}
195+ rx={4}
196+ fill={surface}
197+ fillOpacity={0.94}
198+ stroke={ink.axis}
199+ strokeWidth={1}
200+ />
201+ {lines.map((l, i) => (
202+ <text
203+ key={i}
204+ x={bx + pad}
205+ y={by + pad + lh * i + size}
206+ fill={i === 0 ? color : ink.secondary}
207+ fontSize={size}
208+ style={i === 0 ? undefined : { fontVariantNumeric: 'tabular-nums' }}
209+ >
210+ {l}
211+ </text>
212+ ))}
213+ </g>
214+ )
215+}
216+
156217 /** A label placed at the end of a curve, so the reader need not use the legend. */
157218 export function EndLabel({
158219 f,
src/plot/Legend.tsxmodified+19−4View file
@@ -15,11 +15,26 @@ export interface LegendItem {
1515 * plot appears here too, and the ones the reader most needs are also labelled
1616 * directly on the plot.
1717 */
18-export default function Legend({ items }: { items: LegendItem[] }) {
18+export default function Legend({
19+ items,
20+ onHover,
21+ active,
22+}: {
23+ items: LegendItem[]
24+ /** hovering an entry reports its index, so the plot can pick that curve out */
25+ onHover?: (k: number | null) => void
26+ /** the entry to show as picked out; the others are dimmed */
27+ active?: number | null
28+}) {
1929 return (
20- <div className="legend">
21- {items.map((it) => (
22- <span className="legend-item" key={it.label}>
30+ <div className="legend" onPointerLeave={onHover ? () => onHover(null) : undefined}>
31+ {items.map((it, k) => (
32+ <span
33+ className={active != null && active !== k ? 'legend-item dim' : 'legend-item'}
34+ key={it.label}
35+ onPointerEnter={onHover ? () => onHover(k) : undefined}
36+ style={onHover ? { cursor: 'default' } : undefined}
37+ >
2338 <svg width={18} height={10} aria-hidden="true">
2439 {it.shape === 'dot' ? (
2540 <circle cx={9} cy={5} r={3.5} fill={it.color} stroke={ink.grid} strokeWidth={1} />
src/plot/Plot.tsxmodified+28−3View file
@@ -26,6 +26,19 @@ export interface PlotProps {
2626 hoverX?: number | null
2727 /** rendered as an absolutely positioned box over the plot */
2828 overlay?: ReactNode
29+ /**
30+ * Pointer position in frame pixels, or null when it leaves. The caller is
31+ * expected to keep it in state and use the frame it is handed in `over` to
32+ * turn it back into data coordinates, which is how a plot picks the point
33+ * nearest the cursor without either side owning the scales twice.
34+ */
35+ onPointer?: (p: { px: number; py: number } | null) => void
36+}
37+
38+/** Pointer position in frame pixels: what `onPointer` reports. */
39+export interface PointerPx {
40+ px: number
41+ py: number
2942 }
3043
3144 const DEFAULT_MARGIN = { l: 54, r: 14, t: 10, b: 30 }
@@ -66,6 +79,7 @@ export default function Plot(props: PlotProps) {
6679 onHoverX,
6780 hoverX,
6881 overlay,
82+ onPointer,
6983 } = props
7084 const m = { ...DEFAULT_MARGIN, ...props.margin }
7185 const [hostRef, width] = useWidth()
@@ -101,11 +115,21 @@ export default function Plot(props: PlotProps) {
101115 const yt = props.yTicks ?? (yLog ? decadeTicks(yDomain[0], yDomain[1]) : niceTicks(yDomain[0], yDomain[1], 5))
102116
103117 const pointer = (e: React.PointerEvent) => {
104- if (!onHoverX) return
118+ if (!onHoverX && !onPointer) return
105119 const rect = svgRef.current?.getBoundingClientRect()
106120 if (!rect) return
107121 const px = e.clientX - rect.left - m.l
108- onHoverX(px < -4 || px > iw + 4 ? null : f.ix(Math.min(iw, Math.max(0, px))))
122+ const py = e.clientY - rect.top - m.t
123+ onHoverX?.(px < -4 || px > iw + 4 ? null : f.ix(Math.min(iw, Math.max(0, px))))
124+ // a little slack outside the frame, so a marker sitting on an axis is
125+ // still reachable
126+ const out = px < -10 || px > iw + 10 || py < -10 || py > ih + 10
127+ onPointer?.(out ? null : { px, py })
128+ }
129+
130+ const leave = () => {
131+ onHoverX?.(null)
132+ onPointer?.(null)
109133 }
110134
111135 return (
@@ -116,7 +140,8 @@ export default function Plot(props: PlotProps) {
116140 height={height}
117141 role="img"
118142 onPointerMove={pointer}
119- onPointerLeave={() => onHoverX?.(null)}
143+ onPointerLeave={leave}
144+ style={onPointer ? { cursor: 'crosshair' } : undefined}
120145 >
121146 <defs>
122147 <clipPath id={clipId}>