An interactive tutorial on the Dirichlet problem, accuracy against cost
Ten chapters, in one and two dimensions, from the three-point difference
scheme to a double-layer boundary integral equation on a curved domain. Every
problem has an exact solution known by construction, so every error is exact
and every claim is measured rather than asserted; each method is a MATLAB
script the reader can edit and run in the page through numbl.
The drivers are written around numbl 0.4.18's compiler: sweeps are recursions
rather than loops, no square root is taken in MATLAB where a loop can see it,
and the multigrid hierarchy lives in a flat vector. CLAUDE.md records why, with
the measurements behind each rule.
86 changed files+11629−0
.gitignoreadded+5−0View file
@@ -0,0 +1,5 @@
1+node_modules
2+dist
3+.DS_Store
4+*.local
5+tmp
CLAUDE.mdadded+180−0View file
@@ -0,0 +1,180 @@
1+# CLAUDE.md
2+
3+Notes for anyone, human or agent, working in this repository. Read
4+[README.md](README.md) first for what the tutorial contains.
5+
6+## The shape of the thing
7+
8+A Vite + React app that embeds [numbl](https://numbl.org) as a library
9+(`numbl/browser`, `createNumblSession`). This is the seqlab / barycentric-rational
10+pattern, not the `numbl-project.json` + site-viewer pattern: there is no project
11+file and the deploy is a plain `vite build`.
12+
13+```
14+src/
15+ App.tsx sidebar, hash routing, chapter switching
16+ ui/ Lab (one runnable section), ChapterView, Controls, Markdown
17+ engine/ the numbl session wrapper
18+ plot/ a small hand-rolled SVG layer, plus a canvas field plot
19+ matlab/ the shared .m library, shown to the reader on its own page
20+ chapters/
21+ onedim/driver.m the driver for chapters 1-4, plus their shared panels
22+ twodim/driver.m the driver for chapters 5-7, plus their shared panels
23+ chNN/ per-chapter text.md, editable .m script, chapter.ts
24+```
25+
26+A **chapter** is prose in one `text.md` with `<!--lab:id-->` markers where its labs
27+go. A **lab** is a driver (fixed, hidden), an editable script (the method), a list of
28+declarative controls, and a render function. Each lab owns one numbl session for as
29+long as it is on screen: moving a control rewrites `params.json` and re-runs, while
30+editing the script reboots the session, because the script is compiled into `main.m`
31+at boot.
32+
33+`main.m` is `driver.m` + `"\n"` + the editable script, so the script's functions
34+become **local functions of the driver** and the driver can call them. That is the
35+whole of the contract; each driver's header comment states which functions its
36+script must define.
37+
38+## Gotchas in numbl 0.4.18
39+
40+These were all found the hard way, by measurement. Most of them cost a factor of
41+twenty in speed rather than producing a wrong answer, which makes them easy to miss
42+in a project whose subject is speed.
43+
44+### Invoke the script with `run('main.m')`, not `main;`
45+
46+numbl mis-binds the arguments of a script's local functions when the script is
47+invoked by name from the REPL, which is what `session.execute` gives us: the callee
48+sees its parameters as undefined. `engine/runner.ts` goes through `run()` for this
49+reason.
50+
51+### The compiler's fast path, and the eight ways to fall off it
52+
53+numbl compiles MATLAB to JavaScript and infers types as it goes. A loop whose body
54+it cannot type is run through a generic path that is 20 to 70 times slower. The
55+following all cause that, and none of them are errors:
56+
57+1. **Cell or struct indexing inside a loop.** `C{k}` or `p.tol` in a loop body
58+ deoptimises the whole loop *and every function it calls*. Hoist the read above
59+ the loop.
60+2. **A multiple-output call inside a loop.** `[a, b] = f(...)`. Return one matrix
61+ and unpack it afterwards.
62+3. **A slice assignment inside a loop.** `A(k, :) = row`. Do the copy inside a
63+ called function, or with a scalar loop.
64+4. **Bare `tic` / `toc` inside a loop.** Use handles: `t = tic; ...; toc(t)`.
65+5. **Growing an array inside a loop**, that is, `v(k) = ...` for a `v` that was
66+ never preallocated.
67+6. **A loop variable with more than one kind of assignment.** This is the subtle
68+ one. `Y = V; for ... Y = f(Y); end` is slow, and `Y = f(V); for ... Y = f(Y); end`
69+ is fast, because in the first case the type of `Y` at the top of the loop is the
70+ merge of "a parameter" and "whatever `f` returns". It applies to dead branches
71+ too: an `if` arm that assigns `Y = reshape(...)` slows the loop down even when it
72+ never runs. Give each loop variable one source, or use separate names.
73+7. **Passing a loop-written array into a function call, even after the loop.**
74+ `for ... res(k) = ...; end; H = pack(res, ...)` is slow, while `H = [res; ...]`
75+ is fast. This is why the iterative chapter hands the page its *squared*
76+ residuals and lets JavaScript take the square roots.
77+8. **A call inside a loop, to a function that itself loops and calls further
78+ functions.** The chain "loop → solver → kernel with loops" is compiled
79+ generically even when every individual piece is fine. Recursion in place of the
80+ outer loop is compiled normally, which is what `drv_res`, `drv_methods` and their
81+ relatives in the drivers are for. This one was worth a factor of 30 on the
82+ L-shaped conjugate gradient solve.
83+
84+The rule of thumb that follows: **hot loops should touch only preallocated numeric
85+arrays and scalars, nothing they write should later be handed to a function, and any
86+loop that drives a timed solve should be a recursion.** When in doubt, measure: time
87+the same call inside and outside a loop and compare.
88+
89+### `sqrt`, `log` and `abs` widen a scalar's type to complex
90+
91+`v = v + r*r` in a loop followed by `v = sqrt(v)` after it makes the loop 20 times
92+slower, because `sqrt` in numbl can return a complex number and the inferred type of
93+`v` is widened accordingly. Accumulate into one variable and take the root into
94+another, or better, return the sum of squares and take the root in the caller. Every
95+`res_sumsq` in this repository exists for that reason.
96+
97+### `sqrt(NaN)` is `0 + NaN i`, and `log(NaN)` is complex
98+
99+MATLAB returns `NaN`. numbl returns a complex value whose **real part is zero**, so
100+`real(sqrt(NaN))` is `0`, not `NaN`. Square-rooting an array with a NaN tail silently
101+turns the tail into zeros, which in an error history reads as the iteration having
102+converged perfectly. The iterative chapter therefore never takes a square root in
103+MATLAB at all, and `lp_orders` guards its ratio before taking a logarithm.
104+
105+### `jsondecode` of a JSON array of arrays gives a cell of column vectors
106+
107+MATLAB gives a matrix when the rows are the same length. The L-shaped driver unpacks
108+`p.nsm` by hand.
109+
110+### `jsonencode` refuses complex values, and flattens single-row matrices
111+
112+`jsonencode([1 2 3])` is `[1,2,3]`, not `[[1,2,3]]`, so anything matrix-shaped goes
113+through `lp_rows`, which converts it to a cell array of rows first. NaN and Inf cross
114+as `null`, which is why the TypeScript side uses `Num = number | null` everywhere and
115+every plot helper skips nulls.
116+
117+### Backslash densifies sparse matrices
118+
119+`A \ b` converts a sparse `A` to dense before calling LAPACK, so there is no sparse
120+direct solver here. That is why Chapter 5 stops at $N = 48$ and why Chapters 6 to 10
121+are about methods that do not need one. The text says so where it matters.
122+
123+### Scalar loops are often faster than vectorised slices
124+
125+A hand-written double loop over a stencil beats the equivalent `A(2:end-1,...)`
126+expression by about a factor of three, because the vectorised form allocates several
127+temporaries per sweep. The kernels in `ch08/iterative.m` are written as loops for
128+that reason, which also happens to make them readable.
129+
130+### The browser is slower than the CLI
131+
132+The CLI uses a native LAPACK addon; the browser falls back to a TypeScript
133+implementation, and it announces which one on `console.error`. Dense solves are
134+noticeably slower in the page than in `npm run test:matlab`. Nothing in the tutorial
135+depends on the absolute numbers, but do not quote a CLI timing as if it were what a
136+reader will see.
137+
138+## Measuring time honestly
139+
140+Timings come from `tic`/`toc` **inside** the MATLAB layer, around the solve only, so
141+they exclude session boot, JSON encoding and rendering. Each method gets a warm-up
142+call before its first timed solve, because the first call also pays for compilation,
143+and any solve under 5 ms is repeated a few times with the minimum kept. Below about
144+0.1 ms the measurement is dominated by fixed overhead and the work-precision curves
145+flatten out on the left; that is the clock, not the method.
146+
147+## Testing
148+
149+```
150+npm run test:matlab # ~2 min: the .m layer through the numbl CLI
151+npm run build && npm run test:browser [--shots]
152+```
153+
154+Every assertion in `test:matlab` is a number stated in the prose: an observed order,
155+an error at a given resolution, a ratio between two methods. If one fails, either the
156+code or the text is wrong, and the text is not allowed to win by default.
157+
158+`test:browser` opens every chapter, presses run on the labs that wait to be asked,
159+and checks that nothing throws and that the figures contain data. It cannot judge
160+whether a picture is any good; `--shots` writes full-page screenshots into `tmp/` so
161+that a human, or an agent with vision, can look.
162+
163+## Plotting
164+
165+`src/plot/` is a small SVG layer with no chart library. `palette.ts` records which
166+colour does which job and the validator run that justifies the set; if you add a
167+series, re-run the validator rather than picking a hex that looks nice. Fields are
168+drawn into a canvas at one pixel per grid cell and scaled, so a 512-by-512 grid costs
169+no more to draw than a 16-by-16 one, and points outside the domain are passed as
170+`null` and left as background.
171+
172+## Things deliberately not done
173+
174+- A sparse direct solver, for the reason above.
175+- Graded meshes in Chapter 9. The chapter names local refinement and singular
176+ enrichment as the standard remedies and then demonstrates a third, because the
177+ third is the one that is short enough to read.
178+- Close-evaluation quadrature in Chapter 10. The failure is shown and quantified; the
179+ fix is a chapter of its own.
180+- Neumann and mixed boundary conditions, anywhere.
README.mdadded+97−0View file
@@ -0,0 +1,97 @@
1+# The Dirichlet problem: accuracy against cost
2+
3+An interactive tutorial on solving the Dirichlet problem for Laplace's and Poisson's
4+equations, in one and two dimensions, with every method written as a MATLAB-syntax
5+script that you can edit and run in the page. It runs in your browser through
6+[numbl](https://numbl.org); there is no server, nothing to install, and nothing is
7+sent anywhere.
8+
9+The organising idea is that every problem in it has an exact solution known in
10+closed form, by construction, so the error of every computation is known at every
11+point to full precision. What the tutorial then measures is not whether a method
12+works but what accuracy it delivers for what cost, and under which conditions that
13+answer changes.
14+
15+## The problems
16+
17+In one dimension, Laplace's equation is $u'' = 0$ and its solution is a straight
18+line, so Part I uses Poisson's equation $-u'' = f$ with the source manufactured from
19+a chosen solution. Five solutions are available, differing in how hard they are to
20+resolve: an entire function of moderate amplitude, an oscillatory one, a narrow
21+Gaussian, a boundary layer, and one that is $C^2$ but not $C^3$.
22+
23+In two dimensions the equation is Laplace's, and the exact solutions are harmonic
24+functions written down in closed form: $e^x \sin y$, $\operatorname{Re}(x+iy)^k$,
25+$e^{ky}\sin kx$, and $\log|z - z_0|$ with the singularity just outside the domain.
26+Their boundary traces are the data, so the only error in the computation is the
27+error of the method.
28+
29+## The chapters
30+
31+**One dimension**
32+
33+1. **Ground truth by construction** — manufactured solutions, the three-point
34+ scheme, the tridiagonal solve, and what the error field looks like.
35+2. **Second order, to its limit** — the condition number $\kappa \approx 0.4N^2$,
36+ the roundoff floor near $10^{-8.5}$, and Richardson extrapolation.
37+3. **Two more orders, almost free** — Numerov's compact fourth-order scheme, which
38+ changes only the right-hand side, and how to read a work-precision diagram.
39+4. **Chebyshev collocation** — a dense matrix, geometric convergence to $10^{-13}$
40+ with 64 points, and the two situations that take that away.
41+
42+**Two dimensions**
43+
44+5. **The five-point Laplacian** — harmonic ground truth, second order, and the
45+ $O(N^6)$ wall of a dense solve.
46+6. **The fast Poisson solver** — the same system diagonalised by the discrete sine
47+ transform: a million unknowns in half a second, identical answers.
48+7. **Sixth order for the same price** — the compact nine-point stencil, which is
49+ sixth-order accurate on Laplace's equation and uses the same fast solver, and
50+ two-dimensional Chebyshev collocation by matrix decomposition.
51+8. **Iterating instead of solving** — Jacobi, Gauss-Seidel, SOR, conjugate
52+ gradients and a multigrid V-cycle on one problem, and the question of when to
53+ stop.
54+
55+**Beyond the square**
56+
57+9. **A corner, and what it costs** — the L-shaped domain, the loss of second order
58+ to the $r^{2/3}$ singularity, and the method of particular solutions, which puts
59+ the singularity into the basis and reaches $10^{-14}$.
60+10. **Boundary integral equations** — a smooth curved domain, the double-layer
61+ potential, a second-kind Fredholm equation, the spectral accuracy of the
62+ periodic trapezoid rule, and the close-evaluation problem that comes with it.
63+
64+A **shared MATLAB library** page lists the small function files every lab boots
65+with.
66+
67+## Using it
68+
69+Each lab has controls above the figures and, behind the "show script" button, the
70+method itself. Edit the script and press run: the session reboots with your version
71+compiled in, and the figures are whatever your code produces. Ctrl-Enter (Cmd-Enter)
72+runs from inside the editor, and "reset" restores the original.
73+
74+Some labs run as soon as the chapter opens; the heavier ones wait for the run
75+button, which the caption says.
76+
77+## Development
78+
79+```
80+npm install
81+npm run dev # vite dev server
82+npm run build # typecheck and bundle to dist/
83+npm run test:matlab # the .m layer, through the numbl CLI (needs a numbl checkout)
84+npm run test:browser # the built app in headless Chrome
85+```
86+
87+`test:matlab` needs a numbl clone; set `NUMBL_DIR` if it is not at `~/src/numbl`.
88+Every assertion in it is a number claimed somewhere in the prose, so a failure means
89+either the code or the text is wrong. `test:browser` needs `npm run build` first, and
90+takes `--shots` to write a full-page screenshot of each chapter into `tmp/`.
91+
92+See [CLAUDE.md](CLAUDE.md) for how the pieces fit together, and in particular for the
93+numbl performance rules that the drivers are written around.
94+
95+## Licence
96+
97+Apache-2.0.
index.htmladded+17−0View file
@@ -0,0 +1,17 @@
1+<!doctype html>
2+<html lang="en">
3+ <head>
4+ <meta charset="UTF-8" />
5+ <link rel="icon" type="image/svg+xml" href="./favicon.svg" />
6+ <meta name="viewport" content="width=device-width, initial-scale=1.0" />
7+ <meta
8+ name="description"
9+ content="An interactive tutorial on the Dirichlet problem for Laplace's equation in one and two dimensions: manufactured solutions with exact ground truth, finite differences, spectral collocation, fast Poisson solvers, multigrid, corner singularities and boundary integral equations, measured for accuracy against cost. Every method is an editable MATLAB-syntax script that runs in your browser."
10+ />
11+ <title>The Dirichlet problem: accuracy and cost</title>
12+ </head>
13+ <body>
14+ <div id="root"></div>
15+ <script type="module" src="/src/main.tsx"></script>
16+ </body>
17+</html>
package-lock.jsonadded+4952−0View file
This diff is 4,957 lines long and is not shown.
package.jsonadded+38−0View file
@@ -0,0 +1,38 @@
1+{
2+ "name": "laplace-dirichlet",
3+ "private": true,
4+ "version": "0.1.0",
5+ "type": "module",
6+ "description": "An interactive tutorial on solving Laplace and Poisson Dirichlet problems in 1-D and 2-D, with every method written as an editable MATLAB-syntax script that runs in the browser through numbl",
7+ "license": "Apache-2.0",
8+ "scripts": {
9+ "dev": "vite",
10+ "build": "tsc -b && vite build",
11+ "preview": "vite preview",
12+ "typecheck": "tsc -b",
13+ "test:matlab": "node scripts/matlab-test.mjs",
14+ "test:browser": "node scripts/browser-test.mjs"
15+ },
16+ "dependencies": {
17+ "@codemirror/language": "^6.11.3",
18+ "@codemirror/legacy-modes": "^6.5.2",
19+ "@codemirror/state": "^6.5.2",
20+ "@codemirror/theme-one-dark": "^6.1.3",
21+ "@codemirror/view": "^6.38.6",
22+ "codemirror": "^6.0.2",
23+ "katex": "^0.16.11",
24+ "marked": "^15.0.0",
25+ "numbl": "^0.4.18",
26+ "react": "^19.2.0",
27+ "react-dom": "^19.2.0"
28+ },
29+ "devDependencies": {
30+ "@types/katex": "^0.16.7",
31+ "@types/react": "^19.2.5",
32+ "@types/react-dom": "^19.2.3",
33+ "@vitejs/plugin-react": "^5.1.1",
34+ "puppeteer": "^24.0.0",
35+ "typescript": "~5.9.3",
36+ "vite": "^7.2.4"
37+ }
38+}
public/favicon.svgadded+10−0View file
@@ -0,0 +1,10 @@
1+<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 64 64">
2+ <rect width="64" height="64" rx="12" fill="#1a1a19" />
3+ <!-- the square domain, with its boundary data marked on one side -->
4+ <rect x="12" y="12" width="40" height="40" rx="2" fill="none" stroke="#3c3c38" stroke-width="3" />
5+ <path d="M12 12 H52" stroke="#c98500" stroke-width="4" stroke-linecap="round" />
6+ <!-- level sets of a harmonic function, hottest near that side -->
7+ <path d="M12 21 Q32 29 52 21" fill="none" stroke="#3987e5" stroke-width="3" stroke-linecap="round" />
8+ <path d="M12 31 Q32 41 52 31" fill="none" stroke="#3987e5" stroke-width="3" stroke-linecap="round" opacity="0.72" />
9+ <path d="M12 41 Q32 49 52 41" fill="none" stroke="#3987e5" stroke-width="3" stroke-linecap="round" opacity="0.45" />
10+</svg>
scripts/browser-test.mjsadded+163−0View file
@@ -0,0 +1,163 @@
1+#!/usr/bin/env node
2+// Drives the built app in a headless browser: every chapter is opened, every
3+// lab is run, and the page is checked for console errors and for figures that
4+// actually contain data.
5+//
6+// npm run build && node scripts/browser-test.mjs [--shots] [--only <id>]
7+//
8+// --shots writes a full-page screenshot per chapter to tmp/, which is the only
9+// way to judge layout; the assertions here cover behaviour, not appearance.
10+import { spawn } from 'node:child_process'
11+import { mkdirSync } from 'node:fs'
12+import { dirname, join } from 'node:path'
13+import { fileURLToPath } from 'node:url'
14+import puppeteer from 'puppeteer'
15+
16+const root = join(dirname(fileURLToPath(import.meta.url)), '..')
17+const PORT = 5197
18+const URL = `http://localhost:${PORT}/`
19+const SHOTS = process.argv.includes('--shots')
20+const ONLY = process.argv.includes('--only') ? process.argv[process.argv.indexOf('--only') + 1] : null
21+
22+let failures = 0
23+const check = (name, ok, detail = '') => {
24+ console.log(`${ok ? ' ok ' : ' FAIL '} ${name}${detail ? ` ${detail}` : ''}`)
25+ if (!ok) failures++
26+}
27+
28+const server = spawn('npx', ['vite', 'preview', '--port', String(PORT), '--strictPort'], {
29+ cwd: root,
30+ stdio: ['ignore', 'pipe', 'pipe'],
31+ detached: true,
32+})
33+let stopped = false
34+const stop = () => {
35+ if (stopped) return
36+ stopped = true
37+ try {
38+ process.kill(-server.pid, 'SIGKILL')
39+ } catch {
40+ /* already gone */
41+ }
42+}
43+process.on('exit', stop)
44+process.on('SIGINT', () => {
45+ stop()
46+ process.exit(130)
47+})
48+
49+{
50+ const deadline = Date.now() + 30000
51+ let up = false
52+ while (Date.now() < deadline) {
53+ try {
54+ const res = await fetch(URL)
55+ if (res.ok) {
56+ up = true
57+ break
58+ }
59+ } catch {
60+ /* not listening yet */
61+ }
62+ await new Promise((r) => setTimeout(r, 250))
63+ }
64+ if (!up) {
65+ stop()
66+ throw new Error(`vite preview never answered on ${URL}`)
67+ }
68+}
69+
70+mkdirSync(join(root, 'tmp'), { recursive: true })
71+const browser = await puppeteer.launch({ headless: true, args: ['--no-sandbox'] })
72+const page = await browser.newPage()
73+await page.setViewport({ width: 1400, height: 1100 })
74+
75+const problems = []
76+// numbl announces which LAPACK bridge it picked on console.error; that is
77+// information, not a fault
78+const BENIGN = /using bridge:/
79+page.on('console', (m) => {
80+ if (m.type() === 'error' && !BENIGN.test(m.text())) problems.push(`console.error: ${m.text()}`)
81+})
82+page.on('pageerror', (e) => problems.push(`pageerror: ${e.message}`))
83+page.on('requestfailed', (r) => problems.push(`requestfailed: ${r.url()}`))
84+
85+const settle = (timeout = 240000) =>
86+ page.waitForFunction(
87+ () => {
88+ const labs = [...document.querySelectorAll('.lab')]
89+ if (labs.length === 0) return document.querySelector('.chapter') !== null
90+ return labs.every((l) => {
91+ const busy = [...l.querySelectorAll('button')].some((b) => /running/.test(b.textContent))
92+ return !busy && !l.querySelector('.lab-placeholder')
93+ })
94+ },
95+ { timeout, polling: 400 },
96+ )
97+
98+try {
99+ await page.goto(URL, { waitUntil: 'networkidle2' })
100+ const chapters = await page.evaluate(() =>
101+ [...document.querySelectorAll('.nav-item')].map((b) => b.textContent.replace(/^\d+/, '').trim()),
102+ )
103+ check('the sidebar lists chapters', chapters.length > 1, `${chapters.length} entries`)
104+
105+ for (let i = 0; i < chapters.length; i++) {
106+ const name = chapters[i]
107+ if (ONLY && !name.toLowerCase().includes(ONLY.toLowerCase())) continue
108+ console.log(`\n${name}`)
109+ await page.evaluate((k) => {
110+ document.querySelectorAll('.nav-item')[k].click()
111+ }, i)
112+ await new Promise((r) => setTimeout(r, 300))
113+ // labs that wait to be asked: ask them
114+ await page.evaluate(() => {
115+ for (const lab of document.querySelectorAll('.lab')) {
116+ if (!lab.querySelector('.lab-placeholder')) continue
117+ const b = [...lab.querySelectorAll('button')].find((x) => /^run$/.test(x.textContent.trim()))
118+ b?.click()
119+ }
120+ })
121+ await new Promise((r) => setTimeout(r, 300))
122+ try {
123+ await settle()
124+ } catch {
125+ check(`${name}: labs finish`, false, 'timed out')
126+ continue
127+ }
128+ const info = await page.evaluate(() => {
129+ const labs = [...document.querySelectorAll('.lab')]
130+ return {
131+ labs: labs.length,
132+ errors: [...document.querySelectorAll('.lab-error')].map((e) => e.textContent.slice(0, 300)),
133+ drawn: [...document.querySelectorAll('.lab-figures svg path')].filter(
134+ (p) => (p.getAttribute('d') ?? '').length > 12,
135+ ).length,
136+ canvases: document.querySelectorAll('.lab-figures canvas').length,
137+ stats: [...document.querySelectorAll('.stat-value')].map((s) => s.textContent),
138+ }
139+ })
140+ check(`${name}: no MATLAB error`, info.errors.length === 0, info.errors.join(' | '))
141+ if (info.labs > 0) {
142+ check(
143+ `${name}: figures carry data`,
144+ info.drawn + info.canvases > 0,
145+ `${info.drawn} paths, ${info.canvases} canvases`,
146+ )
147+ }
148+ if (info.labs > 0) console.log(` ${info.stats.join(' ')}`)
149+ if (SHOTS) {
150+ const file = join(root, 'tmp', `ch-${String(i + 1).padStart(2, '0')}.png`)
151+ await page.screenshot({ path: file, fullPage: true })
152+ console.log(` shot: ${file}`)
153+ }
154+ }
155+
156+ check('no console errors', problems.length === 0, problems.slice(0, 4).join(' | '))
157+} finally {
158+ await browser.close()
159+ stop()
160+}
161+
162+console.log(failures === 0 ? '\nall checks passed' : `\n${failures} check(s) failed`)
163+process.exit(failures === 0 ? 0 : 1)
scripts/matlab-test.mjsadded+235−0View file
@@ -0,0 +1,235 @@
1+#!/usr/bin/env node
2+// Checks the MATLAB layer against the numbers the chapters claim, through the
3+// numbl CLI, with no browser involved.
4+//
5+// npm run test:matlab (NUMBL_DIR defaults to ~/src/numbl)
6+//
7+// Every assertion here is a statement made in the prose, so if one of these
8+// fails either the code or the text is wrong.
9+import { join, dirname } from 'node:path'
10+import { fileURLToPath } from 'node:url'
11+import { runLab } from './runlab.mjs'
12+
13+const root = join(dirname(fileURLToPath(import.meta.url)), '..')
14+const src = (p) => join(root, 'src', p)
15+
16+let failures = 0
17+const check = (name, ok, detail = '') => {
18+ console.log(`${ok ? ' ok ' : ' FAIL '} ${name}${detail ? ` ${detail}` : ''}`)
19+ if (!ok) failures++
20+}
21+const near = (a, b, tol) => Math.abs(a - b) <= tol
22+const last = (row) => {
23+ for (let i = row.length - 1; i >= 0; i--) if (row[i] != null) return row[i]
24+ return null
25+}
26+const best = (row) => Math.min(...row.filter((v) => v != null))
27+
28+const ONE = src('chapters/onedim/driver.m')
29+const TWO = src('chapters/twodim/driver.m')
30+
31+console.log('\nChapter 1: three-point differences')
32+{
33+ const { out } = runLab(ONE, src('chapters/ch01/fd2.m'), {
34+ kind: 'smooth',
35+ mode: 'sweep',
36+ ns: [16, 32, 64, 128, 256, 512],
37+ nmaxs: [1e9],
38+ timed: 0,
39+ })
40+ const q = last(out.orders[0])
41+ check('observed order is 2', near(q, 2, 0.02), `q = ${q.toFixed(3)}`)
42+ const { out: low } = runLab(ONE, src('chapters/ch01/fd2.m'), {
43+ kind: 'lowreg',
44+ mode: 'sweep',
45+ ns: [64, 128, 256, 512],
46+ nmaxs: [1e9],
47+ timed: 0,
48+ })
49+ const ql = last(low.orders[0])
50+ check('the low-regularity solution converges more slowly', ql < 1.7, `q = ${ql.toFixed(3)}`)
51+}
52+
53+console.log('\nChapter 2: the roundoff floor and Richardson')
54+{
55+ const { out } = runLab(ONE, src('chapters/ch01/fd2.m'), {
56+ kind: 'smooth',
57+ mode: 'sweep',
58+ ns: [1024, 4096, 16384, 65536, 262144],
59+ nmaxs: [1e9],
60+ timed: 0,
61+ })
62+ const e = out.E[0]
63+ const bi = e.indexOf(best(e))
64+ check('the error stops improving before N = 262144', bi < e.length - 1, `best at N = ${out.ns[bi]}`)
65+ check('the best error is a few times 1e-9', best(e) < 3e-8 && best(e) > 1e-10, `${best(e).toExponential(2)}`)
66+
67+ const { out: r } = runLab(ONE, src('chapters/ch02/richardson.m'), {
68+ kind: 'smooth',
69+ mode: 'sweep',
70+ ns: [32, 64, 128, 256, 512],
71+ nmaxs: [1e9, 1e9],
72+ timed: 0,
73+ })
74+ check('Richardson gives order 4', near(last(r.orders[1]), 4, 0.1), `q = ${last(r.orders[1]).toFixed(3)}`)
75+ check('and it beats the plain scheme by >1e4', last(r.E[0]) / last(r.E[1]) > 1e4)
76+
77+ const { out: rl } = runLab(ONE, src('chapters/ch02/richardson.m'), {
78+ kind: 'lowreg',
79+ mode: 'sweep',
80+ ns: [64, 128, 256, 512],
81+ nmaxs: [1e9, 1e9],
82+ timed: 0,
83+ })
84+ check(
85+ 'on lowreg the extrapolation buys no extra order',
86+ near(last(rl.orders[1]), last(rl.orders[0]), 0.1) && last(rl.orders[1]) < 1.7,
87+ `q = ${last(rl.orders[1]).toFixed(2)} against ${last(rl.orders[0]).toFixed(2)}`,
88+ )
89+}
90+
91+console.log('\nChapter 3: the compact fourth-order scheme')
92+{
93+ const { out } = runLab(ONE, src('chapters/ch03/highorder.m'), {
94+ kind: 'smooth',
95+ mode: 'sweep',
96+ ns: [64, 256, 1024],
97+ nmaxs: [1e9, 1e9],
98+ timed: 0,
99+ })
100+ check('Numerov is fourth order', near(last(out.orders[1]), 4, 0.05), `q = ${last(out.orders[1]).toFixed(3)}`)
101+ check('at N = 1024 it is below 1e-9', last(out.E[1]) < 1e-9, last(out.E[1]).toExponential(2))
102+}
103+
104+console.log('\nChapter 4: Chebyshev collocation')
105+{
106+ const { out } = runLab(ONE, src('chapters/ch04/spectral.m'), {
107+ kind: 'smooth',
108+ mode: 'sweep',
109+ ns: [16, 32, 64, 128],
110+ nmaxs: [1e9, 1e9, 1e9],
111+ timed: 0,
112+ })
113+ check('spectral reaches 1e-13 by N = 64', out.E[2][2] < 1e-12, out.E[2][2].toExponential(2))
114+ const { out: low } = runLab(ONE, src('chapters/ch04/spectral.m'), {
115+ kind: 'lowreg',
116+ mode: 'sweep',
117+ ns: [64, 128, 256],
118+ nmaxs: [1e9, 1e9, 1e9],
119+ timed: 0,
120+ })
121+ check('on lowreg it is no better than the three-point scheme', low.E[2][2] > low.E[0][2])
122+}
123+
124+console.log('\nChapters 5 and 6: the five-point Laplacian, two solvers')
125+{
126+ const { out } = runLab(TWO, src('chapters/ch06/fastpoisson.m'), {
127+ kind: 'expsin',
128+ par: 1,
129+ mode: 'sweep',
130+ ns: [8, 16, 32, 128, 512],
131+ nmaxs: [32, 512],
132+ timed: 0,
133+ })
134+ check('the five-point scheme is second order', near(last(out.orders[1]), 2, 0.02), `q = ${last(out.orders[1]).toFixed(3)}`)
135+ for (let j = 0; j < 3; j++) {
136+ check(
137+ `dense and transform agree at N = ${out.ns[j]}`,
138+ Math.abs(out.E[0][j] - out.E[1][j]) < 1e-9 * out.E[0][j],
139+ `${out.E[0][j].toExponential(4)} vs ${out.E[1][j].toExponential(4)}`,
140+ )
141+ }
142+ check('at N = 512 the error is 4e-8', near(last(out.E[1]) / 3.99e-8, 1, 0.05), last(out.E[1]).toExponential(3))
143+}
144+
145+console.log('\nChapter 7: the compact nine-point stencil and 2-D Chebyshev')
146+{
147+ const { out } = runLab(TWO, src('chapters/ch07/highorder2d.m'), {
148+ kind: 'expsin',
149+ par: 1,
150+ mode: 'sweep',
151+ ns: [8, 16, 32, 64],
152+ nmaxs: [512, 512, 64],
153+ timed: 0,
154+ })
155+ check('nine-point is sixth order on Laplace', near(out.orders[1][1], 6, 0.1), `q = ${out.orders[1][1].toFixed(2)}`)
156+ check('nine-point reaches 1e-13 by N = 32', out.E[1][2] < 1e-13, out.E[1][2].toExponential(2))
157+ check('Chebyshev reaches 1e-13 by N = 16', out.E[2][1] < 1e-13, out.E[2][1].toExponential(2))
158+}
159+
160+console.log('\nChapter 8: iterative solvers')
161+{
162+ const { out } = runLab(src('chapters/ch08/driver.m'), src('chapters/ch08/iterative.m'), {
163+ kind: 'expsin',
164+ par: 1,
165+ N: 64,
166+ nit: 60,
167+ tol: 1e-11,
168+ })
169+ const mg = out.R2[4].filter((v) => v != null)
170+ check('multigrid converges in fewer than 15 cycles', mg.length < 16, `${mg.length - 1} cycles`)
171+ const mgErr = out.E2[4].filter((v) => v != null).map(Math.sqrt)
172+ check(
173+ 'and it lands on the discretisation error',
174+ near(mgErr[mgErr.length - 1] / out.discErr, 1, 0.02),
175+ `${mgErr[mgErr.length - 1].toExponential(3)} vs ${out.discErr.toExponential(3)}`,
176+ )
177+ const jac = out.R2[0].filter((v) => v != null).map(Math.sqrt)
178+ check('Jacobi is nowhere near after 60 sweeps', jac[jac.length - 1] > 1e-2, jac[jac.length - 1].toExponential(2))
179+}
180+
181+console.log('\nChapter 9: the L-shaped domain')
182+{
183+ const NSM = [
184+ [16, 32, 64],
185+ [8, 16, 32],
186+ ]
187+ const { out: b } = runLab(src('chapters/ch09/driver.m'), src('chapters/ch09/lshape.m'), {
188+ kind: 'both',
189+ mode: 'sweep',
190+ nsm: NSM,
191+ })
192+ check('the difference scheme loses second order', near(last(b.orders[0]), 1.37, 0.1), `q = ${last(b.orders[0]).toFixed(3)}`)
193+ check('the harmonic basis reaches 1e-11', last(b.E[1]) < 1e-11, last(b.E[1]).toExponential(2))
194+
195+ const { out: s } = runLab(src('chapters/ch09/driver.m'), src('chapters/ch09/lshape.m'), {
196+ kind: 'smooth',
197+ mode: 'sweep',
198+ nsm: NSM,
199+ })
200+ check('with smooth data the corner costs nothing', near(last(s.orders[0]), 2, 0.02), `q = ${last(s.orders[0]).toFixed(3)}`)
201+
202+ const { out: g } = runLab(src('chapters/ch09/driver.m'), src('chapters/ch09/lshape.m'), {
203+ kind: 'sing',
204+ mode: 'sweep',
205+ nsm: NSM,
206+ })
207+ check('the singular function is in the basis exactly', g.E[1][0] < 1e-14, g.E[1][0].toExponential(2))
208+}
209+
210+console.log('\nChapter 10: boundary integral equations')
211+{
212+ const { out } = runLab(src('chapters/ch10/driver.m'), src('chapters/ch10/bie.m'), {
213+ a: 0.3,
214+ nlobe: 5,
215+ par: 2.2,
216+ mode: 'sweep',
217+ ns: [64, 128, 256, 512],
218+ })
219+ check('Nyström reaches machine precision by n = 512', last(out.E[0]) < 1e-14, last(out.E[0]).toExponential(2))
220+ check('convergence is geometric, not algebraic', out.E[0][0] / out.E[0][1] > 100)
221+ check(
222+ 'close evaluation is orders of magnitude worse',
223+ last(out.Enear[0]) / last(out.E[0]) > 1e6,
224+ `${last(out.Enear[0]).toExponential(2)} at 98% radius`,
225+ )
226+ check('fundamental solutions reach 1e-13', best(out.E[1]) < 1e-13, best(out.E[1]).toExponential(2))
227+ check(
228+ 'and they have no close-evaluation problem',
229+ best(out.Enear[1]) < 1e-12,
230+ best(out.Enear[1]).toExponential(2),
231+ )
232+}
233+
234+console.log(failures === 0 ? '\nall checks passed' : `\n${failures} check(s) failed`)
235+process.exit(failures === 0 ? 0 : 1)
scripts/runlab.mjsadded+48−0View file
@@ -0,0 +1,48 @@
1+// Run one lab's MATLAB layer through the numbl CLI, outside the browser.
2+//
3+// node scripts/runlab.mjs <driver.m> <script.m> '<params json>' [--show key,key]
4+//
5+// It composes main.m exactly as the app does (driver, newline, script), copies
6+// the shared library into a scratch directory, runs it, and prints out.json.
7+// NUMBL_DIR points at a numbl checkout (default ~/src/numbl).
8+import { execFileSync } from 'node:child_process'
9+import { mkdtempSync, mkdirSync, readFileSync, writeFileSync, readdirSync, copyFileSync } from 'node:fs'
10+import { tmpdir } from 'node:os'
11+import { join, dirname, resolve } from 'node:path'
12+import { fileURLToPath } from 'node:url'
13+
14+const here = dirname(fileURLToPath(import.meta.url))
15+const root = resolve(here, '..')
16+const numbl = process.env.NUMBL_DIR ?? join(process.env.HOME, 'src', 'numbl')
17+const cli = join(numbl, 'dist-cli', 'cli.js')
18+
19+export function runLab(driverPath, scriptPath, params, opts = {}) {
20+ const dir = mkdtempSync(join(tmpdir(), 'lapl-'))
21+ const lib = join(root, 'src', 'matlab')
22+ for (const f of readdirSync(lib)) if (f.endsWith('.m')) copyFileSync(join(lib, f), join(dir, f))
23+ const driver = readFileSync(driverPath, 'utf8')
24+ const script = readFileSync(scriptPath, 'utf8')
25+ writeFileSync(join(dir, 'main.m'), `${driver}\n${script}`)
26+ writeFileSync(join(dir, 'params.json'), JSON.stringify(params))
27+ const t0 = Date.now()
28+ const stdout = execFileSync('node', [cli, 'run', 'main.m'], { cwd: dir, encoding: 'utf8' })
29+ const ms = Date.now() - t0
30+ const out = JSON.parse(readFileSync(join(dir, 'out.json'), 'utf8'))
31+ if (opts.verbose && stdout.trim()) process.stdout.write(stdout)
32+ return { out, stdout, ms, dir }
33+}
34+
35+if (process.argv[1] === fileURLToPath(import.meta.url)) {
36+ const [driverPath, scriptPath, paramsJson] = process.argv.slice(2)
37+ if (!driverPath || !scriptPath) {
38+ console.error('usage: node scripts/runlab.mjs <driver.m> <script.m> \'<params json>\'')
39+ process.exit(2)
40+ }
41+ const { out, stdout, ms } = runLab(driverPath, scriptPath, JSON.parse(paramsJson ?? '{}'))
42+ if (stdout.trim()) process.stdout.write(stdout)
43+ const show = process.argv.includes('--full')
44+ console.log(`[${ms} ms]`)
45+ console.log(JSON.stringify(out, null, show ? 1 : 0).slice(0, show ? 1e9 : 4000))
46+}
47+
48+mkdirSync(join(root, 'tmp'), { recursive: true })
src/App.tsxadded+92−0View file
@@ -0,0 +1,92 @@
1+import { useCallback, useEffect, useState } from 'react'
2+import ChapterView from './ui/ChapterView.tsx'
3+import LibraryView from './ui/LibraryView.tsx'
4+import { CHAPTERS } from './chapters/index.ts'
5+
6+const PARTS = ['One dimension', 'Two dimensions', 'Beyond the square'] as const
7+
8+function readHash(): string {
9+ const h = window.location.hash.replace(/^#\/?/, '')
10+ return h === '' ? CHAPTERS[0].id : h
11+}
12+
13+export default function App() {
14+ const [route, setRoute] = useState(readHash)
15+ useEffect(() => {
16+ const on = () => setRoute(readHash())
17+ window.addEventListener('hashchange', on)
18+ return () => window.removeEventListener('hashchange', on)
19+ }, [])
20+
21+ const go = useCallback((id: string) => {
22+ window.location.hash = `/${id}`
23+ window.scrollTo({ top: 0 })
24+ }, [])
25+
26+ const chapter = CHAPTERS.find((c) => c.id === route)
27+ const index = CHAPTERS.findIndex((c) => c.id === route)
28+
29+ return (
30+ <div className="app">
31+ <nav className="sidebar">
32+ <div className="brand">
33+ <div className="brand-title">The Dirichlet problem</div>
34+ <div className="brand-sub">accuracy against cost, in the browser</div>
35+ </div>
36+ {PARTS.map((part) => (
37+ <div className="nav-part" key={part}>
38+ <div className="nav-part-title">{part}</div>
39+ <ol>
40+ {CHAPTERS.filter((c) => c.part === part).map((c) => (
41+ <li key={c.id}>
42+ <button className={c.id === route ? 'nav-item active' : 'nav-item'} onClick={() => go(c.id)}>
43+ <span className="nav-num">{CHAPTERS.indexOf(c) + 1}</span>
44+ {c.nav}
45+ </button>
46+ </li>
47+ ))}
48+ </ol>
49+ </div>
50+ ))}
51+ <div className="nav-part">
52+ <div className="nav-part-title">Reference</div>
53+ <ol>
54+ <li>
55+ <button className={route === 'library' ? 'nav-item active' : 'nav-item'} onClick={() => go('library')}>
56+ <span className="nav-num" />
57+ shared MATLAB library
58+ </button>
59+ </li>
60+ </ol>
61+ </div>
62+ <div className="nav-foot">
63+ MATLAB syntax runs in your browser through <a href="https://numbl.org">numbl</a>. Nothing is installed and
64+ nothing is sent anywhere.
65+ </div>
66+ </nav>
67+
68+ <main className="main">
69+ {route === 'library' ? (
70+ <LibraryView />
71+ ) : chapter ? (
72+ <>
73+ <ChapterView chapter={chapter} />
74+ <nav className="pager">
75+ {index > 0 && (
76+ <button onClick={() => go(CHAPTERS[index - 1].id)}>
77+ ← {CHAPTERS[index - 1].nav}
78+ </button>
79+ )}
80+ <span className="lab-bar-gap" />
81+ {index < CHAPTERS.length - 1 && (
82+ <button onClick={() => go(CHAPTERS[index + 1].id)}>{CHAPTERS[index + 1].nav} →</button>
83+ )}
84+ </nav>
85+ </>
86+ ) : (
87+ <p className="lab-error">no such chapter: {route}</p>
88+ )}
89+ </main>
90+ </div>
91+ )
92+}
src/chapters/ch01/chapter.tsadded+65−0View file
@@ -0,0 +1,65 @@
1+import { createElement } from 'react'
2+import type { Chapter } from '../../ui/types.ts'
3+import driver from '../onedim/driver.m?raw'
4+import script from './fd2.m?raw'
5+import text from './text.md?raw'
6+import { Convergence, Solution1D, fmtErr, type SingleOut, type SweepOut } from '../onedim/panels.tsx'
7+import { KIND_CONTROL, NS_1D } from '../shared.ts'
8+
9+export const ch01: Chapter = {
10+ id: 'manufactured',
11+ nav: 'Ground truth by construction',
12+ part: 'One dimension',
13+ title: 'Ground truth by construction',
14+ blurb:
15+ 'The one-dimensional model problem, the method of manufactured solutions, and the three-point difference scheme.',
16+ text,
17+ labs: [
18+ {
19+ id: 'first',
20+ title: 'One solve',
21+ caption: 'Move the resolution and change the solution being manufactured.',
22+ driver,
23+ script,
24+ controls: [
25+ KIND_CONTROL,
26+ {
27+ kind: 'slider',
28+ key: 'N',
29+ label: 'intervals N',
30+ min: 4,
31+ max: 128,
32+ values: [4, 8, 12, 16, 24, 32, 48, 64, 96, 128, 256],
33+ },
34+ ],
35+ defaults: { kind: 'smooth', mode: 'single', mid: 1, N: 16, ngrid: 700 },
36+ editorOpen: true,
37+ render: (d: SingleOut) => createElement(Solution1D, { data: d }),
38+ summary: (d: SingleOut) => [
39+ { label: 'unknowns', value: String(d.n - 1) },
40+ { label: 'max error', value: fmtErr(d.maxerr) },
41+ ],
42+ },
43+ {
44+ id: 'refine',
45+ title: 'Refinement study',
46+ caption: 'The same script, run over a sequence of meshes, with the observed order beside each error.',
47+ driver,
48+ script,
49+ controls: [KIND_CONTROL],
50+ defaults: { kind: 'smooth', mode: 'sweep', ns: NS_1D, nmaxs: [1e9, 1e9], timed: 0 },
51+ render: (d: SweepOut) =>
52+ createElement(Convergence, {
53+ data: d,
54+ style: { slots: ['fd2'], refs: [{ p: -2, label: 'slope -2' }] },
55+ }),
56+ summary: (d: SweepOut) => {
57+ const q = d.orders[0][d.orders[0].length - 1]
58+ return [
59+ { label: 'finest error', value: fmtErr(d.E[0][d.E[0].length - 1]) },
60+ { label: 'observed order', value: q == null ? '--' : q.toFixed(2) },
61+ ]
62+ },
63+ },
64+ ],
65+}
src/chapters/ch01/fd2.madded+47−0View file
@@ -0,0 +1,47 @@
1+% Second-order finite differences for -u'' = f on (0, 1), u(0) = ua, u(1) = ub.
2+%
3+% Replace u''(x_i) by the three-point second difference
4+%
5+% u''(x_i) ~ ( u(i-1) - 2 u(i) + u(i+1) ) / h^2 ,
6+%
7+% which is exact for quadratics and wrong by (h^2/12) u''''(x_i) for anything
8+% else. Collecting the N-1 interior equations gives a tridiagonal system whose
9+% matrix is (1/h^2) * tridiag(-1, 2, -1); the two known boundary values move to
10+% the right-hand side.
11+
12+function names = method_names()
13+names = {'3-point differences'};
14+end
15+
16+function [x, u] = solve1d(mid, N, f, ua, ub)
17+h = 1 / N;
18+x = (0:N)' / N;
19+
20+r = h^2 * f(x(2:N)); % scale by h^2 so the matrix is tridiag(-1, 2, -1)
21+r(1) = r(1) + ua; % the neighbour at x = 0 is known
22+r(end) = r(end) + ub; % and so is the one at x = 1
23+
24+u = [ua; thomas(N-1, r); ub];
25+end
26+
27+function w = thomas(n, d)
28+% Solve tridiag(-1, 2, -1) w = d by Gaussian elimination without pivoting,
29+% which is safe here because the matrix is symmetric positive definite and
30+% diagonally dominant. Two sweeps, O(n) work and O(n) storage: the whole
31+% reason a one-dimensional problem is cheap at any resolution.
32+c = zeros(n, 1);
33+e = zeros(n, 1);
34+b = 2;
35+c(1) = -1 / b;
36+e(1) = d(1) / b;
37+for i = 2:n
38+ b = 2 + c(i-1); % = 2 - (-1)*c(i-1)
39+ c(i) = -1 / b;
40+ e(i) = (d(i) + e(i-1)) / b;
41+end
42+w = zeros(n, 1);
43+w(n) = e(n);
44+for i = n-1:-1:1
45+ w(i) = e(i) - c(i) * w(i+1);
46+end
47+end
src/chapters/ch01/text.mdadded+51−0View file
@@ -0,0 +1,51 @@
1+The Dirichlet problem asks for a function that satisfies a differential equation inside a region and takes prescribed values on its boundary. For Laplace's equation the interior condition is $\Delta u = 0$, and in one dimension that reads $u'' = 0$, whose solution is the straight line through the two boundary values. There is nothing left to compute. The one-dimensional model problem worth solving numerically is therefore Poisson's equation,
2+
3+$$-u''(x) = f(x), \qquad 0 < x < 1, \qquad u(0) = \alpha, \quad u(1) = \beta.$$
4+
5+It has the same second-difference structure, the same conditioning, and the same convergence theory as the two-dimensional problem, at a hundredth of the cost, so it is the right place to develop the tools. From Chapter 5 on we set $f = 0$ and let the second dimension supply the difficulty.
6+
7+### Knowing the answer in advance
8+
9+A convergence claim needs an error, and an error needs a truth to measure against. Reference solutions computed by the same family of methods at higher resolution are the usual substitute, and they are unreliable in exactly the situations one wants to study: near a singularity, at high frequency, or wherever the method being tested is systematically wrong.
10+
11+Here we use manufactured solutions throughout. Choose $u$ first, differentiate it twice by hand, and *define* the data by
12+
13+$$f := -u'', \qquad \alpha := u(0), \qquad \beta := u(1).$$
14+
15+Then $u$ is the exact solution of the resulting boundary value problem, by construction. The error of any computed approximation is known at every point, to full precision, at every resolution, with no reference computation and no extrapolation.
16+
17+The obvious limitation is that we are testing a method against a solution we chose. A method can be flattered by a convenient choice, so the five manufactured solutions in these labs differ deliberately in how hard they are to resolve: an entire function of moderate amplitude, one with eight wavelengths across the interval, a narrow Gaussian, a boundary layer of width $1/50$, and one that is $C^2$ but not $C^3$. Anything that is claimed below should be checked against all five, and the last two are where the claims get interesting.
18+
19+### The three-point scheme
20+
21+Put $N$ equal intervals of length $h = 1/N$ on $[0,1]$, write $u_i$ for the value at $x_i = ih$, and replace the second derivative by the second difference,
22+
23+$$u''(x_i) \;\approx\; \frac{u_{i-1} - 2u_i + u_{i+1}}{h^2}.$$
24+
25+A Taylor expansion of the right-hand side gives $u''(x_i) + \frac{h^2}{12}u''''(x_i) + O(h^4)$, so the substitution is exact for quadratics and wrong by a term proportional to $h^2$ for anything else. Imposing the resulting equation at the $N-1$ interior nodes, and moving the two known boundary values to the right-hand side, gives
26+
27+$$\frac{1}{h^2}\begin{pmatrix} 2 & -1 & & \\ -1 & 2 & \ddots & \\ & \ddots & \ddots & -1 \\ & & -1 & 2\end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \\ \vdots \\ u_{N-1}\end{pmatrix} = \begin{pmatrix} f_1 + \alpha/h^2 \\ f_2 \\ \vdots \\ f_{N-1} + \beta/h^2 \end{pmatrix}.$$
28+
29+The matrix is symmetric positive definite and tridiagonal, so Gaussian elimination without pivoting needs one forward and one backward sweep, $O(N)$ operations and $O(N)$ storage. Nothing about one dimension is expensive.
30+
31+<!--lab:first-->
32+
33+Two things are worth noticing before moving on. The error is not noise: it is a smooth function of $x$, of one sign over long stretches, with a shape that does not change as $N$ grows. Only its amplitude changes. And the amplitude is a strong function of which solution you ask for. At $N = 32$ the smooth case is already accurate to a few parts in $10^4$, while the boundary layer is not resolved at all, because a mesh that never places a point inside the layer cannot know it is there.
34+
35+### Why the error behaves that way
36+
37+Let $e_i = u_i - u(x_i)$ be the nodal error. Applying the discrete operator to $e$ and using the Taylor expansion above,
38+
39+$$\frac{e_{i-1} - 2e_i + e_{i+1}}{h^2} \;=\; \frac{h^2}{12}\,u''''(x_i) + O(h^4),$$
40+
41+with $e_0 = e_N = 0$. The error therefore solves the *same* discrete problem, with the truncation error as its source. Since the discrete solution operator is bounded independently of $h$, the error inherits the size of that source:
42+
43+$$e \;\approx\; \frac{h^2}{12}\, w, \qquad -w'' = u'''', \quad w(0) = w(1) = 0.$$
44+
45+That is the whole of the classical theory in one line, and it predicts two things the lab above shows. The error is a fixed smooth function, $w/12$, scaled by $h^2$, which is why its shape is stable under refinement. And its size is set by $u''''$, which is why the narrow Gaussian and the boundary layer are hard: their fourth derivatives are enormous.
46+
47+<!--lab:refine-->
48+
49+The observed order settles at exactly 2 for the first four solutions. For `lowreg`, where $u''''$ does not exist at $x = 1/2$, the argument above does not apply, and the observed order is lower. This is the first appearance of a theme that runs through the rest of the tutorial: convergence rates are statements about the smoothness of the solution, not about the method alone.
50+
51+Second order is a modest rate. Each extra digit of accuracy costs a factor of $\sqrt{10} \approx 3.2$ in $N$, so three more digits costs a factor of 32 in work. The next chapter asks how far that can be pushed before something other than the discretization gets in the way.
src/chapters/ch02/chapter.tsadded+96−0View file
@@ -0,0 +1,96 @@
1+import { createElement } from 'react'
2+import type { Chapter } from '../../ui/types.ts'
3+import driver from '../onedim/driver.m?raw'
4+import fd2 from '../ch01/fd2.m?raw'
5+import richardson from './richardson.m?raw'
6+import text from './text.md?raw'
7+import { Convergence, fmtErr, type SweepOut } from '../onedim/panels.tsx'
8+import { KIND_CONTROL, NS_1D } from '../shared.ts'
9+
10+export const ch02: Chapter = {
11+ id: 'limits',
12+ nav: 'Second order to its limit',
13+ part: 'One dimension',
14+ title: 'Second order, to its limit',
15+ blurb: 'Conditioning, the roundoff floor at about eight digits, and Richardson extrapolation.',
16+ text,
17+ labs: [
18+ {
19+ id: 'floor',
20+ title: 'Refining until it stops helping',
21+ caption:
22+ 'Nine resolutions from 16 to 262144 intervals. The last few take a second or two each; the whole sweep runs on request.',
23+ driver,
24+ script: fd2,
25+ manual: true,
26+ controls: [
27+ KIND_CONTROL,
28+ {
29+ kind: 'slider',
30+ key: 'nmax',
31+ label: 'finest N',
32+ values: [16384, 65536, 262144, 1048576],
33+ min: 0,
34+ max: 3,
35+ format: (v) => String(v),
36+ },
37+ ],
38+ defaults: { kind: 'smooth', mode: 'sweep', ns: [], nmax: 262144, nmaxs: [1e9], timed: 1 },
39+ derive: (p) => ({
40+ ...p,
41+ ns: [16, 64, 256, 1024, 4096, 16384, 65536, 262144, 1048576].filter((n) => n <= Number(p.nmax)),
42+ }),
43+ render: (d: SweepOut) =>
44+ createElement(Convergence, {
45+ data: d,
46+ style: {
47+ slots: ['fd2'],
48+ refs: [
49+ { p: -2, label: 'slope −2: discretisation' },
50+ { p: 2, label: 'slope +2: roundoff × κ(A)', at: [1e3, 2e-13] },
51+ ],
52+ },
53+ height: 340,
54+ }),
55+ summary: (d: SweepOut) => {
56+ const es = d.E[0]
57+ let bi = 0
58+ es.forEach((v, i) => {
59+ if (v != null && (es[bi] == null || v < (es[bi] as number))) bi = i
60+ })
61+ return [
62+ { label: 'best error', value: fmtErr(es[bi]) },
63+ { label: 'at N', value: String(d.ns[bi]) },
64+ ]
65+ },
66+ },
67+ {
68+ id: 'richardson',
69+ title: 'Two second-order solves, one fourth-order answer',
70+ caption: 'Change the manufactured solution to lowreg and watch the gain evaporate.',
71+ driver,
72+ script: richardson,
73+ editorOpen: true,
74+ controls: [KIND_CONTROL],
75+ defaults: { kind: 'smooth', mode: 'sweep', ns: NS_1D, nmaxs: [1e9, 1e9], timed: 0 },
76+ render: (d: SweepOut) =>
77+ createElement(Convergence, {
78+ data: d,
79+ style: {
80+ slots: ['fd2', 'fd4'],
81+ refs: [
82+ { p: -2, label: 'slope −2' },
83+ { p: -4, label: 'slope −4' },
84+ ],
85+ },
86+ }),
87+ summary: (d: SweepOut) => {
88+ const q = d.orders[1]?.[d.orders[1].length - 1]
89+ return [
90+ { label: 'extrapolated error', value: fmtErr(d.E[1]?.[d.E[1].length - 1]) },
91+ { label: 'observed order', value: q == null ? '--' : q.toFixed(2) },
92+ ]
93+ },
94+ },
95+ ],
96+}
src/chapters/ch02/richardson.madded+56−0View file
@@ -0,0 +1,56 @@
1+% Two methods on the same tridiagonal solver:
2+%
3+% 1. the three-point scheme of the previous chapter;
4+% 2. the same scheme run at N and at 2N, combined by Richardson
5+% extrapolation to cancel the leading h^2 term of the error.
6+%
7+% If the error has an asymptotic expansion u_h = u + c h^2 + O(h^4), then
8+%
9+% (4 u_{h/2} - u_h) / 3 = u + O(h^4)
10+%
11+% at the points the two meshes share. Nothing about the discretisation
12+% changes; the fourth-order accuracy comes out of two second-order solves and
13+% one subtraction. The expansion is what makes it work, and the expansion
14+% needs u to be smooth: on the low-regularity solution there is no h^2 term to
15+% cancel and the extrapolation gains nothing.
16+
17+function names = method_names()
18+names = {'3-point differences', 'Richardson from N and 2N'};
19+end
20+
21+function [x, u] = solve1d(mid, N, f, ua, ub)
22+if mid == 1
23+ [x, u] = fd2(N, f, ua, ub);
24+else
25+ [x, uc] = fd2(N, f, ua, ub);
26+ [~, uf] = fd2(2*N, f, ua, ub);
27+ u = (4 * uf(1:2:end) - uc) / 3;
28+end
29+end
30+
31+function [x, u] = fd2(N, f, ua, ub)
32+h = 1 / N;
33+x = (0:N)' / N;
34+r = h^2 * f(x(2:N));
35+r(1) = r(1) + ua;
36+r(end) = r(end) + ub;
37+u = [ua; thomas(N-1, r); ub];
38+end
39+
40+function w = thomas(n, d)
41+c = zeros(n, 1);
42+e = zeros(n, 1);
43+b = 2;
44+c(1) = -1 / b;
45+e(1) = d(1) / b;
46+for i = 2:n
47+ b = 2 + c(i-1);
48+ c(i) = -1 / b;
49+ e(i) = (d(i) + e(i-1)) / b;
50+end
51+w = zeros(n, 1);
52+w(n) = e(n);
53+for i = n-1:-1:1
54+ w(i) = e(i) - c(i) * w(i+1);
55+end
56+end
src/chapters/ch02/text.mdadded+41−0View file
@@ -0,0 +1,41 @@
1+Second order says that halving $h$ divides the error by four. Taken literally it promises any accuracy at all, given enough points. It does not deliver, and the two reasons it does not are the subject of this chapter: arithmetic is not exact, and the matrix gets worse as the mesh gets finer.
2+
3+### The condition number of the difference matrix
4+
5+The eigenvalues of the $(N-1) \times (N-1)$ matrix $\mathrm{tridiag}(-1,2,-1)$ are known in closed form,
6+
7+$$\lambda_k \;=\; 4\sin^2\!\Big(\frac{k\pi h}{2}\Big), \qquad k = 1, \dots, N-1,$$
8+
9+so after the factor $1/h^2$ the extreme eigenvalues of $A$ are about $\pi^2$ and $4/h^2$, and
10+
11+$$\kappa(A) \;=\; \frac{\lambda_{N-1}}{\lambda_1} \;\approx\; \frac{4}{\pi^2 h^2} \;\approx\; 0.4\,N^2 .$$
12+
13+This is not a defect of the discretisation. The continuous operator $-d^2/dx^2$ is unbounded, and any consistent discretisation of it must have a condition number that grows as the mesh is refined. The practical consequence is a floor: solving in double precision, with unit roundoff $\varepsilon \approx 2.2 \times 10^{-16}$, cannot be expected to give a relative accuracy better than roughly $\varepsilon \kappa(A) \approx 10^{-16} N^2$. Setting that equal to the discretisation error $C h^2 = C N^{-2}$ gives
14+
15+$$N_{\text{best}} \;\sim\; \varepsilon^{-1/4} \approx 10^4, \qquad \text{error}_{\text{best}} \;\sim\; \sqrt{\varepsilon} \approx 10^{-8}.$$
16+
17+About half the digits, and no more, whatever the machine and however long one is willing to wait.
18+
19+<!--lab:floor-->
20+
21+The prediction is quantitative and the lab confirms it. The error falls at exactly second order until $N \approx 10^5$, reaches a little under $10^{-8}$, and then climbs again. The climb is roundoff amplified by the condition number, and its detailed rate depends on how the rounding errors happen to fall, but the shape of the curve is completely reproducible: past the minimum, refining the mesh makes the answer worse and takes longer to do it.
22+
23+Note also the cost. Each row of the table is a solve of $O(N)$ work, so the wall clock in the last column grows linearly, and the whole 20-fold improvement in accuracy from $N = 10^3$ to $N = 10^5$ costs a factor of 100 in time.
24+
25+### Getting a better rate without a better scheme
26+
27+The error expansion of the previous chapter says more than "the error is $O(h^2)$". For smooth $u$ the nodal error has an asymptotic expansion in even powers of $h$,
28+
29+$$u_h(x) = u(x) + c_2(x)\,h^2 + c_4(x)\,h^4 + \cdots,$$
30+
31+with coefficients independent of $h$. Two solves, at $h$ and $h/2$, then give two equations at the points the meshes share, and the combination
32+
33+$$u_R \;=\; \frac{4\,u_{h/2} - u_h}{3}$$
34+
35+cancels $c_2$ exactly and leaves $O(h^4)$. This is Richardson extrapolation. It costs one extra solve, which for a tridiagonal system is a factor of three in total work, and it needs no new discretisation, no new stencil, and no new code.
36+
37+<!--lab:richardson-->
38+
39+On the smooth solution the observed order goes from 2 to 4, and at $N = 512$ the extrapolated answer is about five orders of magnitude better than the one it was built from. Switch to `lowreg` and the gain in *order* disappears entirely: both curves converge at about $1.5$, because $u''''$ does not exist at $x = 1/2$ and the expansion in powers of $h^2$ that justifies the combination is not valid there. The extrapolation still improves the constant, by a factor of about seven, which is a fair return for one extra solve but is not what it was asked to do.
40+
41+That is the honest summary of extrapolation. It is nearly free when it works, it is easy to check whether it worked, since the observed order tells you, and when the assumed expansion is absent it degenerates into a mildly better constant. The next chapter takes the other route: change the discretisation itself.
src/chapters/ch03/chapter.tsadded+52−0View file
@@ -0,0 +1,52 @@
1+import { createElement, Fragment } from 'react'
2+import type { Chapter } from '../../ui/types.ts'
3+import driver from '../onedim/driver.m?raw'
4+import script from './highorder.m?raw'
5+import text from './text.md?raw'
6+import { Convergence, WorkPrecision, fmtErr, type SweepOut } from '../onedim/panels.tsx'
7+import { KIND_CONTROL } from '../shared.ts'
8+
9+const NS = [16, 64, 256, 1024, 4096, 16384]
10+
11+export const ch03: Chapter = {
12+ id: 'higher-order',
13+ nav: 'Two more orders, almost free',
14+ part: 'One dimension',
15+ title: 'Two more orders, almost free',
16+ blurb: 'The compact fourth-order scheme, and how to read a work-precision diagram.',
17+ text,
18+ labs: [
19+ {
20+ id: 'orders',
21+ title: 'Second order against fourth',
22+ caption:
23+ 'Both schemes solve the same tridiagonal system; only the right-hand side differs. The second figure is the one that decides anything.',
24+ driver,
25+ script,
26+ editorOpen: true,
27+ controls: [KIND_CONTROL],
28+ defaults: { kind: 'smooth', mode: 'sweep', ns: NS, nmaxs: [16384, 16384], timed: 1 },
29+ render: (d: SweepOut) =>
30+ createElement(
31+ Fragment,
32+ null,
33+ createElement(Convergence, {
34+ key: 'c',
35+ data: d,
36+ style: {
37+ slots: ['fd2', 'fd4'],
38+ refs: [
39+ { p: -2, label: 'slope −2' },
40+ { p: -4, label: 'slope −4' },
41+ ],
42+ },
43+ }),
44+ createElement(WorkPrecision, { key: 'w', data: d, style: { slots: ['fd2', 'fd4'] } }),
45+ ),
46+ summary: (d: SweepOut) => [
47+ { label: 'best 2nd order', value: fmtErr(Math.min(...(d.E[0].filter((v) => v != null) as number[]))) },
48+ { label: 'best 4th order', value: fmtErr(Math.min(...(d.E[1].filter((v) => v != null) as number[]))) },
49+ ],
50+ },
51+ ],
52+}
src/chapters/ch03/highorder.madded+55−0View file
@@ -0,0 +1,55 @@
1+% Two discretisations that share a matrix.
2+%
3+% 1. the three-point scheme: (u_{i-1} - 2u_i + u_{i+1}) / h^2 = -f_i
4+% 2. the compact fourth-order scheme, due to Numerov:
5+%
6+% (u_{i-1} - 2u_i + u_{i+1}) / h^2 = -(f_{i-1} + 10 f_i + f_{i+1}) / 12
7+%
8+% The left-hand sides are identical. All that changes is that the source is
9+% averaged over the same three points with weights 1/12, 10/12, 1/12, which is
10+% exactly the combination that cancels the h^2 term of the truncation error.
11+% The matrix is still tridiag(-1, 2, -1), the solve is still O(N), and the
12+% method is fourth-order accurate: the extra two orders cost one weighted
13+% average of numbers we already had.
14+%
15+% The catch is that the scheme wants f at the nodes and nothing else, so it
16+% needs f to be smooth as well: the cancellation is a statement about the
17+% Taylor expansion of f, not just of u.
18+
19+function names = method_names()
20+names = {'3-point differences', 'compact 4th order (Numerov)'};
21+end
22+
23+function [x, u] = solve1d(mid, N, f, ua, ub)
24+h = 1 / N;
25+x = (0:N)' / N;
26+fx = f(x);
27+
28+if mid == 1
29+ r = h^2 * fx(2:N);
30+else
31+ r = h^2 * (fx(1:N-1) + 10*fx(2:N) + fx(3:N+1)) / 12;
32+end
33+r(1) = r(1) + ua;
34+r(end) = r(end) + ub;
35+
36+u = [ua; thomas(N-1, r); ub];
37+end
38+
39+function w = thomas(n, d)
40+c = zeros(n, 1);
41+e = zeros(n, 1);
42+b = 2;
43+c(1) = -1 / b;
44+e(1) = d(1) / b;
45+for i = 2:n
46+ b = 2 + c(i-1);
47+ c(i) = -1 / b;
48+ e(i) = (d(i) + e(i-1)) / b;
49+end
50+w = zeros(n, 1);
51+w(n) = e(n);
52+for i = n-1:-1:1
53+ w(i) = e(i) - c(i) * w(i+1);
54+end
55+end
src/chapters/ch03/text.mdadded+29−0View file
@@ -0,0 +1,29 @@
1+Richardson extrapolation improved the answer without changing the discretisation. The alternative is to change the discretisation, and in one dimension there is a version of that which is almost free.
2+
3+### A fourth-order scheme with the same matrix
4+
5+Write $\delta^2 u_i = u_{i-1} - 2u_i + u_{i+1}$. The Taylor expansion of the second difference is
6+
7+$$\frac{\delta^2 u_i}{h^2} \;=\; u''_i + \frac{h^2}{12}u''''_i + \frac{h^4}{360}u''''''_i + \cdots,$$
8+
9+and since $u'' = -f$, the offending term is $-\frac{h^2}{12} f''_i$. Approximating that with the second difference of $f$ itself, which is all that is needed at this order, gives
10+
11+$$\frac{\delta^2 u_i}{h^2} \;=\; -f_i - \frac{h^2}{12}\,\frac{\delta^2 f_i}{h^2} + O(h^4) \;=\; -\frac{f_{i-1} + 10 f_i + f_{i+1}}{12} + O(h^4).$$
12+
13+The scheme this defines is Numerov's, and the important feature is on the left. The matrix is unchanged: still $\mathrm{tridiag}(-1,2,-1)$, still symmetric positive definite, still solved in $O(N)$ by two sweeps. The right-hand side is a weighted average of values of $f$ we already had. Two extra orders of accuracy cost one pass over an array.
14+
15+Schemes with this property, where a higher order is bought with a wider stencil on the *data* but not on the *unknowns*, are called compact. They matter more in two dimensions than in one, and Chapter 7 uses the two-dimensional version.
16+
17+<!--lab:orders-->
18+
19+### Reading a work-precision diagram
20+
21+The convergence plot answers "how does the error behave as $N$ grows", which is a question about the method. The question a user actually has is "what does it cost to get eight digits", and the answer to that is the second figure: the same errors, plotted against the measured time for one solve.
22+
23+The times are measured inside the MATLAB layer, around the solve and nothing else, with a warm-up call first and the best of several repeats for anything under five milliseconds. Below about a tenth of a millisecond what is being measured is the clock and the fixed overhead of a call, not the method, which is why the curves flatten at the left. Any absolute number quoted in the text is from the command-line runtime on the author's machine; in a browser the same code is a few times slower, and your own numbers are in the strip above each figure.
24+
25+Read the diagram by drawing a horizontal line at the accuracy you need and taking the leftmost curve. On the smooth solution, at $10^{-7}$, the compact scheme is faster than the three-point scheme by more than an order of magnitude, and the gap widens as the accuracy demand tightens, because the two curves have different slopes. That is the general shape of these comparisons: a higher order does not make a method faster at low accuracy, it makes the *slope* steeper, and the slope is what decides who wins once the accuracy requirement is serious.
26+
27+Two honest qualifications. First, the fourth-order curve stops improving near $5 \times 10^{-12}$, for the reason established in the last chapter: the condition number of the same matrix is still growing like $N^2$, and roundoff catches up with a fourth-order discretisation earlier, in terms of $N$, than with a second-order one. Higher order reaches the floor sooner, and it reaches it at a much better error.
28+
29+Second, all of this is contingent on smoothness. Run the sweep on `lowreg` and the compact scheme converges at very nearly the same rate as the three-point scheme, at a slightly better constant. Where the function has no fourth derivative, no amount of Taylor-series bookkeeping produces fourth-order accuracy. The next chapter pushes the smooth case to its limit, and then asks what happens to that limit when smoothness is withdrawn.
src/chapters/ch04/chapter.tsadded+83−0View file
@@ -0,0 +1,83 @@
1+import { createElement, Fragment } from 'react'
2+import type { Chapter } from '../../ui/types.ts'
3+import driver from '../onedim/driver.m?raw'
4+import script from './spectral.m?raw'
5+import text from './text.md?raw'
6+import { Convergence, Solution1D, WorkPrecision, fmtErr, type SingleOut, type SweepOut } from '../onedim/panels.tsx'
7+import { KIND_CONTROL } from '../shared.ts'
8+import type { MethodKey } from '../../plot/palette.ts'
9+
10+const NS = [8, 16, 32, 64, 128, 256, 1024, 4096, 16384]
11+const SLOTS: MethodKey[] = ['fd2', 'fd4', 'spectral']
12+
13+export const ch04: Chapter = {
14+ id: 'spectral',
15+ nav: 'Chebyshev collocation',
16+ part: 'One dimension',
17+ title: 'Chebyshev collocation',
18+ blurb: 'A dense matrix, geometric convergence, and the conditions under which that is worth it.',
19+ text,
20+ labs: [
21+ {
22+ id: 'single',
23+ title: 'One spectral solve',
24+ caption:
25+ 'Note where the nodes are. The error is plotted at the nodes, so at N = 32 there are 33 of them and the answer is already near machine precision.',
26+ driver,
27+ script,
28+ controls: [
29+ KIND_CONTROL,
30+ {
31+ kind: 'select',
32+ key: 'mid',
33+ label: 'method',
34+ options: [
35+ { value: 3, label: 'Chebyshev collocation' },
36+ { value: 2, label: 'compact 4th order' },
37+ { value: 1, label: '3-point differences' },
38+ ],
39+ },
40+ { kind: 'slider', key: 'N', label: 'N', values: [4, 8, 12, 16, 24, 32, 48, 64, 96, 128], min: 4, max: 128 },
41+ ],
42+ defaults: { kind: 'smooth', mode: 'single', mid: 3, N: 24, ngrid: 700 },
43+ editorOpen: true,
44+ render: (d: SingleOut, p) =>
45+ createElement(Solution1D, { data: d, slot: SLOTS[Number(p.mid) - 1] ?? 'spectral' }),
46+ summary: (d: SingleOut) => [
47+ { label: 'unknowns', value: String(d.n - 1) },
48+ { label: 'max error', value: fmtErr(d.maxerr) },
49+ ],
50+ },
51+ {
52+ id: 'sweep',
53+ title: 'All three, by resolution and by cost',
54+ caption: 'Chebyshev is capped at N = 256, where the dense solve starts to dominate the timing.',
55+ driver,
56+ script,
57+ controls: [KIND_CONTROL],
58+ defaults: { kind: 'smooth', mode: 'sweep', ns: NS, nmaxs: [16384, 4096, 256], timed: 1 },
59+ render: (d: SweepOut) =>
60+ createElement(
61+ Fragment,
62+ null,
63+ createElement(Convergence, {
64+ key: 'c',
65+ data: d,
66+ style: {
67+ slots: SLOTS,
68+ refs: [
69+ { p: -2, label: 'slope −2' },
70+ { p: -4, label: 'slope −4' },
71+ ],
72+ },
73+ height: 340,
74+ }),
75+ createElement(WorkPrecision, { key: 'w', data: d, style: { slots: SLOTS }, height: 340 }),
76+ ),
77+ summary: (d: SweepOut) => [
78+ { label: 'best spectral', value: fmtErr(Math.min(...(d.E[2].filter((v) => v != null) as number[]))) },
79+ { label: 'best 3-point', value: fmtErr(Math.min(...(d.E[0].filter((v) => v != null) as number[]))) },
80+ ],
81+ },
82+ ],
83+}
src/chapters/ch04/spectral.madded+73−0View file
@@ -0,0 +1,73 @@
1+% Three methods on the same problem: the two finite difference schemes of the
2+% previous chapters, and Chebyshev collocation.
3+%
4+% Collocation drops the idea of a local stencil. Put the N+1 Chebyshev points
5+%
6+% x_j = ( 1 - cos(pi j / N) ) / 2, j = 0, ..., N,
7+%
8+% on [0,1], let p be the polynomial of degree N through the unknown values, and
9+% require -p'' = f at the interior points. Differentiation of a polynomial
10+% given by its values at fixed nodes is a matrix, so the discrete operator is
11+% the dense matrix D^2 rather than a three-term recurrence. Every unknown is
12+% coupled to every other, the solve costs O(N^3) instead of O(N), and in
13+% exchange the error decays faster than any power of N for analytic solutions.
14+%
15+% The clustering of the nodes near the ends is not decoration. It is what makes
16+% polynomial interpolation of high degree stable, and equally spaced points of
17+% the same degree would be useless.
18+
19+function names = method_names()
20+names = {'3-point differences', 'compact 4th order', 'Chebyshev collocation'};
21+end
22+
23+function [x, u] = solve1d(mid, N, f, ua, ub)
24+if mid == 3
25+ [x, u] = cheb_solve(N, f, ua, ub);
26+else
27+ [x, u] = fd_solve(mid, N, f, ua, ub);
28+end
29+end
30+
31+function [x, u] = cheb_solve(N, f, ua, ub)
32+[D, s] = lp_cheb(N);
33+x = (1 - s) / 2; % Chebyshev points, ascending, on [0, 1]
34+Dx = -2 * D; % d/dx after the affine map
35+L = -(Dx * Dx); % the operator -d^2/dx^2, as a dense matrix
36+
37+ii = 2:N; % interior collocation points
38+A = L(ii, ii);
39+b = f(x(ii)) - L(ii, 1) * ua - L(ii, N+1) * ub;
40+u = [ua; A \ b; ub];
41+end
42+
43+function [x, u] = fd_solve(mid, N, f, ua, ub)
44+h = 1 / N;
45+x = (0:N)' / N;
46+fx = f(x);
47+if mid == 1
48+ r = h^2 * fx(2:N);
49+else
50+ r = h^2 * (fx(1:N-1) + 10*fx(2:N) + fx(3:N+1)) / 12;
51+end
52+r(1) = r(1) + ua;
53+r(end) = r(end) + ub;
54+u = [ua; thomas(N-1, r); ub];
55+end
56+
57+function w = thomas(n, d)
58+c = zeros(n, 1);
59+e = zeros(n, 1);
60+b = 2;
61+c(1) = -1 / b;
62+e(1) = d(1) / b;
63+for i = 2:n
64+ b = 2 + c(i-1);
65+ c(i) = -1 / b;
66+ e(i) = (d(i) + e(i-1)) / b;
67+end
68+w = zeros(n, 1);
69+w(n) = e(n);
70+for i = n-1:-1:1
71+ w(i) = e(i) - c(i) * w(i+1);
72+end
73+end
src/chapters/ch04/text.mdadded+33−0View file
@@ -0,0 +1,33 @@
1+A difference scheme of order $p$ has an error $C h^p$: on a log-log plot, a straight line of slope $-p$. Raising $p$ tilts the line. The question of this chapter is what happens if $p$ is not fixed at all, but grows with the mesh.
2+
3+### Collocation with a global polynomial
4+
5+Let $p$ be the polynomial of degree $N$ that takes the unknown values $u_0, \dots, u_N$ at the $N+1$ Chebyshev points
6+
7+$$x_j \;=\; \tfrac{1}{2}\Big(1 - \cos\frac{\pi j}{N}\Big), \qquad j = 0, \dots, N,$$
8+
9+and require $-p''(x_j) = f(x_j)$ at the $N-1$ interior points, with $p$ matching the boundary data at the two ends. Differentiation is linear and the nodes are fixed, so $p'$ at the nodes is a matrix times the vector of nodal values. That matrix $D$ is the Chebyshev differentiation matrix, `lp_cheb` in the shared library; the discrete operator is $-D^2$, dense and of size $N+1$.
10+
11+Two things are being given up. The matrix is dense, so a solve costs $O(N^3)$ rather than $O(N)$. And the accuracy at a point now depends on the data everywhere, so there is no notion of a local truncation error, and the classical Taylor-series analysis does not apply.
12+
13+What is bought is the approximation power of polynomials of high degree. For $u$ analytic in a neighbourhood of $[0,1]$, the best degree-$N$ polynomial approximation converges geometrically, like $\rho^{-N}$ for some $\rho > 1$ set by the size of that neighbourhood, and the collocation solution inherits that rate. The error is not $C N^{-p}$ for any $p$; on a log-log plot it is not a straight line at all, but a curve that falls away ever more steeply.
14+
15+The clustering of the nodes near the endpoints is essential rather than cosmetic. Equally spaced points of the same degree produce the Runge phenomenon, and a collocation method built on them is unusable for $N$ beyond about 20.
16+
17+<!--lab:single-->
18+
19+### Where the curve stops
20+
21+<!--lab:sweep-->
22+
23+On the smooth solution the Chebyshev error falls to $10^{-13}$ at $N = 64$, in a fraction of a millisecond. The three-point scheme, given a hundred times as long and $16384$ intervals, manages $10^{-7}$. This is the pattern that makes spectral methods worth their inconvenience: when the solution is analytic and the geometry is simple, they are not a little better, they are several orders of magnitude better at the same cost.
24+
25+The flat tail past $N = 64$ is the roundoff floor again, and it arrives for the same reason as before. The two-derivative operator $D^2$ has a condition number growing like $N^4$, worse than the $N^2$ of the difference matrices, so a spectral method reaches its floor at a rather small $N$. It reaches it at about $10^{-13}$, which is close to what double precision permits, and there is nothing to be gained by going further.
26+
27+Now change the manufactured solution.
28+
29+- `layer` puts a boundary layer of width $1/50$ at $x = 1$, which is exactly where the Chebyshev points cluster, and the method converges even faster than it does on the smooth case. Grid clustering that happens to match the solution's structure is worth a great deal.
30+- `peak` puts the structure in the middle, where the nodes are sparsest. The method still converges geometrically, but only once $N$ is large enough to resolve the peak, and the curve stays flat until then.
31+- `lowreg` is the interesting failure. The solution is $C^2$ but not $C^3$, and all three methods collapse onto very nearly the same algebraic rate, close to $N^{-3/2}$. Chebyshev collocation is *no better than the three-point scheme* here, and at moderate $N$ it is somewhat worse. Geometric convergence is a statement about analyticity, and when analyticity fails there is nothing left of the advantage.
32+
33+The general rule is worth stating plainly, because it survives into two dimensions and into the harder geometries of Part III: the achievable convergence rate is a property of the solution, and a method can only fail to attain it. High-order methods are how one attains it when the solution is smooth. Chapter 9, where the geometry itself forces a singularity into the solution, is about what to do when it is not.
src/chapters/ch05/chapter.tsadded+72−0View file
@@ -0,0 +1,72 @@
1+import { createElement } from 'react'
2+import type { Chapter } from '../../ui/types.ts'
3+import driver from '../twodim/driver.m?raw'
4+import script from './fivepoint.m?raw'
5+import text from './text.md?raw'
6+import { Fields2D, type Single2DOut } from '../twodim/panels.tsx'
7+import { Convergence, fmtErr, type SweepOut } from '../onedim/panels.tsx'
8+import { EXACT2D_CONTROL } from '../shared.ts'
9+import type { Control } from '../../ui/types.ts'
10+
11+const PAR_CONTROL: Control = {
12+ kind: 'slider',
13+ key: 'par',
14+ label: 'shape parameter',
15+ min: 1,
16+ max: 8,
17+ step: 1,
18+ format: (v) => String(v),
19+}
20+
21+export const ch05: Chapter = {
22+ id: 'five-point',
23+ nav: 'The five-point Laplacian',
24+ part: 'Two dimensions',
25+ title: 'The five-point Laplacian',
26+ blurb: 'Harmonic functions as ground truth, the standard stencil, and the cost of solving the system as it comes.',
27+ text,
28+ labs: [
29+ {
30+ id: 'solve',
31+ title: 'One solve, and the error field',
32+ caption:
33+ 'The shape parameter is the degree k for poly, the wavenumber k for wave, and the distance from the square for log.',
34+ driver,
35+ script,
36+ editorOpen: true,
37+ controls: [
38+ EXACT2D_CONTROL,
39+ PAR_CONTROL,
40+ { kind: 'slider', key: 'N', label: 'N', values: [8, 12, 16, 24, 32, 40], min: 8, max: 40 },
41+ ],
42+ defaults: { kind: 'expsin', par: 1, mode: 'single', mid: 1, N: 24, maxpix: 97 },
43+ render: (d: Single2DOut) => createElement(Fields2D, { data: d }),
44+ summary: (d: Single2DOut) => [
45+ { label: 'unknowns', value: String(d.unknowns) },
46+ { label: 'max error', value: fmtErr(d.maxerr) },
47+ ],
48+ },
49+ {
50+ id: 'cost',
51+ title: 'Second order, at O(N⁶)',
52+ caption: 'Five resolutions. The last one is a 2209 by 2209 dense solve and takes a few seconds in a browser.',
53+ driver,
54+ script,
55+ manual: true,
56+ controls: [EXACT2D_CONTROL, PAR_CONTROL],
57+ defaults: { kind: 'expsin', par: 1, mode: 'sweep', ns: [4, 8, 16, 32, 48], nmaxs: [48], timed: 1 },
58+ render: (d: SweepOut) =>
59+ createElement(Convergence, {
60+ data: d,
61+ style: { slots: ['fd2'], refs: [{ p: -2, label: 'slope −2' }] },
62+ }),
63+ summary: (d: SweepOut) => [
64+ { label: 'finest error', value: fmtErr(d.E[0][d.E[0].length - 1]) },
65+ {
66+ label: 'time for that',
67+ value: `${((d.T[0][d.T[0].length - 1] ?? 0) / 1000).toFixed(2)} s`,
68+ },
69+ ],
70+ },
71+ ],
72+}
src/chapters/ch05/fivepoint.madded+52−0View file
@@ -0,0 +1,52 @@
1+% The five-point Laplacian on the unit square, assembled as a matrix and
2+% solved by Gaussian elimination.
3+%
4+% On a uniform grid of N intervals in each direction, replace
5+%
6+% u_xx + u_yy by ( u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) - 4 u(i,j) ) / h^2
7+%
8+% at the (N-1)^2 interior points. Scaling by h^2 leaves the matrix
9+%
10+% A = kron(I, T) + kron(T, I), T = tridiag(-1, 2, -1),
11+%
12+% which is what the vec of the interior array is multiplied by. Boundary
13+% neighbours are known, so they move to the right-hand side.
14+%
15+% This is the honest, naive version: the matrix is formed in full and handed to
16+% backslash. It is also the version that stops working first. The unknown count
17+% is (N-1)^2, so the dense matrix has (N-1)^4 entries and the elimination costs
18+% (N-1)^6 / 3 operations: N = 32 needs 7 MB and a fraction of a second, N = 64
19+% needs 450 MB, and N = 128 is out of reach on any machine.
20+
21+function names = method_names()
22+names = {'5-point, dense solve'};
23+end
24+
25+function [x, y, U] = solve2d(mid, N, g)
26+h = 1 / N;
27+x = (0:N)' / N;
28+y = x;
29+n = N - 1;
30+
31+% boundary values, and zeros in the interior
32+U = zeros(N+1, N+1);
33+U(1, :) = g(x', y(1));
34+U(N+1, :) = g(x', y(N+1));
35+U(:, 1) = g(x(1), y);
36+U(:, N+1) = g(x(N+1), y);
37+
38+% right-hand side: h^2 f (zero here, since we are solving Laplace's equation)
39+% plus the known neighbours of the points next to the boundary
40+B = zeros(n, n);
41+B(1, :) = B(1, :) + U(1, 2:N);
42+B(n, :) = B(n, :) + U(N+1, 2:N);
43+B(:, 1) = B(:, 1) + U(2:N, 1);
44+B(:, n) = B(:, n) + U(2:N, N+1);
45+
46+T = 2*eye(n) - diag(ones(n-1, 1), 1) - diag(ones(n-1, 1), -1);
47+I = eye(n);
48+A = kron(I, T) + kron(T, I);
49+
50+W = A \ B(:);
51+U(2:N, 2:N) = reshape(W, n, n);
52+end
src/chapters/ch05/text.mdadded+44−0View file
@@ -0,0 +1,44 @@
1+From here on the equation is Laplace's,
2+
3+$$\Delta u = u_{xx} + u_{yy} = 0 \ \ \text{in } \Omega, \qquad u = g \ \ \text{on } \partial\Omega,$$
4+
5+and $\Omega$ is the unit square. Poisson's equation is no harder to discretise, and everything below applies to it unchanged, but with $f = 0$ the only error in the computation is the error of the method, which is what we are trying to measure.
6+
7+### Ground truth in two dimensions
8+
9+The manufactured-solution trick needs an adjustment. Choosing $u$ freely and setting $f := -\Delta u$ would give a Poisson problem, not a Laplace problem. Instead we choose $u$ from the supply of functions that are *already* harmonic, and take its boundary trace as the data. Four are available in these labs:
10+
11+| choice | $u$ | why it is here |
12+|---|---|---|
13+| `expsin` | $e^{x}\sin y$ | entire, gentle, the standard test |
14+| `poly` | $\operatorname{Re}(x+iy)^k$ | a harmonic polynomial; exact for a scheme of order $k$ or higher |
15+| `wave` | $e^{ky}\sin(kx)$ | oscillates along one edge and grows by $e^{k}$ across the square |
16+| `log` | $\log|z - z_0|$, $z_0$ just outside | harmonic inside, but with a singularity a distance $d$ from the boundary |
17+
18+The real and imaginary parts of any analytic function are harmonic, which is where all four come from, and it is worth noticing how much freedom that gives: any solution one can write down in closed form is available as ground truth, including ones that are nearly singular.
19+
20+### The five-point stencil
21+
22+Put a uniform grid of spacing $h = 1/N$ on the square and replace both second derivatives by their three-point differences. The result is the five-point formula
23+
24+$$\frac{u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1} - 4u_{i,j}}{h^2} \;=\; 0,$$
25+
26+one equation for each of the $(N-1)^2$ interior points, with the boundary values moved to the right-hand side. Its truncation error is $\frac{h^2}{12}(u_{xxxx} + u_{yyyy}) + O(h^4)$, so the scheme is second-order accurate, and the argument of Chapter 1 carries over unchanged: the error satisfies the same discrete equation with the truncation error as its source, and therefore inherits its size.
27+
28+Written in terms of the array of interior unknowns, the matrix has the Kronecker form $A = I \otimes T + T \otimes I$ with $T = \mathrm{tridiag}(-1,2,-1)$. That structure will matter a great deal in the next chapter. In this one we ignore it and hand the assembled matrix to Gaussian elimination.
29+
30+<!--lab:solve-->
31+
32+The error field is worth looking at rather than summarising. It is smooth, largest in the interior, and zero on the boundary, because the boundary values are exact by construction. Its shape is again the solution of a Poisson problem with the truncation error as source. Try `poly` with $k = 3$: the fourth derivatives of a cubic vanish, the leading truncation error is identically zero, and the scheme is exact to roundoff. Then try `log` with the singularity close to the edge, where $u_{xxxx}$ is large near that edge and the error concentrates there.
33+
34+### The wall
35+
36+<!--lab:cost-->
37+
38+The convergence is second order, cleanly. The cost is the problem. Dense elimination on $(N-1)^2$ unknowns costs $O(N^6)$ operations and $O(N^4)$ storage, and the measured times follow that faithfully: each doubling of $N$ multiplies the time by about 64. The last row of the table is a 2209 by 2209 dense solve, which buys five and a half correct digits and takes a fraction of a second at the command line, several seconds in a browser. $N = 128$ would need about 30 GB and is not a question of patience.
39+
40+Two remarks on what this does and does not show.
41+
42+The matrix is extremely sparse, with at most five entries per row, and treating it as dense is a waste. A real sparse direct solver with a good ordering, nested dissection, would solve this system in $O(N^3)$ operations with $O(N^2 \log N)$ fill, which is a completely different curve. numbl's backslash converts a sparse matrix to dense before solving, so that option is not available in this page, and the naive version is what the lab measures.
43+
44+But even nested dissection is not the right answer here. The matrix is not merely sparse: it is the Kronecker sum of two identical tridiagonal matrices, whose eigenvectors are known in closed form. The next chapter uses that to solve the same system in $O(N^2 \log N)$ operations, with the same answer to the last bit.
src/chapters/ch06/chapter.tsadded+75−0View file
@@ -0,0 +1,75 @@
1+import { createElement, Fragment } from 'react'
2+import type { Chapter } from '../../ui/types.ts'
3+import driver from '../twodim/driver.m?raw'
4+import script from './fastpoisson.m?raw'
5+import text from './text.md?raw'
6+import { Convergence, WorkPrecision, fmtErr, type SweepOut } from '../onedim/panels.tsx'
7+import { EXACT2D_CONTROL } from '../shared.ts'
8+import type { MethodKey } from '../../plot/palette.ts'
9+
10+const SLOTS: MethodKey[] = ['fd2', 'fast']
11+
12+export const ch06: Chapter = {
13+ id: 'fast-poisson',
14+ nav: 'The fast Poisson solver',
15+ part: 'Two dimensions',
16+ title: 'The fast Poisson solver',
17+ blurb: 'The same discretisation, diagonalised by the sine transform: a million unknowns in half a second.',
18+ text,
19+ labs: [
20+ {
21+ id: 'compare',
22+ title: 'Elimination against the sine transform',
23+ caption:
24+ 'Both curves are the same discretisation, so where they overlap they lie exactly on top of each other; the sine transform is drawn dashed so that the dense solve shows through. The dense solve stops at N = 32 because that is where it stops being possible.',
25+ driver,
26+ script,
27+ manual: true,
28+ editorOpen: true,
29+ controls: [
30+ EXACT2D_CONTROL,
31+ {
32+ kind: 'slider',
33+ key: 'nmax',
34+ label: 'finest N for the fast solver',
35+ values: [128, 256, 512, 1024, 2048],
36+ min: 0,
37+ max: 4,
38+ format: (v) => String(v),
39+ },
40+ ],
41+ defaults: {
42+ kind: 'expsin',
43+ par: 1,
44+ mode: 'sweep',
45+ ns: [],
46+ nmax: 512,
47+ nmaxs: [32, 512],
48+ timed: 1,
49+ },
50+ derive: (p) => ({
51+ ...p,
52+ ns: [8, 16, 32, 64, 128, 256, 512, 1024, 2048].filter((n) => n <= Number(p.nmax)),
53+ nmaxs: [32, Number(p.nmax)],
54+ }),
55+ render: (d: SweepOut) =>
56+ createElement(
57+ Fragment,
58+ null,
59+ createElement(Convergence, {
60+ key: 'c',
61+ data: d,
62+ style: { slots: SLOTS, dashes: [undefined, '7 5'], refs: [{ p: -2, label: 'slope −2' }] },
63+ }),
64+ createElement(WorkPrecision, { key: 'w', data: d, style: { slots: SLOTS, dashes: [undefined, '7 5'] } }),
65+ ),
66+ summary: (d: SweepOut) => {
67+ const last = d.E[1].length - 1
68+ return [
69+ { label: 'finest error', value: fmtErr(d.E[1][last]) },
70+ { label: 'in', value: `${(d.T[1][last] ?? 0).toFixed(0)} ms` },
71+ ]
72+ },
73+ },
74+ ],
75+}
src/chapters/ch06/fastpoisson.madded+52−0View file
@@ -0,0 +1,52 @@
1+% The same five-point system, solved two ways.
2+%
3+% The matrix is A = kron(I,T) + kron(T,I) with T = tridiag(-1,2,-1), and the
4+% eigenvectors of T are known exactly: the sine vectors
5+%
6+% v_k(j) = sin(pi j k / N), T v_k = lambda_k v_k,
7+% lambda_k = 2 (1 - cos(pi k / N)), k = 1, ..., N-1.
8+%
9+% So the eigenvectors of A are the products v_k(i) v_l(j) and its eigenvalues
10+% are the sums lambda_k + lambda_l. Expanding the right-hand side in that basis,
11+% dividing by the eigenvalues and transforming back solves the system exactly,
12+% and the expansion is a discrete sine transform, which one FFT computes in
13+% O(N log N) per line. Total: O(N^2 log N) for N^2 unknowns, no matrix formed
14+% and nothing stored beyond the grid itself.
15+%
16+% This is not an iterative method and not an approximation. It returns the same
17+% vector the dense solve returns, to within roundoff, and it does so for grids
18+% the dense solve cannot represent.
19+
20+function names = method_names()
21+names = {'5-point, dense solve', '5-point, sine transform'};
22+end
23+
24+function [x, y, U] = solve2d(mid, N, g)
25+h = 1 / N;
26+x = (0:N)' / N;
27+y = x;
28+n = N - 1;
29+
30+U = zeros(N+1, N+1);
31+U(1, :) = g(x', y(1));
32+U(N+1, :) = g(x', y(N+1));
33+U(:, 1) = g(x(1), y);
34+U(:, N+1) = g(x(N+1), y);
35+
36+B = zeros(n, n);
37+B(1, :) = B(1, :) + U(1, 2:N);
38+B(n, :) = B(n, :) + U(N+1, 2:N);
39+B(:, 1) = B(:, 1) + U(2:N, 1);
40+B(:, n) = B(:, n) + U(2:N, N+1);
41+
42+if mid == 1
43+ T = 2*eye(n) - diag(ones(n-1, 1), 1) - diag(ones(n-1, 1), -1);
44+ A = kron(eye(n), T) + kron(T, eye(n));
45+ U(2:N, 2:N) = reshape(A \ B(:), n, n);
46+else
47+ lam = 2 * (1 - cos(pi * (1:n)' / N)); % eigenvalues of T
48+ D = lam * ones(1, n) + ones(n, 1) * lam'; % eigenvalues of A
49+ Bh = lp_dst1(lp_dst1(B).').'; % to the sine basis, both ways
50+ U(2:N, 2:N) = lp_dst1(lp_dst1(Bh ./ D).').'; % divide, and back again
51+end
52+end
src/chapters/ch06/text.mdadded+36−0View file
@@ -0,0 +1,36 @@
1+The five-point matrix is not a general sparse matrix that happens to have five diagonals. It is
2+
3+$$A \;=\; I \otimes T + T \otimes I, \qquad T = \mathrm{tridiag}(-1,2,-1),$$
4+
5+and everything about $T$ is known in closed form. Its eigenvectors are the discrete sine vectors
6+
7+$$v_k(j) = \sin\frac{\pi j k}{N}, \qquad T v_k = \lambda_k v_k, \qquad \lambda_k = 2\Big(1 - \cos\frac{\pi k}{N}\Big) = 4\sin^2\frac{\pi k}{2N},$$
8+
9+for $k = 1, \dots, N-1$. A Kronecker sum has the products of the eigenvectors as its eigenvectors and the sums of the eigenvalues as its eigenvalues, so
10+
11+$$A\,\big(v_k \otimes v_l\big) \;=\; (\lambda_k + \lambda_l)\,\big(v_k \otimes v_l\big).$$
12+
13+The system is therefore diagonal in a basis we can write down. Expand the right-hand side in that basis, divide each coefficient by $\lambda_k + \lambda_l$, and transform back. The expansion is a two-dimensional discrete sine transform, and each one-dimensional sine transform is one FFT of the odd extension, so the whole solve costs $O(N^2 \log N)$ for $N^2$ unknowns, with no matrix formed and no storage beyond the grid.
14+
15+This is the fast Poisson solver, and the important thing about it is that it is *not* an approximation. It returns the solution of the same linear system the dense solve returns, to within roundoff, in the same sense that a better route to the same place is not a different place.
16+
17+<!--lab:compare-->
18+
19+Where the two methods overlap, the errors agree to every digit shown, which is the check that matters: the fast solver has not changed the discretisation, only the arithmetic used to invert it. Where they do not overlap, the dense solve is simply unavailable, and the sine transform is solving a million unknowns in about half a second.
20+
21+### What it buys, and what it costs
22+
23+The accuracy at fixed $N$ is unchanged, because the discretisation is unchanged. What changes is which $N$ is affordable, and with second-order accuracy that is worth about four digits: $N = 32$ gives $10^{-5}$, $N = 1024$ gives $10^{-8}$.
24+
25+It also shows the limit of the whole approach. Second order means the error is $C N^{-2}$, so twelve digits would need $N \approx 10^6$, which is $10^{12}$ unknowns. No solver, however fast, rescues a second-order discretisation at that accuracy. If we want twelve digits we have to stop refining and change the stencil, which is Chapter 7.
26+
27+### The conditions
28+
29+The derivation used every one of the following, and each is a real restriction.
30+
31+- **A rectangle**, so that the grid is a tensor product and the matrix is a Kronecker sum.
32+- **A uniform grid**, so that $T$ is the same in every line and its eigenvectors are sine vectors.
33+- **Constant coefficients**, so that the operator separates. A variable coefficient $\nabla\cdot(a\nabla u)$ destroys the structure at once.
34+- **Dirichlet conditions on all four sides.** Neumann or periodic conditions are equally fine, with the cosine or Fourier transform in place of the sine transform, but a mixture along one side is not.
35+
36+Break any of them and this chapter's method disappears, while the previous chapter's method merely gets slower. That trade, between methods that exploit structure and methods that tolerate its absence, is the theme of the rest of the tutorial. Chapter 8 asks what to do when the structure is gone and the matrix is too large for elimination.
src/chapters/ch07/chapter.tsadded+91−0View file
@@ -0,0 +1,91 @@
1+import { createElement, Fragment } from 'react'
2+import type { Chapter } from '../../ui/types.ts'
3+import driver from '../twodim/driver.m?raw'
4+import script from './highorder2d.m?raw'
5+import text from './text.md?raw'
6+import { Fields2D, type Single2DOut } from '../twodim/panels.tsx'
7+import { Convergence, WorkPrecision, fmtErr, type SweepOut } from '../onedim/panels.tsx'
8+import { EXACT2D_CONTROL } from '../shared.ts'
9+import type { MethodKey } from '../../plot/palette.ts'
10+
11+const SLOTS: MethodKey[] = ['fd2', 'fd4', 'spectral']
12+
13+export const ch07: Chapter = {
14+ id: 'high-order-2d',
15+ nav: 'Sixth order for the same price',
16+ part: 'Two dimensions',
17+ title: 'Sixth order for the same price',
18+ blurb: 'The compact nine-point stencil, Chebyshev collocation on the square, and a million-fold difference in unknowns.',
19+ text,
20+ labs: [
21+ {
22+ id: 'fields',
23+ title: 'The same grid, three stencils',
24+ caption: 'Watch the error scale on the right, not the picture on the left: at N = 32 they differ by ten orders of magnitude.',
25+ driver,
26+ script,
27+ controls: [
28+ EXACT2D_CONTROL,
29+ {
30+ kind: 'select',
31+ key: 'mid',
32+ label: 'method',
33+ options: [
34+ { value: 1, label: '5-point' },
35+ { value: 2, label: '9-point compact' },
36+ { value: 3, label: 'Chebyshev' },
37+ ],
38+ },
39+ { kind: 'slider', key: 'N', label: 'N', values: [8, 16, 24, 32, 48, 64], min: 8, max: 64 },
40+ { kind: 'slider', key: 'par', label: 'shape parameter', min: 1, max: 8, step: 1, format: (v) => String(v) },
41+ ],
42+ defaults: { kind: 'expsin', par: 1, mode: 'single', mid: 2, N: 32, maxpix: 97 },
43+ render: (d: Single2DOut) => createElement(Fields2D, { data: d }),
44+ summary: (d: Single2DOut) => [
45+ { label: 'unknowns', value: String(d.unknowns) },
46+ { label: 'max error', value: fmtErr(d.maxerr) },
47+ ],
48+ },
49+ {
50+ id: 'sweep',
51+ title: 'Order six against order two',
52+ caption: 'Chebyshev is capped at N = 64, which is already 3969 unknowns and machine precision.',
53+ driver,
54+ script,
55+ manual: true,
56+ editorOpen: true,
57+ controls: [EXACT2D_CONTROL, { kind: 'slider', key: 'par', label: 'shape parameter', min: 1, max: 8, step: 1, format: (v) => String(v) }],
58+ defaults: {
59+ kind: 'expsin',
60+ par: 1,
61+ mode: 'sweep',
62+ ns: [8, 16, 32, 64, 128, 256, 512],
63+ nmaxs: [512, 512, 64],
64+ timed: 1,
65+ },
66+ render: (d: SweepOut) =>
67+ createElement(
68+ Fragment,
69+ null,
70+ createElement(Convergence, {
71+ key: 'c',
72+ data: d,
73+ style: {
74+ slots: SLOTS,
75+ refs: [
76+ { p: -2, label: 'slope −2' },
77+ { p: -6, label: 'slope −6' },
78+ ],
79+ },
80+ height: 340,
81+ }),
82+ createElement(WorkPrecision, { key: 'w', data: d, style: { slots: SLOTS }, height: 340 }),
83+ ),
84+ summary: (d: SweepOut) => [
85+ { label: 'best 5-point', value: fmtErr(Math.min(...(d.E[0].filter((v) => v != null) as number[]))) },
86+ { label: 'best 9-point', value: fmtErr(Math.min(...(d.E[1].filter((v) => v != null) as number[]))) },
87+ { label: 'best Chebyshev', value: fmtErr(Math.min(...(d.E[2].filter((v) => v != null) as number[]))) },
88+ ],
89+ },
90+ ],
91+}
src/chapters/ch07/highorder2d.madded+98−0View file
@@ -0,0 +1,98 @@
1+% Three discretisations of the same Dirichlet problem.
2+%
3+% 1. the five-point stencil, solved by the sine transform (Chapter 6);
4+% 2. the compact nine-point stencil, solved by the same transform;
5+% 3. Chebyshev collocation in both directions, solved by diagonalising the
6+% one-dimensional operator.
7+%
8+% The nine-point stencil weights the four edge neighbours by 4 and the four
9+% corner neighbours by 1,
10+%
11+% ( 4(u_E + u_W + u_N + u_S) + (u_NE + u_NW + u_SE + u_SW) - 20 u_P ) / (6h^2),
12+%
13+% which is the Mehrstellen operator. Its truncation error is
14+% (h^2/12) Delta^2 u + O(h^4), and Delta^2 u = 0 for a harmonic u, so on
15+% Laplace's equation the leading two error terms both vanish and the scheme is
16+% sixth-order accurate. On Poisson's equation the same stencil is fourth-order
17+% unless the source is also averaged, which is the two-dimensional version of
18+% the trick in Chapter 3.
19+%
20+% It is still diagonalised by the sine basis, because it is still a sum of
21+% tensor products of tridiagonal operators: with a = pi k / N and b = pi l / N,
22+% the eigenvalue is 20 - 8 cos a - 8 cos b - 4 cos a cos b. So the sixth-order
23+% scheme costs exactly what the second-order one costs.
24+
25+function names = method_names()
26+names = {'5-point (sine transform)', '9-point compact (sine transform)', 'Chebyshev collocation'};
27+end
28+
29+function [x, y, U] = solve2d(mid, N, g)
30+if mid == 3
31+ [x, y, U] = cheb2d(N, g);
32+else
33+ [x, y, U] = fd2d(mid, N, g);
34+end
35+end
36+
37+function [x, y, U] = fd2d(mid, N, g)
38+x = (0:N)' / N;
39+y = x;
40+n = N - 1;
41+ii = 2:N;
42+
43+G = zeros(N+1, N+1); % boundary values, zero inside
44+G(1, :) = g(x', y(1));
45+G(N+1, :) = g(x', y(N+1));
46+G(:, 1) = g(x(1), y);
47+G(:, N+1) = g(x(N+1), y);
48+
49+a = pi * (1:n)' / N;
50+ca = cos(a);
51+
52+if mid == 1
53+ R = G(ii-1, ii) + G(ii+1, ii) + G(ii, ii-1) + G(ii, ii+1);
54+ M = 4 - 2*ca*ones(1, n) - 2*ones(n, 1)*ca';
55+else
56+ R = 4 * (G(ii-1, ii) + G(ii+1, ii) + G(ii, ii-1) + G(ii, ii+1)) ...
57+ + (G(ii-1, ii-1) + G(ii-1, ii+1) + G(ii+1, ii-1) + G(ii+1, ii+1));
58+ M = 20 - 8*ca*ones(1, n) - 8*ones(n, 1)*ca' - 4*(ca*ca');
59+end
60+
61+Rh = lp_dst1(lp_dst1(R).').';
62+U = G;
63+U(ii, ii) = lp_dst1(lp_dst1(Rh ./ M).').';
64+end
65+
66+function [x, y, U] = cheb2d(N, g)
67+% Collocation on the tensor product of Chebyshev grids. The discrete problem is
68+% the Sylvester equation A*U + U*A' = F, where A is the one-dimensional
69+% operator on interior points. Diagonalising A once, A = V diag(lam) inv(V),
70+% turns it into a division: the two-dimensional analogue of Chapter 6, with the
71+% eigenvectors computed rather than known.
72+[D, s] = lp_cheb(N);
73+x = (1 - s) / 2;
74+y = x;
75+Dx = -2 * D;
76+L = -(Dx * Dx);
77+n = N - 1;
78+ii = 2:N;
79+
80+G = zeros(N+1, N+1);
81+G(1, :) = g(x', y(1));
82+G(N+1, :) = g(x', y(N+1));
83+G(:, 1) = g(x(1), y);
84+G(:, N+1) = g(x(N+1), y);
85+
86+A = L(ii, ii);
87+F = -(L(ii, 1) * G(1, ii) + L(ii, N+1) * G(N+1, ii)) ...
88+ - (G(ii, 1) * L(ii, 1)' + G(ii, N+1) * L(ii, N+1)');
89+
90+[V, Lam] = eig(A);
91+lam = diag(Lam);
92+Vi = inv(V);
93+Ft = Vi * F * Vi.';
94+Ut = Ft ./ (lam * ones(1, n) + ones(n, 1) * lam.');
95+
96+U = G;
97+U(ii, ii) = real(V * Ut * V.');
98+end
src/chapters/ch07/text.mdadded+48−0View file
@@ -0,0 +1,48 @@
1+A second-order method with a perfect solver still needs $10^{12}$ unknowns for twelve digits. The way out is the same as in Chapter 3: keep the structure that makes the solver fast, and change the stencil.
2+
3+### The compact nine-point stencil
4+
5+Add the four diagonal neighbours, with weights chosen to cancel as much of the truncation error as possible:
6+
7+$$\frac{4(u_E + u_W + u_N + u_S) + (u_{NE} + u_{NW} + u_{SE} + u_{SW}) - 20\,u_P}{6h^2} \;=\; \Delta u + \frac{h^2}{12}\Delta^2 u + O(h^4).$$
8+
9+This is the Mehrstellen operator. For Poisson's equation the $h^2$ term is $\frac{h^2}{12}\Delta f$, which is known and can be subtracted, giving a fourth-order scheme with the same nine points. For **Laplace's** equation, $\Delta^2 u = \Delta(\Delta u) = 0$ identically, so that term is not merely computable, it is zero, and the next one vanishes as well: the scheme is sixth-order accurate.
10+
11+The stencil is still a sum of tensor products of tridiagonal operators, so the sine vectors are still its eigenvectors. With $a = \pi k/N$ and $b = \pi l/N$ the eigenvalue is
12+
13+$$20 - 8\cos a - 8\cos b - 4\cos a \cos b,$$
14+
15+and the solver of the previous chapter applies verbatim, with one line changed. A sixth-order method for the price of a second-order one, on the same grid, in the same time.
16+
17+### Chebyshev in two dimensions
18+
19+The one-dimensional collocation of Chapter 4 extends to the square by taking the tensor product of two Chebyshev grids. Written on the array $U$ of interior values, the discrete problem is a Sylvester equation,
20+
21+$$A U + U A^{\!\top} = F,$$
22+
23+where $A$ is the one-dimensional operator on interior nodes and $F$ collects the boundary contributions. Solving it by forming $I \otimes A + A \otimes I$ and calling backslash costs $O(N^6)$ and is the two-dimensional dense solve all over again. Diagonalising $A = V\Lambda V^{-1}$ once instead turns the equation into a division,
24+
25+$$\widetilde U_{ij} \;=\; \frac{\big(V^{-1} F V^{-\top}\big)_{ij}}{\lambda_i + \lambda_j}, \qquad U = V \widetilde U V^{\!\top},$$
26+
27+at a cost of one eigendecomposition and four matrix products, $O(N^3)$ in total. This is the matrix-decomposition method, and it is the same idea as the sine transform of Chapter 6, with the eigenvectors computed at run time instead of known in advance. Unlike the sine transform it does not care that the operator is a difference operator, only that the problem separates.
28+
29+<!--lab:fields-->
30+
31+<!--lab:sweep-->
32+
33+### What the comparison says
34+
35+On `expsin`, at the default settings:
36+
37+- the five-point scheme reaches $4 \times 10^{-8}$ with $10^{6}$ unknowns;
38+- the nine-point scheme reaches $3 \times 10^{-14}$ with $10^{3}$ unknowns, in under a millisecond;
39+- Chebyshev collocation reaches $10^{-14}$ with $225$ unknowns.
40+
41+Three methods, the same problem, the same machine, and a factor of a million in unknowns between the first and the last. The nine-point curve is a straight line of slope $-6$ until it hits the roundoff floor near $10^{-14}$; the Chebyshev curve does not have a slope, it falls off a cliff.
42+
43+The qualifications matter as much as the numbers.
44+
45+- **The sixth order is specific to Laplace's equation.** Solve Poisson's equation with the same stencil and the order drops to four unless the source term is averaged too.
46+- **Everything here needs the square.** All three methods are tensor-product constructions. Chapter 9 removes one corner from the square and Chapter 10 replaces it with a curve, and none of these three survive that.
47+- **The roundoff floor arrives earlier for higher order.** Both high-order curves flatten near $10^{-13}$ and then creep upward as $N$ grows, exactly as in Chapters 2 and 4. There is no reason to run a sixth-order scheme past the resolution where it reaches the floor.
48+- **Chebyshev's cost grows faster.** $O(N^3)$ against $O(N^2\log N)$, and the eigendecomposition has a large constant. At the accuracies reachable here it never matters, because the required $N$ is so small, but on a problem needing $N = 500$ it would.
src/chapters/ch08/chapter.tsadded+43−0View file
@@ -0,0 +1,43 @@
1+import { createElement } from 'react'
2+import type { Chapter } from '../../ui/types.ts'
3+import driver from './driver.m?raw'
4+import script from './iterative.m?raw'
5+import text from './text.md?raw'
6+import { IterPanels, iterSummary, type IterOut } from './panels.tsx'
7+import { EXACT2D_CONTROL } from '../shared.ts'
8+
9+export const ch08: Chapter = {
10+ id: 'iterative',
11+ nav: 'Iterating instead of solving',
12+ part: 'Two dimensions',
13+ title: 'Iterating instead of solving',
14+ blurb: 'Jacobi to multigrid on the same system, and the question of when to stop.',
15+ text,
16+ labs: [
17+ {
18+ id: 'history',
19+ title: 'Five iterations on one problem',
20+ caption:
21+ 'Raise N and watch which curves move. Only the multigrid one stays where it is.',
22+ driver,
23+ script,
24+ manual: true,
25+ editorOpen: true,
26+ controls: [
27+ EXACT2D_CONTROL,
28+ { kind: 'slider', key: 'N', label: 'N', values: [32, 64, 128, 256], min: 32, max: 256 },
29+ {
30+ kind: 'slider',
31+ key: 'nit',
32+ label: 'iterations',
33+ values: [20, 60, 120, 240, 480],
34+ min: 20,
35+ max: 480,
36+ },
37+ ],
38+ defaults: { kind: 'expsin', par: 1, N: 64, nit: 120, tol: 1e-11 },
39+ render: (d: IterOut) => createElement(IterPanels, { data: d }),
40+ summary: iterSummary,
41+ },
42+ ],
43+}
src/chapters/ch08/driver.madded+87−0View file
@@ -0,0 +1,87 @@
1+% Driver for the iterative-solver chapter.
2+%
3+% The editable script supplies
4+%
5+% names = method_names()
6+% H = run_iter(mid, N, U0, UE, nit, tol)
7+%
8+% where H has, in its three rows and one column per iteration, the squared
9+% relative residual, the squared maximum error, and the elapsed milliseconds.
10+% The squares are not an oversight: taking the square roots in the MATLAB layer
11+% means either a sqrt inside the iteration loop or a call after it, and either
12+% one makes numbl compile that loop through its slow path. The page roots them.
13+%
14+% and the driver supplies the starting array U0, which holds the exact
15+% boundary data on the boundary and zeros inside, and the exact solution UE for
16+% measurement only. For reference it also solves the same system exactly with
17+% the sine transform of Chapter 6, which gives both the discretisation error,
18+% the accuracy no iterative method can beat, and the time a direct solve takes.
19+p = lp_params();
20+kind = p.kind;
21+par = p.par;
22+N = p.N;
23+names = method_names();
24+nm = numel(names);
25+
26+x = (0:N)' / N;
27+[X, Y] = meshgrid(x, x);
28+UE = lp_exact2d(kind, X, Y, par);
29+
30+U0 = zeros(N+1, N+1);
31+U0(1, :) = UE(1, :);
32+U0(N+1, :) = UE(N+1, :);
33+U0(:, 1) = UE(:, 1);
34+U0(:, N+1) = UE(:, N+1);
35+
36+% the exact solution of the discrete system, by the fast Poisson solver
37+t0 = tic;
38+Ud = drv_fast(U0, N);
39+ms_direct = toc(t0) * 1000;
40+disc_err = max(max(abs(Ud - UE)));
41+
42+nit = p.nit;
43+tol = p.tol; % hoisted: a struct field read inside a loop puts numbl's
44+ % compiler into its slow path, and the loop below is the timed one
45+R = NaN(nm, nit+1);
46+E = NaN(nm, nit+1);
47+T = NaN(nm, nit+1);
48+for k = 1:nm
49+ % the copy goes through a function call: a slice assignment written here in
50+ % the loop body would put the loop, and everything it calls, on numbl's
51+ % slow path
52+ H = run_iter(k, N, U0, UE, nit, tol);
53+ R = drv_store(R, H, 1, k);
54+ E = drv_store(E, H, 2, k);
55+ T = drv_store(T, H, 3, k);
56+end
57+
58+out = struct();
59+out.names = names;
60+out.N = N;
61+out.its = 0:nit;
62+out.R2 = lp_rows(R);
63+out.E2 = lp_rows(E);
64+out.T = lp_rows(T);
65+out.discErr = disc_err;
66+out.directMs = ms_direct;
67+out.unknowns = (N-1)^2;
68+lp_write(out);
69+
70+function A = drv_store(A, H, row, k)
71+for c = 1:size(A, 2)
72+ A(k, c) = H(row, c);
73+end
74+end
75+
76+function U = drv_fast(U, N)
77+n = N - 1;
78+B = zeros(n, n);
79+B(1, :) = B(1, :) + U(1, 2:N);
80+B(n, :) = B(n, :) + U(N+1, 2:N);
81+B(:, 1) = B(:, 1) + U(2:N, 1);
82+B(:, n) = B(:, n) + U(2:N, N+1);
83+lam = 2 * (1 - cos(pi * (1:n)' / N));
84+D = lam * ones(1, n) + ones(n, 1) * lam';
85+Bh = lp_dst1(lp_dst1(B).').';
86+U(2:N, 2:N) = lp_dst1(lp_dst1(Bh ./ D).').';
87+end
src/chapters/ch08/iterative.madded+413−0View file
@@ -0,0 +1,413 @@
1+% Five ways to solve the five-point system without factorising anything.
2+%
3+% Every method works only through the action of
4+%
5+% (A u)(i,j) = 4 u(i,j) - u(i-1,j) - u(i+1,j) - u(i,j-1) - u(i,j+1)
6+%
7+% on the interior points, with the boundary values held fixed in the array. No
8+% matrix is formed, so the storage is the grid and nothing else. What separates
9+% the methods is how fast the error goes away: Jacobi and Gauss-Seidel need
10+% O(N^2) sweeps, optimal SOR and conjugate gradients O(N), and multigrid a
11+% number of cycles that does not grow with N at all.
12+%
13+% UE is passed in only so that the true error can be recorded next to the
14+% residual. No method uses it.
15+%
16+% A note on the shape of this file. Each method has its own loop rather than
17+% one loop with a switch inside it, and multigrid keeps every level in one flat
18+% vector rather than a cell array of grids. Both choices are forced by numbl's
19+% compiler, which types a loop variable from the expressions assigned to it and
20+% drops to a slow path when a loop is not uniform. Written the obvious way this
21+% file runs about twenty times slower, which would make the timings meaningless.
22+
23+function names = method_names()
24+names = {'Jacobi', 'Gauss-Seidel', 'SOR', 'conjugate gradients', 'multigrid V-cycles'};
25+end
26+
27+function H = run_iter(mid, N, U, UE, nit, tol)
28+if mid == 1
29+ H = loop_jacobi(N, U, UE, nit, tol);
30+elseif mid == 2
31+ H = loop_sor(N, U, UE, nit, tol, 1);
32+elseif mid == 3
33+ H = loop_sor(N, U, UE, nit, tol, 2 / (1 + sin(pi/N)));
34+elseif mid == 4
35+ H = loop_cg(N, U, UE, nit, tol);
36+else
37+ H = loop_mg(N, U, UE, nit, tol);
38+end
39+end
40+
41+% ------------------------------------------------------------ the methods ---
42+
43+function H = loop_jacobi(N, U, UE, nit, tol)
44+% Every point replaced by the average of its four old neighbours. Error
45+% components with wavelength comparable to the domain are damped by a factor
46+% 1 - O(h^2) per sweep, so the sweep count grows like N^2.
47+res = NaN(1, nit+1);
48+err = NaN(1, nit+1);
49+tms = NaN(1, nit+1);
50+ss0 = res_sumsq(U, N);
51+res(1) = ss0;
52+err(1) = maxsq(U - UE);
53+tms(1) = 0;
54+stop = tol * tol * ss0;
55+tacc = 0;
56+for it = 1:nit
57+ t1 = tic;
58+ U = jacobi_sweep(U, N);
59+ tacc = tacc + toc(t1);
60+ res(it+1) = res_sumsq(U, N);
61+ err(it+1) = maxsq(U - UE);
62+ tms(it+1) = tacc * 1000;
63+ if res(it+1) < stop
64+ break;
65+ end
66+end
67+% The squares are handed back as they are; the page takes the square roots.
68+% Doing it here would mean either a sqrt inside the loop, which widens the
69+% inferred type of the history to complex and costs a factor of twenty, or a
70+% call after it, which does the same thing to the loop that filled the arrays.
71+H = [res / ss0; err; tms];
72+end
73+
74+function H = loop_sor(N, U, UE, nit, tol, omega)
75+% omega = 1 is Gauss-Seidel: new values are used as soon as they are available,
76+% which costs nothing and halves the sweep count. Over-relaxing with the
77+% optimal omega = 2/(1 + sin(pi h)) changes the rate from 1 - O(h^2) to
78+% 1 - O(h). That optimal value is known only because the eigenvalues of this
79+% particular matrix are known.
80+res = NaN(1, nit+1);
81+err = NaN(1, nit+1);
82+tms = NaN(1, nit+1);
83+ss0 = res_sumsq(U, N);
84+res(1) = ss0;
85+err(1) = maxsq(U - UE);
86+tms(1) = 0;
87+stop = tol * tol * ss0;
88+tacc = 0;
89+for it = 1:nit
90+ t1 = tic;
91+ U = sor_sweep(U, N, omega);
92+ tacc = tacc + toc(t1);
93+ res(it+1) = res_sumsq(U, N);
94+ err(it+1) = maxsq(U - UE);
95+ tms(it+1) = tacc * 1000;
96+ if res(it+1) < stop
97+ break;
98+ end
99+end
100+% The squares are handed back as they are; the page takes the square roots.
101+% Doing it here would mean either a sqrt inside the loop, which widens the
102+% inferred type of the history to complex and costs a factor of twenty, or a
103+% call after it, which does the same thing to the loop that filled the arrays.
104+H = [res / ss0; err; tms];
105+end
106+
107+function H = loop_cg(N, U, UE, nit, tol)
108+% Conjugate gradients. A is symmetric positive definite, so CG applies; it
109+% needs one operator application per iteration and four vectors of storage.
110+% Its convergence is governed by sqrt(kappa) rather than kappa, which for this
111+% matrix means O(N) iterations, like optimal SOR, but without having to know
112+% anything about the spectrum in advance.
113+res = NaN(1, nit+1);
114+err = NaN(1, nit+1);
115+tms = NaN(1, nit+1);
116+R = neg_residual(U, N); % the right-hand side is zero at interior points
117+P = R;
118+rs = dot2(R, R);
119+ss0 = rs;
120+res(1) = ss0;
121+err(1) = maxsq(U - UE);
122+tms(1) = 0;
123+stop = tol * tol * ss0;
124+tacc = 0;
125+for it = 1:nit
126+ t1 = tic;
127+ AP = applyA(P, N);
128+ al = rs / dot2(P, AP);
129+ U = U + al * P;
130+ R = R - al * AP;
131+ rsn = dot2(R, R);
132+ P = R + (rsn / rs) * P;
133+ rs = rsn;
134+ tacc = tacc + toc(t1);
135+ res(it+1) = rsn;
136+ err(it+1) = maxsq(U - UE);
137+ tms(it+1) = tacc * 1000;
138+ if res(it+1) < stop
139+ break;
140+ end
141+end
142+% The squares are handed back as they are; the page takes the square roots.
143+% Doing it here would mean either a sqrt inside the loop, which widens the
144+% inferred type of the history to complex and costs a factor of twenty, or a
145+% call after it, which does the same thing to the loop that filled the arrays.
146+H = [res / ss0; err; tms];
147+end
148+
149+function H = loop_mg(N, U, UE, nit, tol)
150+% Geometric multigrid. Gauss-Seidel is a poor solver but an excellent smoother:
151+% two sweeps remove the oscillatory part of the error and leave something
152+% smooth, which is well represented on a grid twice as coarse and much cheaper
153+% to correct there. Recursing on that idea gives a V-cycle whose cost is
154+% proportional to the number of unknowns and whose error reduction per cycle is
155+% a constant independent of N.
156+nlev = 1;
157+while mod(N, 2^nlev) == 0 && N / 2^nlev >= 2
158+ nlev = nlev + 1;
159+end
160+ns = zeros(nlev, 1);
161+offs = zeros(nlev, 1);
162+tot = 0;
163+nn = N;
164+for l = 1:nlev
165+ ns(l) = nn;
166+ offs(l) = tot;
167+ tot = tot + (nn+1)^2;
168+ nn = nn / 2;
169+end
170+
171+% one flat vector holds every level: the iterate in the first half, the
172+% right-hand side in the second
173+W = zeros(2*tot, 1);
174+W = flatten_into(W, U, 0, N);
175+WE = flatten_into(zeros((N+1)^2, 1), UE, 0, N);
176+
177+res = NaN(1, nit+1);
178+err = NaN(1, nit+1);
179+tms = NaN(1, nit+1);
180+ss0 = res_sumsq_flat(W, 0, tot, N);
181+res(1) = ss0;
182+err(1) = maxsq_flat(W, WE, N);
183+tms(1) = 0;
184+stop = tol * tol * ss0;
185+tacc = 0;
186+for it = 1:nit
187+ t1 = tic;
188+ W = vcycle(W, offs, ns, nlev, tot, 2, 2);
189+ tacc = tacc + toc(t1);
190+ res(it+1) = res_sumsq_flat(W, 0, tot, N);
191+ err(it+1) = maxsq_flat(W, WE, N);
192+ tms(it+1) = tacc * 1000;
193+ if res(it+1) < stop
194+ break;
195+ end
196+end
197+% The squares are handed back as they are; the page takes the square roots.
198+% Doing it here would mean either a sqrt inside the loop, which widens the
199+% inferred type of the history to complex and costs a factor of twenty, or a
200+% call after it, which does the same thing to the loop that filled the arrays.
201+H = [res / ss0; err; tms];
202+end
203+
204+% ---------------------------------------------------------------- kernels ---
205+
206+function U = jacobi_sweep(U, N)
207+V = U;
208+for j = 2:N
209+ for i = 2:N
210+ U(i,j) = 0.25 * (V(i-1,j) + V(i+1,j) + V(i,j-1) + V(i,j+1));
211+ end
212+end
213+end
214+
215+function U = sor_sweep(U, N, omega)
216+for j = 2:N
217+ for i = 2:N
218+ U(i,j) = (1-omega)*U(i,j) + 0.25*omega*(U(i-1,j) + U(i+1,j) + U(i,j-1) + U(i,j+1));
219+ end
220+end
221+end
222+
223+function Y = applyA(U, N)
224+Y = zeros(N+1, N+1);
225+for j = 2:N
226+ for i = 2:N
227+ Y(i,j) = 4*U(i,j) - U(i-1,j) - U(i+1,j) - U(i,j-1) - U(i,j+1);
228+ end
229+end
230+end
231+
232+function Y = neg_residual(U, N)
233+Y = zeros(N+1, N+1);
234+for j = 2:N
235+ for i = 2:N
236+ Y(i,j) = -(4*U(i,j) - U(i-1,j) - U(i+1,j) - U(i,j-1) - U(i,j+1));
237+ end
238+end
239+end
240+
241+function s = dot2(A, B)
242+s = 0;
243+n = numel(A);
244+for k = 1:n
245+ s = s + A(k) * B(k);
246+end
247+end
248+
249+function ss = res_sumsq(U, N)
250+% The sum of squares, not the norm: a square root taken here would widen the
251+% inferred type of the accumulator and slow this loop by a factor of 20.
252+ss = 0;
253+for j = 2:N
254+ for i = 2:N
255+ r = 4*U(i,j) - U(i-1,j) - U(i+1,j) - U(i,j-1) - U(i,j+1);
256+ ss = ss + r*r;
257+ end
258+end
259+end
260+
261+function m = maxsq(E)
262+m = max(max(E .* E));
263+end
264+
265+
266+% -------------------------------------------------------------- multigrid ---
267+
268+function W = flatten_into(W, U, off, N)
269+m = N + 1;
270+for j = 1:m
271+ b = off + (j-1)*m;
272+ for i = 1:m
273+ W(b+i) = U(i,j);
274+ end
275+end
276+end
277+
278+function ss = res_sumsq_flat(W, off, tot, N)
279+m = N + 1;
280+ss = 0;
281+for j = 2:N
282+ b = off + (j-1)*m;
283+ for i = 2:N
284+ r = W(tot + b + i) - (4*W(b+i) - W(b+i-1) - W(b+i+1) - W(b+i-m) - W(b+i+m));
285+ ss = ss + r*r;
286+ end
287+end
288+end
289+
290+function m = maxsq_flat(W, WE, N)
291+mm = N + 1;
292+m = 0;
293+for k = 1:mm*mm
294+ d = W(k) - WE(k);
295+ if d*d > m
296+ m = d*d;
297+ end
298+end
299+end
300+
301+function W = vcycle(W, offs, ns, nlev, tot, nu1, nu2)
302+for l = 1:nlev-1
303+ W = mg_down(W, offs(l), offs(l+1), tot, ns(l), nu1);
304+end
305+W = gs_flat(W, offs(nlev), tot, ns(nlev), 8); % the coarsest grid: one unknown
306+for l = nlev-1:-1:1
307+ W = mg_up(W, offs(l), offs(l+1), tot, ns(l), nu2);
308+end
309+end
310+
311+function W = mg_down(W, off, offc, tot, n, nu)
312+% Smooth, then carry the residual to the next grid and start the coarse
313+% iterate at zero.
314+%
315+% The Gauss-Seidel sweeps are written out here rather than delegated to
316+% gs_flat, and the residual is formed inline rather than through a helper, for
317+% the same compiler reason as everywhere else in this file: a loop that writes
318+% W by index runs fast only if W did not arrive from a function call. Calling
319+% gs_flat first and then entering the restriction loop costs a factor of 250 on
320+% this function alone, which was measured, not guessed.
321+m = n + 1;
322+for s = 1:nu
323+ for j = 2:n
324+ b = off + (j-1)*m;
325+ for i = 2:n
326+ W(b+i) = 0.25 * (W(b+i-1) + W(b+i+1) + W(b+i-m) + W(b+i+m) + W(tot+b+i));
327+ end
328+ end
329+end
330+
331+nc = n / 2;
332+mc = nc + 1;
333+for jc = 2:nc
334+ j = 2*jc - 1;
335+ for ic = 2:nc
336+ i = 2*ic - 1;
337+ k = off + (j-1)*m + i;
338+ r00 = W(tot+k) + W(k+1) + W(k-1) + W(k+m) + W(k-m) - 4*W(k);
339+ rm0 = W(tot+k-1) + W(k) + W(k-2) + W(k-1+m) + W(k-1-m) - 4*W(k-1);
340+ rp0 = W(tot+k+1) + W(k+2) + W(k) + W(k+1+m) + W(k+1-m) - 4*W(k+1);
341+ r0m = W(tot+k-m) + W(k-m+1) + W(k-m-1) + W(k) + W(k-2*m) - 4*W(k-m);
342+ r0p = W(tot+k+m) + W(k+m+1) + W(k+m-1) + W(k+2*m) + W(k) - 4*W(k+m);
343+ rmm = W(tot+k-m-1) + W(k-m) + W(k-m-2) + W(k-1) + W(k-2*m-1) - 4*W(k-m-1);
344+ rpm = W(tot+k-m+1) + W(k-m+2) + W(k-m) + W(k+1) + W(k-2*m+1) - 4*W(k-m+1);
345+ rmp = W(tot+k+m-1) + W(k+m) + W(k+m-2) + W(k+2*m-1) + W(k-1) - 4*W(k+m-1);
346+ rpp = W(tot+k+m+1) + W(k+m+2) + W(k+m) + W(k+2*m+1) + W(k+1) - 4*W(k+m+1);
347+ s = 4*r00 + 2*(rm0 + rp0 + r0m + r0p) + (rmm + rpm + rmp + rpp);
348+ % full weighting is s/16; the extra factor of 4 is because the coarse
349+ % operator is written in the same scaled form and so stands for (2h)^2
350+ % times the Laplacian. Leave it out and the correction comes back four
351+ % times too small, which shows up as a V-cycle that reduces the
352+ % residual by a factor of two instead of twenty.
353+ W(tot + offc + (jc-1)*mc + ic) = s / 4;
354+ end
355+end
356+for k = 1:mc*mc
357+ W(offc + k) = 0;
358+end
359+end
360+
361+function W = mg_up(W, off, offc, tot, n, nu)
362+% Interpolate the coarse correction, add it, smooth again. The correction is
363+% built in a scratch array first: interpolating in place would mix it with the
364+% fine iterate that is already there.
365+m = n + 1;
366+nc = n / 2;
367+mc = nc + 1;
368+Z = zeros(m, m);
369+for jc = 1:mc
370+ for ic = 1:mc
371+ Z(2*ic-1, 2*jc-1) = W(offc + (jc-1)*mc + ic);
372+ end
373+end
374+for j = 1:2:m
375+ for i = 2:2:n
376+ Z(i,j) = 0.5 * (Z(i-1,j) + Z(i+1,j));
377+ end
378+end
379+for j = 2:2:n
380+ for i = 1:m
381+ Z(i,j) = 0.5 * (Z(i,j-1) + Z(i,j+1));
382+ end
383+end
384+for j = 1:m
385+ b = off + (j-1)*m;
386+ for i = 1:m
387+ W(b+i) = W(b+i) + Z(i,j);
388+ end
389+end
390+for s = 1:nu
391+ for j = 2:n
392+ b = off + (j-1)*m;
393+ for i = 2:n
394+ W(b+i) = 0.25 * (W(b+i-1) + W(b+i+1) + W(b+i-m) + W(b+i+m) + W(tot+b+i));
395+ end
396+ end
397+end
398+end
399+
400+function W = gs_flat(W, off, tot, n, nu)
401+% Gauss-Seidel on one level of the flat hierarchy. Used only for the coarsest
402+% grid, where there is a single unknown and eight sweeps solve it exactly; the
403+% sweeps on the other levels are written out inside mg_down and mg_up.
404+m = n + 1;
405+for s = 1:nu
406+ for j = 2:n
407+ b = off + (j-1)*m;
408+ for i = 2:n
409+ W(b+i) = 0.25 * (W(b+i-1) + W(b+i+1) + W(b+i-m) + W(b+i+m) + W(tot+b+i));
410+ end
411+ end
412+end
413+end
src/chapters/ch08/panels.tsxadded+177−0View file
@@ -0,0 +1,177 @@
1+// Figures for the iterative-solver chapter: how the residual falls per
2+// iteration, how the true error falls, and how both compare with a direct
3+// solve of the same system.
4+import Plot from '../../plot/Plot.tsx'
5+import Legend from '../../plot/Legend.tsx'
6+import { EndLabel, Line, Markers } from '../../plot/Curves.tsx'
7+import { ink, marker, method, type MethodKey } from '../../plot/palette.ts'
8+import type { Num } from '../../engine/types.ts'
9+import { fmtErr } from '../onedim/panels.tsx'
10+
11+export interface IterOut {
12+ names: string[]
13+ N: number
14+ its: number[]
15+ /** squared relative residual and squared error: the MATLAB layer hands back
16+ * the squares, for the reason its driver explains */
17+ R2: Num[][]
18+ E2: Num[][]
19+ T: Num[][]
20+ discErr: number
21+ directMs: number
22+ unknowns: number
23+}
24+
25+const rootOf = (rows: Num[][]): Num[][] =>
26+ rows.map((r) => r.map((v) => (v == null || !isFinite(v) || v < 0 ? null : Math.sqrt(v))))
27+
28+// one palette slot per method, in the order the script lists them
29+const SLOTS: MethodKey[] = ['fd2', 'fd4', 'fast', 'spectral', 'iterative']
30+
31+function lowHigh(rows: Num[][], floor = 1e-14): [number, number] {
32+ let lo = Infinity
33+ let hi = -Infinity
34+ for (const r of rows) {
35+ for (const v of r) {
36+ if (v == null || !isFinite(v) || v <= 0) continue
37+ if (v < lo) lo = v
38+ if (v > hi) hi = v
39+ }
40+ }
41+ if (!isFinite(lo)) return [floor, 1]
42+ return [Math.max(lo / 3, floor), hi * 3]
43+}
44+
45+function lastIndex(row: Num[]): number {
46+ for (let i = row.length - 1; i >= 0; i--) if (row[i] != null && isFinite(row[i] as number)) return i
47+ return 0
48+}
49+
50+export function IterPanels({ data }: { data: IterOut }) {
51+ const its = data.its
52+ const R = rootOf(data.R2)
53+ const E = rootOf(data.E2)
54+ const xd: [number, number] = [0, its[its.length - 1]]
55+ const rd = lowHigh(R)
56+ const ed = lowHigh([...E, [data.discErr]])
57+ const td = lowHigh(data.T.map((r) => r.map((v) => (v != null && v > 0 ? v : null))))
58+ td[0] = Math.max(td[0], 0.05)
59+
60+ const legend = data.names.map((n, k) => ({ label: n, color: method[SLOTS[k]] }))
61+
62+ return (
63+ <>
64+ <div className="plot-title">relative residual against iteration count</div>
65+ <Plot height={300} xDomain={xd} yDomain={rd} yLog xLabel="iteration" yLabel="‖r‖ / ‖r₀‖">
66+ {(f) => (
67+ <>
68+ {R.map((row, k) => (
69+ <g key={k}>
70+ <Line x={its} y={row} f={f} color={method[SLOTS[k]]} />
71+ {row.length <= 40 && (
72+ <Markers x={its} y={row} f={f} color={method[SLOTS[k]]} shape={marker[SLOTS[k]]} r={3} />
73+ )}
74+ </g>
75+ ))}
76+ </>
77+ )}
78+ </Plot>
79+ <Legend items={legend} />
80+
81+ <div className="plot-title">
82+ true error against iteration count; the dashed line is the discretisation error, which is where the exact
83+ solution of this linear system sits
84+ </div>
85+ <Plot
86+ height={300}
87+ xDomain={xd}
88+ yDomain={ed}
89+ yLog
90+ xLabel="iteration"
91+ yLabel="max |error|"
92+ margin={{ r: 30 }}
93+ over={(f) => (
94+ <>
95+ {E.map((row, k) => {
96+ const i = lastIndex(row)
97+ return (
98+ <EndLabel
99+ key={k}
100+ f={f}
101+ x={its[i]}
102+ y={row[i] as number}
103+ text={data.names[k]}
104+ color={method[SLOTS[k]]}
105+ anchor="end"
106+ dx={-6}
107+ dy={-7}
108+ />
109+ )
110+ })}
111+ </>
112+ )}
113+ >
114+ {(f) => (
115+ <>
116+ <Line x={xd} y={[data.discErr, data.discErr]} f={f} color={ink.muted} dash="6 4" width={1.5} />
117+ {E.map((row, k) => (
118+ <Line key={k} x={its} y={row} f={f} color={method[SLOTS[k]]} />
119+ ))}
120+ </>
121+ )}
122+ </Plot>
123+ <Legend
124+ items={[...legend, { label: 'discretisation error', color: ink.muted, dash: '6 4' }]}
125+ />
126+
127+ <div className="plot-title">
128+ true error against milliseconds of solver time; the marked point is the sine-transform direct solve of the same
129+ system
130+ </div>
131+ <Plot height={300} xDomain={td} yDomain={ed} xLog yLog xLabel="milliseconds" yLabel="max |error|">
132+ {(f) => (
133+ <>
134+ {E.map((row, k) => (
135+ <Line
136+ key={k}
137+ x={(data.T[k] ?? []).map((v) => (v == null || v <= 0 ? NaN : v))}
138+ y={row}
139+ f={f}
140+ color={method[SLOTS[k]]}
141+ />
142+ ))}
143+ <Markers
144+ x={[data.directMs]}
145+ y={[data.discErr]}
146+ f={f}
147+ color={ink.node}
148+ shape="star"
149+ r={7}
150+ />
151+ <EndLabel
152+ f={f}
153+ x={data.directMs}
154+ y={data.discErr}
155+ text="direct solve"
156+ color={ink.secondary}
157+ dx={-12}
158+ dy={-10}
159+ anchor="end"
160+ />
161+ </>
162+ )}
163+ </Plot>
164+ <Legend items={[...legend, { label: 'direct solve', color: ink.node, shape: 'cross' }]} />
165+ </>
166+ )
167+}
168+
169+export function iterSummary(d: IterOut) {
170+ const mg = d.E2[d.E2.length - 1]
171+ const i = lastIndex(mg)
172+ return [
173+ { label: 'unknowns', value: String(d.unknowns) },
174+ { label: 'discretisation error', value: fmtErr(d.discErr) },
175+ { label: 'multigrid cycles used', value: String(d.its[i]) },
176+ ]
177+}
src/chapters/ch08/text.mdadded+39−0View file
@@ -0,0 +1,39 @@
1+The sine transform solved the five-point system in $O(N^2\log N)$, but it needed a rectangle, a uniform grid and constant coefficients. Without those, and with $N^2$ unknowns too many to eliminate, the remaining option is to iterate: start from a guess and improve it, using only the ability to apply the operator.
2+
3+This chapter is about how much difference the choice of iteration makes, and about a question that only appears once one iterates: when to stop.
4+
5+### Five iterations
6+
7+**Jacobi.** Replace every value by the average of its four old neighbours. It is the simplest thing that could work, and it converges: the eigenvalues of the iteration matrix are $\frac{1}{2}(\cos\pi k h + \cos \pi l h)$, whose largest magnitude is $1 - \pi^2h^2/2 + O(h^4)$. Reducing the error by a fixed factor therefore takes $O(N^2)$ sweeps, and since a sweep costs $O(N^2)$, the total is $O(N^4)$: worse than the dense solve is not, but not by much.
8+
9+**Gauss-Seidel.** Use each new value as soon as it is computed. No extra storage, no extra work, and the asymptotic rate is exactly the square of Jacobi's, which halves the number of sweeps. Still $O(N^2)$ sweeps.
10+
11+**SOR.** Overshoot the Gauss-Seidel update by a factor $\omega$. For this matrix the optimal choice is known,
12+
13+$$\omega_{\text{opt}} = \frac{2}{1 + \sin \pi h}, \qquad \rho = \omega_{\text{opt}} - 1 = 1 - 2\pi h + O(h^2),$$
14+
15+which turns $O(N^2)$ sweeps into $O(N)$. That is a genuine change of order, and it is also the method's weakness: $\omega_{\text{opt}}$ is known here only because the spectrum of this particular matrix is known in closed form, and the rate degrades quickly if $\omega$ is wrong in the wrong direction.
16+
17+**Conjugate gradients.** The matrix is symmetric positive definite, so CG applies. It needs one operator application and a few vector operations per step, and its convergence depends on $\sqrt{\kappa} = O(N)$ rather than on $\kappa$, giving $O(N)$ iterations like optimal SOR, but with no parameter to choose and no knowledge of the spectrum required. With a good preconditioner it does better; the labs here use none.
18+
19+**Multigrid.** Gauss-Seidel is a bad solver but a good *smoother*: two sweeps remove the oscillatory part of the error and leave something smooth. Smooth functions are well represented on a grid twice as coarse, where the remaining problem costs a quarter as much to work on, and the same argument applies again there. The resulting V-cycle costs $O(N^2)$, proportional to the number of unknowns, and reduces the error by a constant factor, independent of $N$, each time. It is the only method here whose iteration count does not grow with the mesh.
20+
21+<!--lab:history-->
22+
23+### Reading the three figures
24+
25+The first shows the residual. On $N = 64$ the classical relaxations have barely started after a hundred and twenty sweeps, SOR is doing what its rate predicts, CG is somewhat better than SOR, and multigrid has reached $10^{-12}$ in ten cycles, at which point it stops because it has hit the tolerance. Raise $N$ to 128 and everything except multigrid gets worse, exactly in proportion to the rates quoted above; multigrid does not notice.
26+
27+The second shows the true error, which is the quantity anybody actually cares about, together with the dashed line marking the error of the *exact* solution of the discrete system. Every curve levels off at that line and cannot go below it. That is the practical point of the chapter:
28+
29+> Iterating past the discretisation error is wasted work. It buys a more accurate solution of the wrong problem.
30+
31+A useful stopping rule follows. The discretisation error here is $O(h^2)$, so a relative residual tolerance of about $h^2$ is enough, and anything smaller is spent for nothing. The labs stop at a fixed relative residual so that the methods can be compared, and the second figure shows how much of that work each method wasted.
32+
33+The third figure puts the same errors against measured milliseconds, and adds the sine-transform solve of Chapter 6 as a single point. On the square, with a uniform grid, the direct solver wins outright and no iteration is worth considering. That is the honest comparison: the iterative methods in this chapter earn their place on the problems where the direct solver is unavailable, which is where the next two chapters go.
34+
35+### What carries over, and what does not
36+
37+Of these five, only three would survive a change of domain or of coefficients unchanged. Jacobi, Gauss-Seidel and CG need nothing but the ability to apply the operator, so they run on an L-shaped grid, on a variable-coefficient problem, on an unstructured mesh. SOR needs a spectrum it will not have in closed form. Multigrid needs a hierarchy of grids and transfer operators between them, which is straightforward on a rectangle, real work on a complicated domain, and the subject of a large literature (algebraic multigrid) when there is no geometry to hierarchise.
38+
39+The convergence *rates* also carry over, with the caveat that they were derived for this operator on this domain. A re-entrant corner, for instance, changes the constants without changing the orders, which the next chapter makes visible in a different way.
src/chapters/ch09/chapter.tsadded+103−0View file
@@ -0,0 +1,103 @@
1+import { createElement, Fragment } from 'react'
2+import type { Chapter, Control } from '../../ui/types.ts'
3+import driver from './driver.m?raw'
4+import script from './lshape.m?raw'
5+import text from './text.md?raw'
6+import { Fields2D, type Single2DOut } from '../twodim/panels.tsx'
7+import { Convergence, WorkPrecision, fmtErr, type SweepOut } from '../onedim/panels.tsx'
8+import type { MethodKey } from '../../plot/palette.ts'
9+
10+const SLOTS: MethodKey[] = ['fd2', 'mps']
11+
12+const KIND: Control = {
13+ kind: 'select',
14+ key: 'kind',
15+ label: 'ground truth',
16+ options: [
17+ { value: 'both', label: 'singular + smooth (generic)' },
18+ { value: 'sing', label: 'singular only: r^(2/3) sin(2θ/3)' },
19+ { value: 'smooth', label: 'smooth only: log|z − z₀|' },
20+ ],
21+}
22+
23+export const ch09: Chapter = {
24+ id: 'corner',
25+ nav: 'A corner, and what it costs',
26+ part: 'Beyond the square',
27+ title: 'A corner, and what it costs',
28+ blurb: 'The L-shaped domain, the loss of second order, and a basis that contains the singularity.',
29+ text,
30+ labs: [
31+ {
32+ id: 'field',
33+ title: 'The solution and where the error lives',
34+ caption: 'Compare the error picture on the generic case with the one on the smooth case.',
35+ driver,
36+ script,
37+ manual: true,
38+ controls: [
39+ KIND,
40+ {
41+ kind: 'select',
42+ key: 'mid',
43+ label: 'method',
44+ options: [
45+ { value: 1, label: '5-point + CG' },
46+ { value: 2, label: 'harmonic basis' },
47+ ],
48+ },
49+ { kind: 'slider', key: 'n', label: 'resolution', values: [32, 64, 96, 128], min: 32, max: 128 },
50+ ],
51+ defaults: { kind: 'both', mode: 'field', mid: 1, n: 64 },
52+ derive: (p) => ({ ...p, n: Number(p.mid) === 2 ? Math.max(8, Math.round(Number(p.n) / 4)) : p.n }),
53+ render: (d: Single2DOut) =>
54+ createElement(Fields2D, { data: d, xDomain: [-1, 1], yDomain: [-1, 1] }),
55+ summary: (d: Single2DOut, p) => [
56+ {
57+ label: 'resolution',
58+ value:
59+ Number(p.mid) === 2
60+ ? `${Math.max(8, Math.round(Number(p.n) / 4))} basis terms`
61+ : `${Number(p.n)} per unit`,
62+ },
63+ { label: 'max error', value: fmtErr(d.maxerr) },
64+ ],
65+ },
66+ {
67+ id: 'sweep',
68+ title: 'A grid against a basis',
69+ caption:
70+ 'The two methods have unrelated resolution parameters, so the first plot is only meaningful curve by curve; the second, against time, compares them directly.',
71+ driver,
72+ script,
73+ manual: true,
74+ editorOpen: true,
75+ controls: [KIND],
76+ defaults: {
77+ kind: 'both',
78+ mode: 'sweep',
79+ nsm: [
80+ [16, 32, 64, 128],
81+ [8, 16, 24, 40],
82+ ],
83+ },
84+ render: (d: SweepOut) =>
85+ createElement(
86+ Fragment,
87+ null,
88+ createElement(Convergence, {
89+ key: 'c',
90+ data: d,
91+ style: { slots: SLOTS, refs: [{ p: -4 / 3, label: 'slope −4/3' }] },
92+ height: 320,
93+ xLabel: 'grid intervals per unit, or basis terms',
94+ }),
95+ createElement(WorkPrecision, { key: 'w', data: d, style: { slots: SLOTS }, height: 320 }),
96+ ),
97+ summary: (d: SweepOut) => [
98+ { label: 'best difference', value: fmtErr(Math.min(...(d.E[0].filter((v) => v != null) as number[]))) },
99+ { label: 'best basis', value: fmtErr(Math.min(...(d.E[1].filter((v) => v != null) as number[]))) },
100+ ],
101+ },
102+ ],
103+}
src/chapters/ch09/driver.madded+119−0View file
@@ -0,0 +1,119 @@
1+% Driver for the L-shaped domain.
2+%
3+% The editable script supplies
4+%
5+% names = method_names()
6+% v = solve_l(mid, n, g, px, py)
7+%
8+% and must return the computed solution at the points (px, py), whatever grid
9+% or basis it used internally. Comparing two methods at a common set of points,
10+% rather than each on its own grid, is what makes a finite difference scheme
11+% and a basis expansion comparable at all.
12+%
13+% g(x, y) is the Dirichlet data, which here is the exact solution restricted to
14+% the boundary. n is the method's own resolution parameter: intervals per unit
15+% length for the difference scheme, number of basis terms for the expansion.
16+p = lp_params();
17+kind = p.kind;
18+names = method_names();
19+nm = numel(names);
20+g = @(xx, yy) lp_lexact(kind, xx, yy);
21+
22+out = struct();
23+out.names = names;
24+
25+if strcmp(p.mode, 'field')
26+ % a picture on a fixed grid of spacing 1/32, masked outside the domain
27+ m = 32;
28+ xs = (-m:m)' / m;
29+ [X, Y] = meshgrid(xs, xs);
30+ inside = ~(X > 1e-12 & Y < -1e-12);
31+ px = X(inside);
32+ py = Y(inside);
33+ t0 = tic;
34+ v = solve_l(p.mid, p.n, g, px, py);
35+ out.ms = toc(t0) * 1000;
36+ ve = lp_lexact(kind, px, py);
37+ out.maxerr = max(abs(v - ve));
38+
39+ U = NaN(2*m+1, 2*m+1);
40+ E = NaN(2*m+1, 2*m+1);
41+ U(inside) = v;
42+ E(inside) = abs(v - ve);
43+ out.x = xs';
44+ out.y = xs';
45+ out.U = lp_rows(U);
46+ out.Err = lp_rows(E);
47+else
48+ % errors at the nodes of a fixed grid of spacing 1/8, well away from the
49+ % plotting grid and coarse enough that every difference grid contains it
50+ m = 8;
51+ xs = (-m:m)' / m;
52+ [X, Y] = meshgrid(xs, xs);
53+ inside = ~(X > 1e-12 & Y < -1e-12) & abs(X) < 1 - 1e-12 & abs(Y) < 1 - 1e-12;
54+ px = X(inside);
55+ py = Y(inside);
56+ ve = lp_lexact(kind, px, py);
57+
58+ % each method gets its own list of resolutions: intervals per unit length
59+ % for the difference scheme, basis terms for the expansion
60+ % numbl's jsondecode turns a JSON array of arrays into a cell of column
61+ % vectors, where MATLAB would give a matrix; unpack it by hand
62+ NSc = p.nsm;
63+ nn = numel(NSc{1});
64+ NS = zeros(nm, nn);
65+ for k = 1:nm
66+ row = NSc{k};
67+ NS(k, :) = row(:)';
68+ end
69+ E = zeros(nm, nn);
70+ T = zeros(nm, nn);
71+ Q = zeros(nm, nn);
72+ ETall = drv_methodsl(1, NaN(2*nm, nn), NS, g, px, py, ve, nm);
73+ for k = 1:nm
74+ E(k, :) = ETall(2*k-1, :);
75+ T(k, :) = ETall(2*k, :);
76+ end
77+ for k = 1:nm
78+ Q(k, :) = lp_orders(E(k, :), NS(k, :));
79+ end
80+ out.ns = NS(1, :);
81+ out.NS = lp_rows(NS);
82+ out.E = lp_rows(E);
83+ out.T = lp_rows(T);
84+ out.orders = lp_rows(Q);
85+end
86+
87+lp_write(out);
88+
89+function ET = drv_methodsl(k, ET, NS, g, px, py, ve, nm)
90+% Recursion rather than a loop: see the note in the one-dimensional driver.
91+if k > nm
92+ return;
93+end
94+vw = solve_l(k, NS(k, 1), g, px(1), py(1)); % warm-up
95+R = drv_resl(1, NaN(2, size(NS, 2)), k, NS(k, :), g, px, py, ve);
96+ET(2*k-1, :) = R(1, :);
97+ET(2*k, :) = R(2, :);
98+ET = drv_methodsl(k+1, ET, NS, g, px, py, ve, nm);
99+end
100+
101+function R = drv_resl(j, R, mid, ns, g, px, py, ve)
102+if j > numel(ns)
103+ return;
104+end
105+n = ns(j);
106+if isfinite(n) && n > 0
107+ t0 = tic;
108+ v = solve_l(mid, n, g, px, py);
109+ R(2, j) = toc(t0) * 1000;
110+ R(1, j) = max(abs(v - ve));
111+end
112+R = drv_resl(j+1, R, mid, ns, g, px, py, ve);
113+end
114+
115+function A = drv_store(A, ET, row, k)
116+for c = 1:size(A, 2)
117+ A(k, c) = ET(row, c);
118+end
119+end
src/chapters/ch09/lshape.madded+188−0View file
@@ -0,0 +1,188 @@
1+% Two ways to solve the Dirichlet problem on the L-shaped domain.
2+%
3+% 1. the five-point scheme on a uniform grid, solved by conjugate gradients
4+% until the residual has fallen by 1e-9, which is far below the
5+% discretisation error, so that what is measured is the discretisation
6+% and not the solver;
7+% 2. an expansion in harmonic functions, fitted to the boundary data in the
8+% least-squares sense: the method of particular solutions.
9+%
10+% The second method's basis is the interesting part. Harmonic polynomials
11+% Re z^m and Im z^m handle everything that is smooth up to the boundary, and
12+% they alone would converge slowly, because the solution is not smooth at the
13+% re-entrant corner. Adding the corner's own singular functions
14+%
15+% r^(2k/3) sin(2k theta / 3), k = 1, 2, 4, 5, 7, ...
16+%
17+% (the multiples of 3 are omitted because they are already harmonic
18+% polynomials) puts the singularity into the basis instead of trying to
19+% approximate it. This is the Fox-Henrici-Moler construction, and it is why a
20+% few dozen basis functions can beat a grid of a hundred thousand points.
21+
22+function names = method_names()
23+names = {'5-point + CG', 'harmonic basis, least squares'};
24+end
25+
26+function v = solve_l(mid, n, g, px, py)
27+if mid == 1
28+ v = fd_l(n, g, px, py);
29+else
30+ v = mps_l(n, g, px, py);
31+end
32+end
33+
34+% --------------------------------------------------- finite differences -----
35+
36+function v = fd_l(N, g, px, py)
37+% N intervals per unit length, so the grid is (-1 : 1/N : 1) in both
38+% directions and every corner of the domain is a grid point.
39+h = 1 / N;
40+m = 2*N + 1;
41+xs = (-N:N)' * h;
42+[X, Y] = meshgrid(xs, xs);
43+
44+inL = ~(X > 1e-12 & Y < -1e-12); % the closed domain
45+% a node is an unknown when it lies in the open domain: inside the square, in
46+% the L, and not on either edge of the slit
47+onslit = (abs(Y) < 1e-12 & X > -1e-12) | (abs(X) < 1e-12 & Y < 1e-12);
48+unknown = inL & abs(X) < 1 - 1e-12 & abs(Y) < 1 - 1e-12 & ~onslit;
49+mask = double(unknown);
50+
51+U = zeros(m, m);
52+U(~unknown) = g(X(~unknown), Y(~unknown)); % boundary values
53+U(~inL) = 0; % outside: never read
54+
55+% Note the fresh name: assigning the result of cg_solve back into U, which was
56+% built here with zeros(), would give U two different inferred types in this
57+% function, and numbl would then pass it to cg_solve as an untyped value and
58+% run the whole solve about twenty times slower.
59+S = cg_solve(U, mask, m, 1e-9, 4000);
60+
61+% the requested points are grid points, so sampling is a lookup
62+v = zeros(numel(px), 1);
63+for k = 1:numel(px)
64+ i = round(py(k)/h) + N + 1;
65+ j = round(px(k)/h) + N + 1;
66+ v(k) = S(i, j);
67+end
68+end
69+
70+function U = cg_solve(U, mask, m, tol, maxit)
71+R = neg_residual(U, mask, m);
72+P = R;
73+rs = dot2(R, R);
74+stop = tol * tol * rs;
75+for it = 1:maxit
76+ AP = applyA(P, mask, m);
77+ al = rs / dot2(P, AP);
78+ U = U + al * P;
79+ R = R - al * AP;
80+ rsn = dot2(R, R);
81+ if rsn < stop
82+ break;
83+ end
84+ P = R + (rsn / rs) * P;
85+ rs = rsn;
86+end
87+end
88+
89+function Y = applyA(V, mask, m)
90+Y = zeros(m, m);
91+for j = 2:m-1
92+ for i = 2:m-1
93+ Y(i,j) = mask(i,j) * (4*V(i,j) - V(i-1,j) - V(i+1,j) - V(i,j-1) - V(i,j+1));
94+ end
95+end
96+end
97+
98+function Y = neg_residual(V, mask, m)
99+Y = zeros(m, m);
100+for j = 2:m-1
101+ for i = 2:m-1
102+ Y(i,j) = -mask(i,j) * (4*V(i,j) - V(i-1,j) - V(i+1,j) - V(i,j-1) - V(i,j+1));
103+ end
104+end
105+end
106+
107+function s = dot2(A, B)
108+s = 0;
109+n = numel(A);
110+for k = 1:n
111+ s = s + A(k) * B(k);
112+end
113+end
114+
115+% ------------------------------------------- the method of particular sols --
116+
117+function v = mps_l(n, g, px, py)
118+M = n; % harmonic polynomials up to degree M
119+K = n; % corner functions r^(2k/3) sin(2k theta/3)
120+
121+nb = (2*M + 1) + K;
122+[bx, by] = boundary_points(max(8, ceil(nb / 2))); % about three points per unknown
123+A = basis_matrix(bx, by, M, K);
124+c = A \ g(bx, by); % least squares: many more points than unknowns
125+
126+v = basis_matrix(px, py, M, K) * c;
127+end
128+
129+function A = basis_matrix(x, y, M, K)
130+np = numel(x);
131+nb = (2*M + 1) + numel(corner_orders(K));
132+A = zeros(np, nb);
133+z = x(:) + 1i*y(:);
134+col = 1;
135+p = ones(np, 1);
136+for mdeg = 0:M
137+ A(:, col) = real(p);
138+ col = col + 1;
139+ if mdeg > 0
140+ A(:, col) = imag(p);
141+ col = col + 1;
142+ end
143+ p = p .* z;
144+end
145+th = atan2(y(:), x(:));
146+th(th < -1e-13) = th(th < -1e-13) + 2*pi;
147+r = abs(z);
148+ks = corner_orders(K);
149+for q = 1:numel(ks)
150+ k = ks(q);
151+ A(:, col) = r.^(2*k/3) .* sin(2*k*th/3);
152+ col = col + 1;
153+end
154+end
155+
156+function ks = corner_orders(K)
157+% k = 1, 2, 4, 5, 7, ...: the multiples of 3 are harmonic polynomials already
158+% in the basis, and including them twice would make the least-squares matrix
159+% exactly rank deficient.
160+ks = [];
161+k = 1;
162+while numel(ks) < K
163+ if mod(k, 3) ~= 0
164+ ks = [ks, k];
165+ end
166+ k = k + 1;
167+end
168+end
169+
170+function [bx, by] = boundary_points(nper)
171+% Points along the six edges of the boundary, walking round from the origin:
172+%
173+% 1: y = 0, x from 0 to 1 the upper edge of the slit
174+% 2: x = 1, y from 0 to 1
175+% 3: y = 1, x from 1 to -1
176+% 4: x = -1, y from 1 to -1
177+% 5: y = -1, x from -1 to 0
178+% 6: x = 0, y from -1 to 0 the left edge of the slit
179+%
180+% The two edges that meet at the re-entrant corner carry points clustered
181+% towards it, where the data varies fastest.
182+t = ((1:nper)' - 0.5) / nper;
183+c = t .^ 2;
184+z = zeros(nper, 1);
185+o = ones(nper, 1);
186+bx = [c; o; 1 - 2*t; -o; -1 + t; z];
187+by = [z; t; o; 1 - 2*t; -o; -c];
188+end
src/chapters/ch09/text.mdadded+59−0View file
@@ -0,0 +1,59 @@
1+Every method so far has been a tensor product on a square, and every convergence rate has come from a Taylor expansion of a smooth solution. This chapter removes one quarter of the square and watches both assumptions fail.
2+
3+The domain is the L,
4+
5+$$\Omega = (-1,1)^2 \setminus \{x > 0,\ y < 0\},$$
6+
7+with a re-entrant corner at the origin of interior angle $3\pi/2$. It is the standard test problem for corner singularities, and its behaviour is representative of any domain with a non-convex corner.
8+
9+### What the corner does
10+
11+Near a corner of interior angle $\alpha$, harmonic functions with zero data on both edges behave like $r^{\pi/\alpha}$. For $\alpha = 3\pi/2$ the exponent is $2/3$, so the local behaviour is spanned by
12+
13+$$s_k(r,\theta) = r^{2k/3}\sin\frac{2k\theta}{3}, \qquad k = 1, 2, 3, \dots$$
14+
15+measuring $\theta$ from the edge $y = 0$, $x > 0$. Each $s_k$ is harmonic and vanishes on both edges of the corner. The first one, $s_1 = r^{2/3}\sin(2\theta/3)$, is continuous but has an infinite gradient at the origin and no second derivative there.
16+
17+This matters because a generic Dirichlet datum on $\partial\Omega$ produces a solution containing $s_1$. Not approximately: the solution genuinely has that term, with some coefficient determined non-locally by the data on the whole boundary. The finite difference scheme is then being asked to approximate a function that has no fourth derivative, no third, and not even a bounded first derivative, and the argument of Chapter 1 collapses.
18+
19+Three ground truths are available in the labs, so that the effect can be isolated:
20+
21+- `sing`, the singular function $s_1$ alone;
22+- `smooth`, $\log|z - z_0|$ with $z_0$ outside, which is harmonic and smooth right up to and including the corner;
23+- `both`, the sum, which is the generic case.
24+
25+<!--lab:field-->
26+
27+Look at the error picture on `both` with the difference method. The error is not spread over the domain; it is concentrated in a small neighbourhood of the origin, and refining the grid shrinks that neighbourhood without reducing the error in it very much. The `smooth` case, on the same domain and the same grid, has a smooth error of the usual shape. The corner is not the problem; the singular function is.
28+
29+### The rate
30+
31+<!--lab:sweep-->
32+
33+Three numbers from the sweep, on the maximum error over interior points at a distance of at least $1/8$ from the corner:
34+
35+| ground truth | observed order of the difference scheme |
36+|---|---|
37+| `smooth` | 2.00 |
38+| `both` | about 1.36 |
39+| `sing` | about 1.40 |
40+
41+The classical result for the five-point scheme on this domain is $O(h^{4/3})$ away from the corner, which is what the last two rows show. Right at the corner the error is worse, $O(h^{2/3})$. Either way, second order is gone, and no solver recovers it: the discretisation, not the linear algebra, is the limit.
42+
43+The standard remedies are worth naming even though the labs here take a different route. **Local refinement** with a grid graded like $h_j \sim j^{-3/2}$ towards the corner restores $O(h^2)$ at the cost of an unstructured or hierarchical grid, and is what an adaptive finite element code does automatically. **Singular enrichment** adds $s_1$ to the finite element space directly. Both are more work than the method below, on this problem.
44+
45+### Putting the singularity in the basis
46+
47+If the solution is a smooth harmonic function plus a few multiples of $s_k$, the obvious thing to do is to build an approximation out of exactly that. Take
48+
49+$$u \;\approx\; \sum_{m=0}^{M} \big(a_m \operatorname{Re} z^m + b_m \operatorname{Im} z^m\big) \;+\; \sum_{k} c_k\, r^{2k/3}\sin\frac{2k\theta}{3},$$
50+
51+every term of which is harmonic, so any combination satisfies the equation *exactly* in the interior. Nothing is left but the boundary condition, and that is imposed by least squares at a few hundred boundary points. This is the method of particular solutions, in the form Fox, Henrici and Moler used for the L-shaped eigenvalue problem in 1967, and it is the natural method for a problem like this one.
52+
53+The results in the sweep are of a different kind from anything else in this tutorial. On `sing`, the exact solution *is* the first basis function, so the method is exact to roundoff with eight terms. On `both` it converges geometrically, reaching $1.6\times10^{-14}$ with 121 basis functions in about ten milliseconds, while the difference scheme with 49000 unknowns and four seconds of conjugate gradients is at $4\times10^{-4}$. That is ten orders of magnitude and a factor of four hundred in time, in the same direction.
54+
55+Three honest qualifications, because a gap that size always has conditions attached.
56+
57+- **It needs the equation to be homogeneous and linear with known solutions.** Every basis function must be harmonic. With a source term, a variable coefficient, or a nonlinearity, this method does not exist in this form.
58+- **It needs to know where the singularities are.** The exponents $2k/3$ come from the corner angle. On a domain with several re-entrant corners each one contributes its own family, which is manageable; on a problem where the singular behaviour is not known analytically it is not.
59+- **The least-squares matrix is ill-conditioned**, with a condition number that grows quickly in the number of terms. It is usable to about $10^{-14}$ here because the boundary data is smooth and the system is solved by QR; pushed further it stalls, and the stall is not visible in the residual of the fit. Adding terms past the point where the error stops improving makes things slowly worse.
src/chapters/ch10/bie.madded+105−0View file
@@ -0,0 +1,105 @@
1+% Two boundary methods for the interior Dirichlet problem on a smooth curve.
2+%
3+% 1. NYSTROM, double layer. Look for the solution as the potential of a dipole
4+% density sigma spread over the curve,
5+%
6+% u(x) = (1/2pi) INT_Gamma (n_y . (x - y)) / |x - y|^2 sigma(y) ds(y),
7+%
8+% which is harmonic inside for any sigma. Letting x approach the curve from
9+% inside picks up a jump, and the boundary condition becomes a Fredholm
10+% equation of the second kind,
11+%
12+% -sigma(x)/2 + INT_Gamma K(x, y) sigma(y) ds(y) = g(x).
13+%
14+% Second kind is the important phrase: the operator is identity plus
15+% compact, its eigenvalues cluster away from zero, and the discrete system
16+% is well conditioned no matter how fine the discretisation. Compare the
17+% finite difference matrices, whose condition number grows like N^2.
18+%
19+% The kernel K looks singular but is not: as y approaches x it tends to
20+% -kappa(x)/(4pi), where kappa is the curvature. So the integrand is smooth
21+% and periodic, and the plain trapezoid rule, which is spectrally accurate
22+% for smooth periodic integrands, converges geometrically. That single fact
23+% is what makes this method fast.
24+%
25+% 2. FUNDAMENTAL SOLUTIONS. Place point sources on a circle of radius par
26+% outside the curve and fit their combined field to the boundary data by
27+% least squares. Every basis function is harmonic, so the equation is
28+% satisfied exactly and only the boundary condition is approximated. It is
29+% three lines of code and it also converges geometrically, but the matrix is
30+% exponentially ill conditioned in the number of sources, and the accuracy
31+% it reaches depends on where the sources are put.
32+
33+function names = method_names()
34+names = {'double-layer Nystrom', 'fundamental solutions'};
35+end
36+
37+function v = solve_bie(mid, n, a, nlobe, par, px, py)
38+if mid == 1
39+ v = nystrom(n, a, nlobe, px, py);
40+else
41+ v = mfs(n, a, nlobe, par, px, py);
42+end
43+end
44+
45+function v = nystrom(n, a, nlobe, px, py)
46+t = 2*pi*(0:n-1)' / n;
47+h = 2*pi / n;
48+C = lp_curve(a, nlobe, t);
49+x = C(:,1);
50+y = C(:,2);
51+xp = C(:,3);
52+yp = C(:,4);
53+xpp = C(:,5);
54+ypp = C(:,6);
55+
56+sp = sqrt(xp.^2 + yp.^2); % speed |gamma'(t)|
57+nx = yp ./ sp; % outward normal, counterclockwise curve
58+ny = -xp ./ sp;
59+kappa = (xp .* ypp - yp .* xpp) ./ sp.^3;
60+
61+% the Nystrom matrix: A = -I/2 + h * K * diag(speed)
62+DX = x * ones(1, n) - ones(n, 1) * x';
63+DY = y * ones(1, n) - ones(n, 1) * y';
64+R2 = DX.^2 + DY.^2;
65+K = ((ones(n, 1) * nx') .* DX + (ones(n, 1) * ny') .* DY) ./ (2*pi*R2);
66+K(1:n+1:end) = -kappa / (4*pi); % the diagonal limit, not a singularity
67+A = -0.5*eye(n) + h * K .* (ones(n, 1) * sp');
68+
69+sigma = A \ lp_bexact(x, y);
70+
71+% evaluate the layer potential at the target points
72+np = numel(px);
73+EX = px(:) * ones(1, n) - ones(np, 1) * x';
74+EY = py(:) * ones(1, n) - ones(np, 1) * y';
75+KE = ((ones(np, 1) * nx') .* EX + (ones(np, 1) * ny') .* EY) ./ (2*pi*(EX.^2 + EY.^2));
76+v = (h * KE .* (ones(np, 1) * sp')) * sigma;
77+end
78+
79+function v = mfs(n, a, nlobe, par, px, py)
80+% n sources on a circle of radius par, fitted at 3n boundary points.
81+ts = 2*pi*(0:n-1)' / n;
82+sx = par * cos(ts);
83+sy = par * sin(ts);
84+
85+nb = 3 * n;
86+tb = 2*pi*(0:nb-1)' / nb;
87+C = lp_curve(a, nlobe, tb);
88+bx = C(:,1);
89+by = C(:,2);
90+
91+A = source_matrix(bx, by, sx, sy);
92+c = A \ lp_bexact(bx, by);
93+v = source_matrix(px(:), py(:), sx, sy) * c;
94+end
95+
96+function A = source_matrix(x, y, sx, sy)
97+% log|z - z_j| for each source, plus a constant column
98+np = numel(x);
99+n = numel(sx);
100+A = zeros(np, n + 1);
101+for j = 1:n
102+ A(:, j) = 0.5 * log((x - sx(j)).^2 + (y - sy(j)).^2);
103+end
104+A(:, n+1) = 1;
105+end
src/chapters/ch10/chapter.tsadded+123−0View file
@@ -0,0 +1,123 @@
1+import { createElement, Fragment } from 'react'
2+import type { Chapter, Control } from '../../ui/types.ts'
3+import driver from './driver.m?raw'
4+import script from './bie.m?raw'
5+import text from './text.md?raw'
6+import { Fields2D, type Single2DOut } from '../twodim/panels.tsx'
7+import { Convergence, WorkPrecision, fmtErr, type SweepOut } from '../onedim/panels.tsx'
8+import type { MethodKey } from '../../plot/palette.ts'
9+import type { Num } from '../../engine/types.ts'
10+
11+interface BieOut extends SweepOut {
12+ Enear: Num[][]
13+ cx: number[]
14+ cy: number[]
15+}
16+
17+const SLOTS: MethodKey[] = ['bie', 'bie', 'mps', 'mps']
18+const DASHES = [undefined, '5 4', undefined, '5 4']
19+
20+/** Interleave the interior and near-boundary errors into one four-series set. */
21+function split(d: BieOut): SweepOut {
22+ return {
23+ names: [
24+ 'Nyström, interior',
25+ 'Nyström, at 98% radius',
26+ 'fundamental solutions, interior',
27+ 'fundamental solutions, at 98% radius',
28+ ],
29+ ns: d.ns,
30+ E: [d.E[0], d.Enear[0], d.E[1], d.Enear[1]],
31+ T: [d.T[0], d.T[0], d.T[1], d.T[1]],
32+ orders: [d.orders[0], d.orders[0], d.orders[1], d.orders[1]],
33+ }
34+}
35+
36+const GEOM: Control[] = [
37+ { kind: 'slider', key: 'a', label: 'lobe amplitude a', min: 0, max: 0.35, step: 0.05, format: (v) => v.toFixed(2) },
38+ { kind: 'slider', key: 'nlobe', label: 'number of lobes', min: 2, max: 8, step: 1, format: (v) => String(v) },
39+]
40+
41+export const ch10: Chapter = {
42+ id: 'boundary-integral',
43+ nav: 'Boundary integral equations',
44+ part: 'Beyond the square',
45+ title: 'Boundary integral equations',
46+ blurb:
47+ 'A curved domain, an unknown that lives only on the boundary, and full double precision from a few hundred points.',
48+ text,
49+ labs: [
50+ {
51+ id: 'field',
52+ title: 'The solution inside the curve',
53+ caption:
54+ 'The error picture is on a logarithmic colour scale covering whatever range the data spans. Note where it is worst.',
55+ driver,
56+ script,
57+ manual: true,
58+ controls: [
59+ ...GEOM,
60+ {
61+ kind: 'select',
62+ key: 'mid',
63+ label: 'method',
64+ options: [
65+ { value: 1, label: 'double-layer Nyström' },
66+ { value: 2, label: 'fundamental solutions' },
67+ ],
68+ },
69+ { kind: 'slider', key: 'n', label: 'boundary unknowns', values: [16, 32, 64, 128, 256], min: 16, max: 256 },
70+ { kind: 'slider', key: 'par', label: 'MFS source radius', min: 1.6, max: 4, step: 0.2, format: (v) => v.toFixed(1) },
71+ ],
72+ defaults: { a: 0.3, nlobe: 5, mode: 'field', mid: 1, n: 128, par: 2.2 },
73+ render: (d: Single2DOut & { cx: number[]; cy: number[] }) =>
74+ createElement(Fields2D, {
75+ data: d,
76+ xDomain: [d.x[0], d.x[d.x.length - 1]],
77+ yDomain: [d.y[0], d.y[d.y.length - 1]],
78+ curves: [{ x: d.cx, y: d.cy }],
79+ }),
80+ summary: (d: Single2DOut) => [
81+ { label: 'unknowns', value: String(d.unknowns) },
82+ { label: 'max error', value: fmtErr(d.maxerr) },
83+ ],
84+ },
85+ {
86+ id: 'sweep',
87+ title: 'Geometric convergence, and the close-evaluation wall',
88+ caption:
89+ 'Solid curves are errors at interior points; dashed curves the same solution evaluated at 98 per cent of the radius.',
90+ driver,
91+ script,
92+ manual: true,
93+ editorOpen: true,
94+ controls: [
95+ ...GEOM,
96+ { kind: 'slider', key: 'par', label: 'MFS source radius', min: 1.6, max: 4, step: 0.2, format: (v) => v.toFixed(1) },
97+ ],
98+ defaults: { a: 0.3, nlobe: 5, mode: 'sweep', par: 2.2, ns: [32, 64, 128, 256, 384, 512, 768] },
99+ render: (d: BieOut) =>
100+ createElement(
101+ Fragment,
102+ null,
103+ createElement(Convergence, {
104+ key: 'c',
105+ data: split(d),
106+ style: { slots: SLOTS, dashes: DASHES },
107+ height: 340,
108+ xLabel: 'boundary unknowns n',
109+ }),
110+ createElement(WorkPrecision, {
111+ key: 'w',
112+ data: split(d),
113+ style: { slots: SLOTS, dashes: DASHES },
114+ height: 340,
115+ }),
116+ ),
117+ summary: (d: BieOut) => [
118+ { label: 'best Nyström', value: fmtErr(Math.min(...(d.E[0].filter((v) => v != null) as number[]))) },
119+ { label: 'best MFS', value: fmtErr(Math.min(...(d.E[1].filter((v) => v != null) as number[]))) },
120+ ],
121+ },
122+ ],
123+}
src/chapters/ch10/driver.madded+141−0View file
@@ -0,0 +1,141 @@
1+% Driver for the boundary integral chapter.
2+%
3+% The editable script supplies
4+%
5+% names = method_names()
6+% v = solve_bie(mid, n, a, nlobe, par, px, py)
7+%
8+% returning the computed solution at the points (px, py) inside the curve
9+% r(t) = 1 + a cos(nlobe t). n is the method's resolution: boundary nodes for
10+% the Nystrom method, source points for the method of fundamental solutions.
11+% par is the auxiliary parameter each method needs (the radius of the source
12+% circle, for the second one).
13+%
14+% The Dirichlet data is the trace of lp_bexact, so the exact solution is known
15+% everywhere inside and the error needs no reference computation.
16+p = lp_params();
17+a = p.a;
18+nlobe = p.nlobe;
19+names = method_names();
20+nm = numel(names);
21+
22+out = struct();
23+out.names = names;
24+
25+% the curve itself, for the pictures
26+tt = linspace(0, 2*pi, 401)';
27+C = lp_curve(a, nlobe, tt);
28+out.cx = C(:,1)';
29+out.cy = C(:,2)';
30+
31+if strcmp(p.mode, 'field')
32+ % a Cartesian grid over the bounding box, masked to the interior
33+ m = 96;
34+ rmax = 1 + abs(a) + 0.05;
35+ xs = linspace(-rmax, rmax, m)';
36+ [X, Y] = meshgrid(xs, xs);
37+ th = atan2(Y, X);
38+ rr = sqrt(X.^2 + Y.^2);
39+ rb = 1 + a * cos(nlobe * th);
40+ inside = rr < rb - 1e-9;
41+ px = X(inside);
42+ py = Y(inside);
43+
44+ t0 = tic;
45+ v = solve_bie(p.mid, p.n, a, nlobe, p.par, px, py);
46+ out.ms = toc(t0) * 1000;
47+ ve = lp_bexact(px, py);
48+ out.maxerr = max(abs(v - ve));
49+
50+ U = NaN(m, m);
51+ E = NaN(m, m);
52+ U(inside) = v;
53+ E(inside) = abs(v - ve);
54+ out.x = xs';
55+ out.y = xs';
56+ out.U = lp_rows(U);
57+ out.Err = lp_rows(E);
58+ out.unknowns = p.n;
59+else
60+ % errors on a polar lattice: four rings, the last of them close enough to
61+ % the boundary to show what close evaluation costs
62+ nth = 60;
63+ th = 2*pi*(0:nth-1)'/nth;
64+ rb = 1 + a * cos(nlobe * th);
65+ rho = [0.25; 0.5; 0.75];
66+ px = [rho(1)*rb .* cos(th); rho(2)*rb .* cos(th); rho(3)*rb .* cos(th)];
67+ py = [rho(1)*rb .* sin(th); rho(2)*rb .* sin(th); rho(3)*rb .* sin(th)];
68+ ve = lp_bexact(px, py);
69+ % a ring at 98 per cent of the radius, where a layer potential is hardest
70+ % to evaluate accurately
71+ qx = 0.98 * rb .* cos(th);
72+ qy = 0.98 * rb .* sin(th);
73+ qe = lp_bexact(qx, qy);
74+
75+ ns = p.ns(:)';
76+ nn = numel(ns);
77+ E = zeros(nm, nn);
78+ T = zeros(nm, nn);
79+ Q = zeros(nm, nn);
80+ EN = zeros(nm, nn);
81+ ET = drv_methodsb(1, NaN(3*nm, nn), ns, a, nlobe, p.par, px, py, ve, qx, qy, qe, nm);
82+ for k = 1:nm
83+ E(k, :) = ET(3*k-2, :);
84+ EN(k, :) = ET(3*k-1, :);
85+ T(k, :) = ET(3*k, :);
86+ Q(k, :) = lp_orders(E(k, :), ns);
87+ end
88+ out.Enear = lp_rows(EN);
89+ out.ns = ns;
90+ out.E = lp_rows(E);
91+ out.T = lp_rows(T);
92+ out.orders = lp_rows(Q);
93+end
94+
95+lp_write(out);
96+
97+function ET = drv_methodsb(k, ET, ns, a, nlobe, par, px, py, ve, qx, qy, qe, nm)
98+% Recursion rather than a loop, as in the other drivers: a call made from
99+% inside a loop to a routine that itself loops is compiled through numbl's slow
100+% path, and these are timed solves.
101+if k > nm
102+ return;
103+end
104+vw = solve_bie(k, 16, a, nlobe, par, px(1), py(1)); % warm-up
105+R = drv_resb(1, NaN(3, numel(ns)), k, ns, a, nlobe, par, px, py, ve, qx, qy, qe);
106+ET(3*k-2, :) = R(1, :);
107+ET(3*k-1, :) = R(2, :);
108+ET(3*k, :) = R(3, :);
109+ET = drv_methodsb(k+1, ET, ns, a, nlobe, par, px, py, ve, qx, qy, qe, nm);
110+end
111+
112+function R = drv_resb(j, R, mid, ns, a, nlobe, par, px, py, ve, qx, qy, qe)
113+if j > numel(ns)
114+ return;
115+end
116+n = ns(j);
117+t0 = tic;
118+v = solve_bie(mid, n, a, nlobe, par, px, py);
119+el = toc(t0);
120+if el < 0.005
121+ el = drv_bestb(4, el, mid, n, a, nlobe, par, px, py);
122+end
123+q = solve_bie(mid, n, a, nlobe, par, qx, qy);
124+R(1, j) = max(abs(v - ve));
125+R(2, j) = max(abs(q - qe));
126+R(3, j) = el * 1000;
127+R = drv_resb(j+1, R, mid, ns, a, nlobe, par, px, py, ve, qx, qy, qe);
128+end
129+
130+function el = drv_bestb(k, el, mid, n, a, nlobe, par, px, py)
131+if k == 0
132+ return;
133+end
134+t0 = tic;
135+v = solve_bie(mid, n, a, nlobe, par, px, py);
136+e2 = toc(t0);
137+if e2 < el
138+ el = e2;
139+end
140+el = drv_bestb(k-1, el, mid, n, a, nlobe, par, px, py);
141+end
src/chapters/ch10/text.mdadded+60−0View file
@@ -0,0 +1,60 @@
1+Every method up to Chapter 9 discretised the *area*: a grid of unknowns filling the domain, an operator relating each to its neighbours, and an error controlled by how finely the area is covered. This chapter does the opposite. The domain is a smooth closed curve
2+
3+$$r(t) = 1 + a\cos(n_{\text{lobe}}\, t),$$
4+
5+which no tensor-product grid fits, and the unknown lives on the boundary alone.
6+
7+### Reformulating on the boundary
8+
9+Look for the solution as the field of a dipole layer spread over the curve with density $\sigma$:
10+
11+$$u(x) \;=\; \frac{1}{2\pi}\int_\Gamma \frac{n_y \cdot (x-y)}{|x-y|^2}\,\sigma(y)\,ds(y).$$
12+
13+This is the double-layer potential. It is harmonic at every $x$ off the curve, whatever $\sigma$ is, because the integrand is harmonic in $x$; the equation is therefore satisfied *exactly*, and only the boundary condition remains. Letting $x$ approach the curve from inside, the potential jumps by $-\sigma/2$, and matching the data $g$ gives a Fredholm integral equation of the second kind,
14+
15+$$-\tfrac{1}{2}\sigma(x) + \int_\Gamma K(x,y)\,\sigma(y)\,ds(y) \;=\; g(x), \qquad x \in \Gamma.$$
16+
17+Two properties of this equation are what the chapter is about.
18+
19+**It is well conditioned.** The operator is the identity plus a compact operator, so its eigenvalues cluster at $-1/2$ and stay away from zero. The condition number of the discrete system is a modest constant, independent of the discretisation. Contrast the five-point matrix, whose condition number grows like $N^2$ and whose accuracy floor followed from that.
20+
21+**Its kernel is not singular.** $K(x,y)$ looks like it blows up as $y \to x$, and for a general kernel it would, but for the double layer on a smooth curve the numerator vanishes at the same rate as the denominator and the limit is finite:
22+
23+$$\lim_{y \to x} K(x,y) \;=\; -\frac{\kappa(x)}{4\pi},$$
24+
25+with $\kappa$ the curvature. So the integrand is a smooth periodic function of the parameter $t$, and for smooth periodic integrands the plain trapezoid rule is spectrally accurate. Discretising with $n$ equally spaced points in $t$ and the trapezoid rule gives an $n \times n$ dense system whose solution converges geometrically. This is the Nyström method, and on this geometry it is very hard to beat.
26+
27+For comparison the lab also runs the **method of fundamental solutions**: place $n$ point sources on a circle of radius $R$ outside the curve and fit their combined field to the boundary data by least squares. It is three lines of code, every basis function is harmonic, and it also converges geometrically.
28+
29+<!--lab:field-->
30+
31+<!--lab:sweep-->
32+
33+### What the numbers say
34+
35+On the default geometry, measured at interior points out to three quarters of the radius:
36+
37+| method | error | boundary unknowns | time |
38+|---|---|---|---|
39+| double-layer Nyström | $5 \times 10^{-16}$ | 512 | 13 ms |
40+| fundamental solutions | $2 \times 10^{-14}$ | 128 | 12 ms |
41+| five-point on the square (Ch. 6) | $1 \times 10^{-8}$ | $10^{6}$ | 440 ms |
42+
43+The first row is the point of the whole tutorial in one line: full double precision, on a curved domain, from a few hundred unknowns and a few milliseconds, because the formulation put the unknown where the information is and the quadrature matched the smoothness of the integrand.
44+
45+### Three things this does not do
46+
47+**Close evaluation.** The dashed curve in the convergence plot is the error at points at 98 per cent of the radius, and it is catastrophically worse: $10^{-1}$ where the interior error is $10^{-16}$. The reason is that the integrand, smooth for a target well inside, develops a near-singularity as the target approaches the curve, and the trapezoid rule stops being spectrally accurate. Nothing is wrong with the density $\sigma$: it is the *evaluation* that fails. Fixing it needs a different quadrature near the boundary, and the literature on that is substantial. The method of fundamental solutions does not have this problem, because its basis functions are smooth all the way to the curve, and the lab shows its near-boundary error tracking its interior error exactly.
48+
49+**Conditioning of the fitting approach.** The method of fundamental solutions wins on unknowns here, but its least-squares matrix is exponentially ill conditioned in the number of sources, and its accuracy depends on the radius $R$: too small and the sources cannot represent the field, too large and the matrix is numerically rank deficient. Move the slider and watch the error change by orders of magnitude for the same amount of work. There is no comparable knob on the Nyström method. That difference, a method with a free parameter that must be tuned against one without, matters more as problems get harder.
50+
51+**Density and dimension.** The Nyström matrix is dense, so the solve is $O(n^3)$ and the matrix is $O(n^2)$ in memory. In two dimensions with $n$ in the hundreds this is irrelevant. In three dimensions, with $n$ in the hundreds of thousands, it is fatal, and boundary integral methods only became practical there with fast summation: the fast multipole method and its relatives, which apply the operator in $O(n)$ and turn the dense solve into a preconditioned iteration. The well-conditioning of the second-kind equation is what makes that iteration converge in a handful of steps.
52+
53+### Where this leaves the tutorial
54+
55+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.
56+
57+- **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$.
58+- **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$.
59+- **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.
60+- **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.
src/chapters/index.tsadded+13−0View file
@@ -0,0 +1,13 @@
1+import type { Chapter } from '../ui/types.ts'
2+import { ch01 } from './ch01/chapter.ts'
3+import { ch02 } from './ch02/chapter.ts'
4+import { ch03 } from './ch03/chapter.ts'
5+import { ch04 } from './ch04/chapter.ts'
6+import { ch05 } from './ch05/chapter.ts'
7+import { ch06 } from './ch06/chapter.ts'
8+import { ch07 } from './ch07/chapter.ts'
9+import { ch08 } from './ch08/chapter.ts'
10+import { ch09 } from './ch09/chapter.ts'
11+import { ch10 } from './ch10/chapter.ts'
12+
13+export const CHAPTERS: Chapter[] = [ch01, ch02, ch03, ch04, ch05, ch06, ch07, ch08, ch09, ch10]
src/chapters/onedim/driver.madded+121−0View file
@@ -0,0 +1,121 @@
1+% Driver for the one-dimensional labs.
2+%
3+% The editable script supplies two functions:
4+%
5+% names = method_names() a cell array of method labels
6+% [x, u] = solve1d(mid, N, f, ua, ub) method mid, resolution N
7+%
8+% solve1d picks its own grid x (uniform, graded, Chebyshev, whatever the method
9+% wants), returns the solution values on that grid, and the driver measures the
10+% error against the manufactured solution at exactly those points. f is a
11+% function handle giving the source term, so a method is free to evaluate it
12+% anywhere.
13+%
14+% The refinement sweep is written as a recursion rather than a for loop. That
15+% is not a stylistic choice: numbl compiles a call made from inside a loop
16+% through a generic path whenever the callee itself contains loops that call
17+% further functions, and the solves being timed here are exactly of that shape.
18+% The penalty is a factor of five to twenty, it falls unevenly on different
19+% methods, and it would therefore corrupt every timing in this tutorial.
20+% Recursion is compiled normally. See CLAUDE.md.
21+p = lp_params();
22+kind = p.kind;
23+names = method_names();
24+nm = numel(names);
25+ub = lp_exact1d(kind, [0; 1]);
26+f = @(xx) drv_rhs(kind, xx);
27+
28+out = struct();
29+out.names = names;
30+out.kind = kind;
31+
32+if strcmp(p.mode, 'single')
33+ t0 = tic;
34+ [x, u] = solve1d(p.mid, p.N, f, ub(1), ub(2));
35+ out.ms = toc(t0) * 1000;
36+ uex = lp_exact1d(kind, x(:));
37+ err = u(:) - uex(:);
38+ out.x = x(:)';
39+ out.u = u(:)';
40+ out.uex = uex(:)';
41+ out.err = err';
42+ out.maxerr = max(abs(err));
43+ out.n = numel(x) - 1;
44+ xf = linspace(0, 1, p.ngrid)';
45+ out.xf = xf';
46+ out.uf = lp_exact1d(kind, xf)';
47+else
48+ ns = p.ns(:)';
49+ nn = numel(ns);
50+ E = zeros(nm, nn);
51+ T = zeros(nm, nn);
52+ Q = zeros(nm, nn);
53+ nmax = p.nmaxs(:)';
54+ timed = p.timed;
55+ ET = drv_methods(1, NaN(2*nm, nn), ns, f, ub(1), ub(2), kind, timed, nmax, nm);
56+ for k = 1:nm
57+ E(k, :) = ET(2*k-1, :);
58+ T(k, :) = ET(2*k, :);
59+ Q(k, :) = lp_orders(E(k, :), ns);
60+ end
61+ out.ns = ns;
62+ out.E = lp_rows(E);
63+ out.T = lp_rows(T);
64+ out.orders = lp_rows(Q);
65+end
66+
67+lp_write(out);
68+
69+function ET = drv_methods(k, ET, ns, f, ua, ub, kind, timed, nmax, nm)
70+% One method per level of the recursion; row 2k-1 of ET collects its errors and
71+% row 2k its milliseconds.
72+if k > nm
73+ return;
74+end
75+% a warm-up call, so that the first timed solve is not also paying for the
76+% compilation of the method
77+[xw, uw] = solve1d(k, 8, f, ua, ub);
78+R = drv_res(1, NaN(2, numel(ns)), k, ns, f, ua, ub, kind, timed, nmax(k));
79+ET(2*k-1, :) = R(1, :);
80+ET(2*k, :) = R(2, :);
81+ET = drv_methods(k+1, ET, ns, f, ua, ub, kind, timed, nmax, nm);
82+end
83+
84+function R = drv_res(j, R, mid, ns, f, ua, ub, kind, timed, nmax)
85+% One resolution per level of the recursion.
86+if j > numel(ns)
87+ return;
88+end
89+N = ns(j);
90+if N <= nmax
91+ t0 = tic;
92+ [x, u] = solve1d(mid, N, f, ua, ub);
93+ el = toc(t0);
94+ if timed > 0 && el < 0.005
95+ % a sub-millisecond solve says more about the clock than about the
96+ % method, so take the best of a few repeats
97+ el = drv_best(4, el, mid, N, f, ua, ub);
98+ end
99+ uex = lp_exact1d(kind, x(:));
100+ R(1, j) = max(abs(u(:) - uex(:)));
101+ R(2, j) = el * 1000;
102+end
103+R = drv_res(j+1, R, mid, ns, f, ua, ub, kind, timed, nmax);
104+end
105+
106+function el = drv_best(k, el, mid, N, f, ua, ub)
107+if k == 0
108+ return;
109+end
110+t0 = tic;
111+[x, u] = solve1d(mid, N, f, ua, ub);
112+e2 = toc(t0);
113+if e2 < el
114+ el = e2;
115+end
116+el = drv_best(k-1, el, mid, N, f, ua, ub);
117+end
118+
119+function y = drv_rhs(kind, x)
120+[~, y] = lp_exact1d(kind, x);
121+end
src/chapters/onedim/panels.tsxadded+289−0View file
@@ -0,0 +1,289 @@
1+// Figures for the one-dimensional labs. Every lab in Part I sends back one of
2+// these two shapes, so the same two panels serve all of them.
3+import { Fragment } from 'react'
4+import Plot from '../../plot/Plot.tsx'
5+import Legend from '../../plot/Legend.tsx'
6+import { EndLabel, Line, Markers, RefSlope } from '../../plot/Curves.tsx'
7+import { ink, marker, method, type MethodKey } from '../../plot/palette.ts'
8+import { extent, padDomain } from '../../plot/scales.ts'
9+import type { Num } from '../../engine/types.ts'
10+
11+export interface SingleOut {
12+ names: string[]
13+ ms: number
14+ n: number
15+ x: number[]
16+ u: Num[]
17+ uex: number[]
18+ err: Num[]
19+ maxerr: Num
20+ xf: number[]
21+ uf: number[]
22+}
23+
24+export interface SweepOut {
25+ names: string[]
26+ ns: number[]
27+ /** per-method resolutions, when the methods do not share a ladder */
28+ NS?: Num[][]
29+ E: Num[][]
30+ T: Num[][]
31+ orders: Num[][]
32+}
33+
34+/** The x values for series k: its own ladder if it has one. */
35+function xsOf(data: SweepOut, k: number): number[] {
36+ const row = data.NS?.[k]
37+ if (!row) return data.ns
38+ return row.map((v) => (v == null ? NaN : v))
39+}
40+
41+export function fmtErr(v: Num): string {
42+ if (v == null || !isFinite(v)) return '--'
43+ if (v === 0) return '0'
44+ return v.toExponential(2)
45+}
46+
47+/** The computed solution over the exact one, and the error underneath. */
48+export function Solution1D({ data, slot = 'fd2' }: { data: SingleOut; slot?: MethodKey }) {
49+ const color = method[slot]
50+ const yd = padDomain(extent(data.uf, data.u))
51+ const ed = padDomain(extent(data.err))
52+ const showNodes = data.x.length <= 130
53+ return (
54+ <>
55+ <div className="plot-title">the solution: exact, and computed at the grid points</div>
56+ <Plot height={220} xDomain={[0, 1]} yDomain={yd} xLabel="x" yLabel="u">
57+ {(f) => (
58+ <>
59+ <Line x={data.xf} y={data.uf} f={f} color={ink.reference} width={2} />
60+ <Line x={data.x} y={data.u} f={f} color={color} width={1.5} dash="4 3" />
61+ {showNodes && <Markers x={data.x} y={data.u} f={f} color={color} shape={marker[slot]} r={3.2} />}
62+ </>
63+ )}
64+ </Plot>
65+ <Legend
66+ items={[
67+ { label: 'exact solution', color: ink.reference },
68+ { label: `computed (${data.n} intervals)`, color, dash: '4 3' },
69+ ]}
70+ />
71+ <div className="plot-title">the error, computed minus exact, at the grid points</div>
72+ <Plot height={170} xDomain={[0, 1]} yDomain={ed} xLabel="x" yLabel="error">
73+ {(f) => (
74+ <>
75+ <Line x={[0, 1]} y={[0, 0]} f={f} color={ink.axis} width={1} />
76+ <Line x={data.x} y={data.err} f={f} color={color} width={2} />
77+ {showNodes && <Markers x={data.x} y={data.err} f={f} color={color} shape={marker[slot]} r={3.2} />}
78+ </>
79+ )}
80+ </Plot>
81+ </>
82+ )
83+}
84+
85+export interface SweepStyle {
86+ /** palette slot per method, in the order the script lists them */
87+ slots: MethodKey[]
88+ /** dashed reference lines: exponent, label, and an optional anchor point */
89+ refs?: { p: number; label: string; at?: [number, number] }[]
90+ /** per-series dash pattern, when two series share a colour on purpose */
91+ dashes?: (string | undefined)[]
92+}
93+
94+function seriesExtent(rows: Num[][]): [number, number] {
95+ let lo = Infinity
96+ let hi = -Infinity
97+ for (const r of rows) {
98+ for (const v of r) {
99+ if (v == null || !isFinite(v) || v <= 0) continue
100+ if (v < lo) lo = v
101+ if (v > hi) hi = v
102+ }
103+ }
104+ return isFinite(lo) ? [lo / 2, hi * 2] : [1e-16, 1]
105+}
106+
107+/** Max error against resolution, on log-log axes. */
108+export function Convergence({
109+ data,
110+ style,
111+ height = 300,
112+ xLabel = 'N',
113+}: {
114+ data: SweepOut
115+ style: SweepStyle
116+ height?: number
117+ xLabel?: string
118+}) {
119+ const yd = seriesExtent(data.E)
120+ const allx = data.NS ? data.NS.flat().filter((v): v is number => v != null && v > 0) : data.ns
121+ const xd: [number, number] = [Math.min(...allx) / 1.3, Math.max(...allx) * 1.3]
122+ // reference lines are anchored a factor of 5 below the first data point, so
123+ // that the reader compares slopes rather than trying to see through them
124+ const e0 = data.E[0]?.find((v) => v != null && isFinite(v) && v > 0) ?? 1
125+ return (
126+ <>
127+ <div className="plot-title">max error against resolution</div>
128+ <Plot height={height} xDomain={xd} yDomain={yd} xLog yLog xLabel={xLabel} yLabel="max |error|">
129+ {(f) => (
130+ <>
131+ {style.refs?.map((r) => (
132+ <RefSlope
133+ key={r.label}
134+ f={f}
135+ p={r.p}
136+ through={r.at ?? [data.ns[0], (e0 as number) / 5]}
137+ span={1}
138+ label={r.label}
139+ />
140+ ))}
141+ {data.E.map((row, k) => {
142+ const slot = style.slots[k] ?? 'fd2'
143+ return (
144+ <g key={k}>
145+ <Line x={xsOf(data, k)} y={row} f={f} color={method[slot]} dash={style.dashes?.[k]} />
146+ <Markers x={xsOf(data, k)} y={row} f={f} color={method[slot]} shape={marker[slot]} />
147+ </g>
148+ )
149+ })}
150+ </>
151+ )}
152+ </Plot>
153+ <Legend
154+ items={data.names.map((n, k) => ({
155+ label: n,
156+ color: method[style.slots[k] ?? 'fd2'],
157+ dash: style.dashes?.[k],
158+ }))}
159+ />
160+ <OrderTable data={data} />
161+ </>
162+ )
163+}
164+
165+export function OrderTable({ data }: { data: SweepOut }) {
166+ return (
167+ <table className="numtable">
168+ <thead>
169+ <tr>
170+ <th>{data.NS ? 'resolution' : 'N'}</th>
171+ {data.names.map((n) => (
172+ <th key={n} colSpan={2}>
173+ {n}
174+ </th>
175+ ))}
176+ </tr>
177+ <tr>
178+ <th />
179+ {data.names.map((n) => (
180+ <Fragment key={n}>
181+ <th>error</th>
182+ <th>order</th>
183+ </Fragment>
184+ ))}
185+ </tr>
186+ </thead>
187+ <tbody>
188+ {data.ns.map((n, j) => (
189+ <tr key={n}>
190+ <td>{data.NS ? data.NS.map((r) => r[j] ?? '--').join(' / ') : n}</td>
191+ {data.names.map((nm, k) => (
192+ <Fragment key={nm}>
193+ <td>{fmtErr(data.E[k][j])}</td>
194+ <td style={{ color: ink.muted }}>
195+ {data.orders[k][j] == null ? '' : (data.orders[k][j] as number).toFixed(2)}
196+ </td>
197+ </Fragment>
198+ ))}
199+ </tr>
200+ ))}
201+ </tbody>
202+ </table>
203+ )
204+}
205+
206+/** Error against wall clock: the only comparison that decides anything. */
207+export function WorkPrecision({
208+ data,
209+ style,
210+ height = 320,
211+ label = 'milliseconds for one solve',
212+}: {
213+ data: SweepOut
214+ style: SweepStyle
215+ height?: number
216+ label?: string
217+}) {
218+ const yd = seriesExtent(data.E)
219+ const xd = seriesExtent(data.T)
220+ return (
221+ <>
222+ <div className="plot-title">accuracy against cost; down and to the left is better</div>
223+ <Plot
224+ height={height}
225+ xDomain={xd}
226+ yDomain={yd}
227+ xLog
228+ yLog
229+ xLabel={label}
230+ yLabel="max |error|"
231+ 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
243+ }
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}
251+ 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}
259+ />
260+ )
261+ })}
262+ </>
263+ )}
264+ >
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+ )}
279+ </Plot>
280+ <Legend
281+ items={data.names.map((n, k) => ({
282+ label: n,
283+ color: method[style.slots[k] ?? 'fd2'],
284+ dash: style.dashes?.[k],
285+ }))}
286+ />
287+ </>
288+ )
289+}
src/chapters/shared.tsadded+32−0View file
@@ -0,0 +1,32 @@
1+// Controls and resolution ladders reused across chapters, so that "the smooth
2+// case" means the same thing everywhere and two chapters can be compared.
3+import type { Control } from '../ui/types.ts'
4+
5+export const KIND_CONTROL: Control = {
6+ kind: 'select',
7+ key: 'kind',
8+ label: 'manufactured solution',
9+ options: [
10+ { value: 'smooth', label: 'smooth: exp(sin 3πx)' },
11+ { value: 'oscil', label: 'oscillatory: sin 8πx + x' },
12+ { value: 'peak', label: 'narrow peak at x = 0.4' },
13+ { value: 'layer', label: 'boundary layer, width 1/50' },
14+ { value: 'lowreg', label: 'low regularity: |x − ½|^{5/2}' },
15+ ],
16+}
17+
18+export const EXACT2D_CONTROL: Control = {
19+ kind: 'select',
20+ key: 'kind',
21+ label: 'harmonic solution',
22+ options: [
23+ { value: 'expsin', label: 'exp(x) sin(y)' },
24+ { value: 'poly', label: 'Re (x + iy)^5' },
25+ { value: 'wave', label: 'exp(ky) sin(kx), k = 8' },
26+ { value: 'log', label: 'log|z − z₀|, z₀ just outside' },
27+ ],
28+}
29+
30+export const NS_1D = [8, 16, 32, 64, 128, 256, 512, 1024]
31+export const NS_1D_LONG = [16, 64, 256, 1024, 4096, 16384, 65536, 262144, 1048576]
32+export const NS_2D = [8, 16, 32, 64, 128, 256]
src/chapters/twodim/driver.madded+107−0View file
@@ -0,0 +1,107 @@
1+% Driver for the two-dimensional labs on the unit square.
2+%
3+% The editable script supplies
4+%
5+% names = method_names()
6+% [x, y, U] = solve2d(mid, N, g)
7+%
8+% where g(x, y) is the Dirichlet data. Because every exact solution here is
9+% harmonic on the whole plane, g happens to be defined everywhere, but a method
10+% is only entitled to evaluate it on the boundary of the square; anything else
11+% is looking up the answer.
12+%
13+% U is the solution on the tensor grid x by y, with U(i, j) at (x(j), y(i)), and
14+% the driver measures the error against the exact solution on exactly that grid.
15+p = lp_params();
16+kind = p.kind;
17+par = p.par;
18+names = method_names();
19+nm = numel(names);
20+g = @(xx, yy) lp_exact2d(kind, xx, yy, par);
21+
22+out = struct();
23+out.names = names;
24+
25+if strcmp(p.mode, 'single')
26+ t0 = tic;
27+ [x, y, U] = solve2d(p.mid, p.N, g);
28+ out.ms = toc(t0) * 1000;
29+ x = x(:);
30+ y = y(:);
31+ [X, Y] = meshgrid(x, y);
32+ UE = lp_exact2d(kind, X, Y, par);
33+ Err = U - UE;
34+ out.maxerr = max(max(abs(Err)));
35+ out.unknowns = (numel(x) - 2) * (numel(y) - 2);
36+ [Us, xs, ys] = lp_thin(U, x, y, p.maxpix);
37+ Es = lp_thin(abs(Err), x, y, p.maxpix);
38+ out.x = xs(:)';
39+ out.y = ys(:)';
40+ out.U = lp_rows(Us);
41+ out.Err = lp_rows(Es);
42+else
43+ ns = p.ns(:)';
44+ nn = numel(ns);
45+ E = zeros(nm, nn);
46+ T = zeros(nm, nn);
47+ Q = zeros(nm, nn);
48+ nmax = p.nmaxs(:)';
49+ timed = p.timed;
50+ ET = drv_methods2(1, NaN(2*nm, nn), ns, g, kind, par, nmax, timed, nm);
51+ for k = 1:nm
52+ E(k, :) = ET(2*k-1, :);
53+ T(k, :) = ET(2*k, :);
54+ Q(k, :) = lp_orders(E(k, :), ns);
55+ end
56+ out.ns = ns;
57+ out.E = lp_rows(E);
58+ out.T = lp_rows(T);
59+ out.orders = lp_rows(Q);
60+end
61+
62+lp_write(out);
63+
64+function ET = drv_methods2(k, ET, ns, g, kind, par, nmax, timed, nm)
65+% One method per level of the recursion. See the note in the one-dimensional
66+% driver: a for loop here would put every timed solve on numbl's slow path.
67+if k > nm
68+ return;
69+end
70+[xw, yw, Uw] = solve2d(k, 8, g); % warm-up, so the first timing is honest
71+R = drv_res2(1, NaN(2, numel(ns)), k, ns, g, kind, par, nmax(k), timed);
72+ET(2*k-1, :) = R(1, :);
73+ET(2*k, :) = R(2, :);
74+ET = drv_methods2(k+1, ET, ns, g, kind, par, nmax, timed, nm);
75+end
76+
77+function R = drv_res2(j, R, mid, ns, g, kind, par, nmax, timed)
78+if j > numel(ns)
79+ return;
80+end
81+N = ns(j);
82+if N <= nmax
83+ t0 = tic;
84+ [x, y, U] = solve2d(mid, N, g);
85+ el = toc(t0);
86+ if timed > 0 && el < 0.005
87+ el = drv_best2(4, el, mid, N, g);
88+ end
89+ [X, Y] = meshgrid(x(:), y(:));
90+ R(1, j) = max(max(abs(U - lp_exact2d(kind, X, Y, par))));
91+ R(2, j) = el * 1000;
92+end
93+R = drv_res2(j+1, R, mid, ns, g, kind, par, nmax, timed);
94+end
95+
96+function el = drv_best2(k, el, mid, N, g)
97+if k == 0
98+ return;
99+end
100+t0 = tic;
101+[x, y, U] = solve2d(mid, N, g);
102+e2 = toc(t0);
103+if e2 < el
104+ el = e2;
105+end
106+el = drv_best2(k-1, el, mid, N, g);
107+end
src/chapters/twodim/panels.tsxadded+83−0View file
@@ -0,0 +1,83 @@
1+// Figures for the two-dimensional labs: the computed field, and the pointwise
2+// error beside it on a logarithmic colour scale, which is the only way to see
3+// both the interior and the corners of an error that spans decades.
4+import Field from '../../plot/Field.tsx'
5+import type { Num } from '../../engine/types.ts'
6+import { fmtErr } from '../onedim/panels.tsx'
7+
8+export interface Single2DOut {
9+ names: string[]
10+ ms: number
11+ maxerr: Num
12+ unknowns: number
13+ x: number[]
14+ y: number[]
15+ U: Num[][]
16+ Err: Num[][]
17+}
18+
19+/** A signed field gets the diverging ramp; a one-signed field the sequential
20+ * one, so that the whole colour range is used either way. */
21+function signOf(z: Num[][]): 'div' | 'seq' {
22+ let lo = Infinity
23+ let hi = -Infinity
24+ for (const row of z) {
25+ for (const v of row) {
26+ if (v == null || !isFinite(v)) continue
27+ if (v < lo) lo = v
28+ if (v > hi) hi = v
29+ }
30+ }
31+ return lo < 0 && hi > 0 ? 'div' : 'seq'
32+}
33+
34+/** Decades to span in an error picture: what the data actually covers. */
35+function decadesOf(z: Num[][]): number {
36+ let lo = Infinity
37+ let hi = 0
38+ for (const row of z) {
39+ for (const v of row) {
40+ if (v == null || !isFinite(v) || v <= 0) continue
41+ if (v < lo) lo = v
42+ if (v > hi) hi = v
43+ }
44+ }
45+ if (!isFinite(lo) || hi <= 0) return 6
46+ return Math.min(14, Math.max(3, Math.ceil(Math.log10(hi / lo))))
47+}
48+
49+export function Fields2D({
50+ data,
51+ xDomain = [0, 1],
52+ yDomain = [0, 1],
53+ curves,
54+}: {
55+ data: Single2DOut
56+ xDomain?: [number, number]
57+ yDomain?: [number, number]
58+ curves?: { x: number[]; y: number[]; color?: string }[]
59+}) {
60+ const emax = data.maxerr == null || !isFinite(data.maxerr) ? 1 : Math.log10(data.maxerr)
61+ return (
62+ <div className="figure-row">
63+ <Field
64+ z={data.U}
65+ xDomain={xDomain}
66+ yDomain={yDomain}
67+ kind={signOf(data.U)}
68+ title="computed solution"
69+ curves={curves}
70+ />
71+ <Field
72+ z={data.Err}
73+ xDomain={xDomain}
74+ yDomain={yDomain}
75+ kind="log"
76+ zmax={Math.ceil(emax)}
77+ decades={decadesOf(data.Err)}
78+ title={`|error|, largest ${fmtErr(data.maxerr)}`}
79+ curves={curves}
80+ />
81+ </div>
82+ )
83+}
src/editor/ScriptEditor.tsxadded+70−0View file
@@ -0,0 +1,70 @@
1+import { useEffect, useRef } from 'react'
2+import { EditorView, basicSetup } from 'codemirror'
3+import { EditorState } from '@codemirror/state'
4+import { keymap } from '@codemirror/view'
5+import { indentUnit, StreamLanguage } from '@codemirror/language'
6+import { octave } from '@codemirror/legacy-modes/mode/octave'
7+import { oneDark } from '@codemirror/theme-one-dark'
8+
9+export interface ScriptEditorProps {
10+ value: string
11+ onChange: (value: string) => void
12+ /** Ctrl/Cmd+Enter */
13+ onRun: () => void
14+}
15+
16+export default function ScriptEditor({ value, onChange, onRun }: ScriptEditorProps) {
17+ const hostRef = useRef<HTMLDivElement>(null)
18+ const viewRef = useRef<EditorView | null>(null)
19+ const onChangeRef = useRef(onChange)
20+ onChangeRef.current = onChange
21+ const onRunRef = useRef(onRun)
22+ onRunRef.current = onRun
23+
24+ useEffect(() => {
25+ if (!hostRef.current) return
26+ const view = new EditorView({
27+ parent: hostRef.current,
28+ state: EditorState.create({
29+ doc: value,
30+ extensions: [
31+ basicSetup,
32+ keymap.of([
33+ {
34+ key: 'Mod-Enter',
35+ run: () => {
36+ onRunRef.current()
37+ return true
38+ },
39+ },
40+ ]),
41+ StreamLanguage.define(octave),
42+ oneDark,
43+ indentUnit.of(' '),
44+ EditorView.updateListener.of((update) => {
45+ if (update.docChanged) onChangeRef.current(update.state.doc.toString())
46+ }),
47+ ],
48+ }),
49+ })
50+ viewRef.current = view
51+ return () => {
52+ view.destroy()
53+ viewRef.current = null
54+ }
55+ // created once; external value changes are synced below
56+ // eslint-disable-next-line react-hooks/exhaustive-deps
57+ }, [])
58+
59+ // sync an externally loaded script (a different method picked) into the editor
60+ useEffect(() => {
61+ const view = viewRef.current
62+ if (!view) return
63+ const current = view.state.doc.toString()
64+ if (current !== value) {
65+ view.dispatch({ changes: { from: 0, to: current.length, insert: value } })
66+ }
67+ }, [value])
68+
69+ return <div className="editor-host" ref={hostRef} />
70+}
src/engine/lib.tsadded+13−0View file
@@ -0,0 +1,13 @@
1+// Shared .m function files, put into every session as ordinary function files.
2+import type { BootFile } from 'numbl/browser'
3+
4+const libRaw = import.meta.glob('../matlab/*.m', {
5+ query: '?raw',
6+ import: 'default',
7+ eager: true,
8+}) as Record<string, string>
9+
10+export const LIB_FILES: BootFile[] = Object.entries(libRaw).map(([path, content]) => ({
11+ path: path.slice(path.lastIndexOf('/') + 1),
12+ content,
13+}))
src/engine/runner.tsadded+0−0View file
Binary file not shown.
src/engine/types.tsadded+15−0View file
@@ -0,0 +1,15 @@
1+// The contract between the app and the MATLAB layer.
2+//
3+// Every lab writes its controls into params.json and reads back out.json, so
4+// the TypeScript side never needs to know more than "some JSON goes in, some
5+// JSON comes out". Individual panels narrow the output to their own shape.
6+
7+/** NaN and +-Inf all cross the JSON boundary as null. */
8+export type Num = number | null
9+
10+export type ParamValue = number | string | boolean | number[] | number[][]
11+export type Params = Record<string, ParamValue>
12+
13+export type RunResult<T> =
14+ | { ok: true; data: T; output: string; ms: number }
15+ | { ok: false; error: string; output: string; ms: number }
src/index.cssadded+527−0View file
@@ -0,0 +1,527 @@
1+:root {
2+ --surface: #1a1a19;
3+ --surface-2: #201f1e;
4+ --surface-3: #262523;
5+ --line: #3c3c38;
6+ --line-soft: #2b2a28;
7+ --text: #e8e6df;
8+ --text-2: #c3c2b7;
9+ --muted: #8a8880;
10+ --accent: #3987e5;
11+ --warn: #c98500;
12+ --bad: #e66767;
13+ color-scheme: dark;
14+}
15+
16+* {
17+ box-sizing: border-box;
18+}
19+
20+html,
21+body {
22+ margin: 0;
23+ padding: 0;
24+ background: var(--surface);
25+ color: var(--text);
26+ font-family: ui-sans-serif, system-ui, -apple-system, 'Segoe UI', Roboto, sans-serif;
27+ font-size: 15px;
28+ line-height: 1.6;
29+}
30+
31+a {
32+ color: var(--accent);
33+}
34+
35+code,
36+pre,
37+.mono {
38+ font-family: ui-monospace, SFMono-Regular, 'SF Mono', Menlo, Consolas, monospace;
39+}
40+
41+.app {
42+ display: flex;
43+ align-items: flex-start;
44+ min-height: 100vh;
45+}
46+
47+/* ---------- sidebar ---------- */
48+
49+.sidebar {
50+ position: sticky;
51+ top: 0;
52+ flex: 0 0 250px;
53+ height: 100vh;
54+ overflow-y: auto;
55+ padding: 20px 14px 24px 20px;
56+ border-right: 1px solid var(--line-soft);
57+ background: var(--surface-2);
58+}
59+
60+.brand-title {
61+ font-size: 16px;
62+ font-weight: 600;
63+ letter-spacing: 0.01em;
64+}
65+
66+.brand-sub {
67+ font-size: 12px;
68+ color: var(--muted);
69+ margin-bottom: 18px;
70+}
71+
72+.nav-part-title {
73+ font-size: 11px;
74+ text-transform: uppercase;
75+ letter-spacing: 0.08em;
76+ color: var(--muted);
77+ margin: 16px 0 4px;
78+}
79+
80+.sidebar ol {
81+ list-style: none;
82+ margin: 0;
83+ padding: 0;
84+}
85+
86+.nav-item {
87+ display: flex;
88+ gap: 8px;
89+ align-items: baseline;
90+ width: 100%;
91+ text-align: left;
92+ background: none;
93+ border: none;
94+ color: var(--text-2);
95+ padding: 4px 6px;
96+ border-radius: 5px;
97+ font: inherit;
98+ font-size: 13.5px;
99+ line-height: 1.35;
100+ cursor: pointer;
101+}
102+
103+.nav-item:hover {
104+ background: var(--surface-3);
105+ color: var(--text);
106+}
107+
108+.nav-item.active {
109+ background: var(--surface-3);
110+ color: var(--text);
111+ box-shadow: inset 2px 0 0 var(--accent);
112+}
113+
114+.nav-num {
115+ display: inline-block;
116+ min-width: 14px;
117+ color: var(--muted);
118+ font-variant-numeric: tabular-nums;
119+ font-size: 12px;
120+}
121+
122+.nav-foot {
123+ margin-top: 26px;
124+ font-size: 11.5px;
125+ line-height: 1.5;
126+ color: var(--muted);
127+}
128+
129+/* ---------- main column ---------- */
130+
131+.main {
132+ flex: 1 1 auto;
133+ min-width: 0;
134+ padding: 28px 32px 80px;
135+}
136+
137+.chapter {
138+ max-width: 860px;
139+ margin: 0 auto;
140+}
141+
142+.chapter-part {
143+ font-size: 11px;
144+ text-transform: uppercase;
145+ letter-spacing: 0.08em;
146+ color: var(--muted);
147+}
148+
149+.chapter-head h2 {
150+ margin: 4px 0 6px;
151+ font-size: 25px;
152+ font-weight: 600;
153+ letter-spacing: -0.01em;
154+}
155+
156+.chapter-blurb {
157+ margin: 0 0 8px;
158+ color: var(--text-2);
159+}
160+
161+.prose {
162+ max-width: 720px;
163+}
164+
165+.prose h3 {
166+ font-size: 17px;
167+ margin: 26px 0 6px;
168+}
169+
170+.prose p {
171+ margin: 12px 0;
172+}
173+
174+.prose ul,
175+.prose ol {
176+ padding-left: 22px;
177+}
178+
179+.prose li {
180+ margin: 5px 0;
181+}
182+
183+.prose code {
184+ background: var(--surface-3);
185+ padding: 1px 4px;
186+ border-radius: 3px;
187+ font-size: 13px;
188+}
189+
190+.prose blockquote {
191+ margin: 14px 0;
192+ padding: 2px 14px;
193+ border-left: 2px solid var(--line);
194+ color: var(--text-2);
195+}
196+
197+.prose table {
198+ border-collapse: collapse;
199+ font-size: 13.5px;
200+ margin: 12px 0;
201+}
202+
203+.prose th,
204+.prose td {
205+ border-bottom: 1px solid var(--line-soft);
206+ padding: 4px 12px 4px 0;
207+ text-align: left;
208+ font-variant-numeric: tabular-nums;
209+}
210+
211+.katex-display {
212+ overflow-x: auto;
213+ overflow-y: hidden;
214+ padding: 2px 0;
215+}
216+
217+/* ---------- labs ---------- */
218+
219+.lab {
220+ margin: 24px 0 30px;
221+ padding: 14px 16px 16px;
222+ border: 1px solid var(--line-soft);
223+ border-radius: 8px;
224+ background: var(--surface-2);
225+}
226+
227+.lab-head h3 {
228+ margin: 0;
229+ font-size: 15px;
230+ font-weight: 600;
231+}
232+
233+.lab-caption {
234+ margin: 2px 0 0;
235+ font-size: 13px;
236+ color: var(--muted);
237+}
238+
239+.controls {
240+ display: flex;
241+ flex-wrap: wrap;
242+ gap: 10px 20px;
243+ margin: 12px 0 4px;
244+}
245+
246+.control {
247+ display: flex;
248+ flex-direction: column;
249+ gap: 2px;
250+ font-size: 12px;
251+}
252+
253+.control-toggle {
254+ flex-direction: row;
255+ align-items: center;
256+ gap: 6px;
257+}
258+
259+.control-label {
260+ color: var(--muted);
261+}
262+
263+.control-value {
264+ color: var(--text);
265+ font-variant-numeric: tabular-nums;
266+}
267+
268+.control input[type='range'] {
269+ width: 150px;
270+ accent-color: var(--accent);
271+}
272+
273+.control select {
274+ background: var(--surface-3);
275+ color: var(--text);
276+ border: 1px solid var(--line);
277+ border-radius: 4px;
278+ padding: 2px 6px;
279+ font: inherit;
280+ font-size: 12.5px;
281+}
282+
283+.lab-bar {
284+ display: flex;
285+ align-items: center;
286+ gap: 14px;
287+ margin: 10px 0 6px;
288+ flex-wrap: wrap;
289+}
290+
291+.lab-bar-gap {
292+ flex: 1 1 auto;
293+}
294+
295+button {
296+ font: inherit;
297+ font-size: 12.5px;
298+ padding: 3px 12px;
299+ border-radius: 5px;
300+ border: 1px solid var(--line);
301+ background: var(--surface-3);
302+ color: var(--text);
303+ cursor: pointer;
304+}
305+
306+button:hover:not(:disabled) {
307+ border-color: var(--muted);
308+}
309+
310+button:disabled {
311+ opacity: 0.5;
312+ cursor: default;
313+}
314+
315+button.primary {
316+ background: #1c5cab;
317+ border-color: #256abf;
318+}
319+
320+button.ghost {
321+ background: none;
322+ border-color: transparent;
323+ color: var(--muted);
324+}
325+
326+button.ghost:hover {
327+ color: var(--text);
328+}
329+
330+.stat {
331+ display: flex;
332+ flex-direction: column;
333+ line-height: 1.25;
334+}
335+
336+.stat-label {
337+ font-size: 10.5px;
338+ text-transform: uppercase;
339+ letter-spacing: 0.06em;
340+ color: var(--muted);
341+}
342+
343+.stat-value {
344+ font-size: 13px;
345+ font-variant-numeric: tabular-nums;
346+ font-family: ui-monospace, SFMono-Regular, Menlo, Consolas, monospace;
347+}
348+
349+.warn {
350+ color: var(--warn);
351+ font-size: 12px;
352+}
353+
354+.lab-error {
355+ color: var(--bad);
356+ background: #2a1d1d;
357+ border: 1px solid #4a2b2b;
358+ border-radius: 6px;
359+ padding: 8px 10px;
360+ font-size: 12.5px;
361+ white-space: pre-wrap;
362+ overflow-x: auto;
363+}
364+
365+.lab-log {
366+ background: var(--surface-3);
367+ border-radius: 6px;
368+ padding: 8px 10px;
369+ font-size: 12px;
370+ max-height: 260px;
371+ overflow: auto;
372+ white-space: pre;
373+ color: var(--text-2);
374+}
375+
376+.lab-placeholder {
377+ color: var(--muted);
378+ font-size: 13px;
379+ padding: 26px 0;
380+ text-align: center;
381+}
382+
383+.lab-figures {
384+ display: flex;
385+ flex-direction: column;
386+ gap: 6px;
387+}
388+
389+.figure-row {
390+ display: flex;
391+ flex-wrap: wrap;
392+ gap: 10px;
393+}
394+
395+.figure-row > * {
396+ flex: 1 1 320px;
397+ min-width: 300px;
398+}
399+
400+.lab-editor {
401+ margin-top: 10px;
402+ border: 1px solid var(--line-soft);
403+ border-radius: 6px;
404+ overflow: hidden;
405+}
406+
407+.lab-editor-bar {
408+ display: flex;
409+ align-items: center;
410+ gap: 10px;
411+ padding: 4px 8px;
412+ font-size: 12px;
413+ color: var(--muted);
414+ background: var(--surface-3);
415+}
416+
417+.editor-host .cm-editor {
418+ max-height: 460px;
419+ font-size: 12.5px;
420+}
421+
422+/* ---------- plots ---------- */
423+
424+.plot-host {
425+ width: 100%;
426+}
427+
428+.plot-title {
429+ font-size: 12px;
430+ color: var(--text-2);
431+ margin: 6px 0 -4px;
432+}
433+
434+.legend {
435+ display: flex;
436+ flex-wrap: wrap;
437+ gap: 4px 16px;
438+ font-size: 12px;
439+ color: var(--text-2);
440+ margin: -4px 0 6px 44px;
441+}
442+
443+.legend-item {
444+ display: inline-flex;
445+ align-items: center;
446+ gap: 5px;
447+}
448+
449+.legend-value {
450+ color: var(--muted);
451+ font-variant-numeric: tabular-nums;
452+}
453+
454+/* ---------- tables of numbers ---------- */
455+
456+.numtable {
457+ border-collapse: collapse;
458+ font-size: 12.5px;
459+ font-variant-numeric: tabular-nums;
460+ font-family: ui-monospace, SFMono-Regular, Menlo, Consolas, monospace;
461+ margin: 4px 0 2px;
462+}
463+
464+.numtable th {
465+ text-align: right;
466+ color: var(--muted);
467+ font-weight: 500;
468+ padding: 2px 10px 2px 0;
469+ border-bottom: 1px solid var(--line-soft);
470+}
471+
472+.numtable td {
473+ text-align: right;
474+ padding: 1px 10px 1px 0;
475+}
476+
477+/* ---------- library ---------- */
478+
479+.lib-file {
480+ border-bottom: 1px solid var(--line-soft);
481+}
482+
483+.lib-name {
484+ background: none;
485+ border: none;
486+ color: var(--text);
487+ font-family: ui-monospace, SFMono-Regular, Menlo, Consolas, monospace;
488+ font-size: 13px;
489+ padding: 6px 0;
490+ cursor: pointer;
491+}
492+
493+.lib-src {
494+ background: var(--surface-2);
495+ border-radius: 6px;
496+ padding: 10px 12px;
497+ font-size: 12.5px;
498+ overflow-x: auto;
499+ color: var(--text-2);
500+}
501+
502+.pager {
503+ display: flex;
504+ max-width: 860px;
505+ margin: 30px auto 0;
506+ padding-top: 14px;
507+ border-top: 1px solid var(--line-soft);
508+}
509+
510+@media (max-width: 900px) {
511+ .app {
512+ flex-direction: column;
513+ }
514+
515+ .sidebar {
516+ position: static;
517+ width: 100%;
518+ height: auto;
519+ flex: none;
520+ border-right: none;
521+ border-bottom: 1px solid var(--line-soft);
522+ }
523+
524+ .main {
525+ padding: 18px 14px 60px;
526+ }
527+}
src/main.tsxadded+10−0View file
@@ -0,0 +1,10 @@
1+import { StrictMode } from 'react'
2+import { createRoot } from 'react-dom/client'
3+import App from './App.tsx'
4+import './index.css'
5+
6+createRoot(document.getElementById('root')!).render(
7+ <StrictMode>
8+ <App />
9+ </StrictMode>,
10+)
src/matlab/lp_bexact.madded+17−0View file
@@ -0,0 +1,17 @@
1+function U = lp_bexact(X, Y)
2+% LP_BEXACT Ground truth for the interior Dirichlet problem on a curve.
3+%
4+% A sum of point sources placed outside the curve,
5+%
6+% u(z) = sum_k c_k log|z - z_k| ,
7+%
8+% which is harmonic everywhere except at the sources and therefore harmonic on
9+% the whole of the interior. Its trace on the curve is the Dirichlet data, and
10+% it is the exact solution by construction, to full precision at every point.
11+zk = [1.8 + 0.9i, -1.5 - 1.7i, 0.2 + 2.4i];
12+ck = [1.0, -0.7, 0.5];
13+U = zeros(size(X));
14+for k = 1:numel(zk)
15+ U = U + ck(k) * 0.5 * log((X - real(zk(k))).^2 + (Y - imag(zk(k))).^2);
16+end
17+end
src/matlab/lp_cheb.madded+13−0View file
@@ -0,0 +1,13 @@
1+function [D, x] = lp_cheb(N)
2+% LP_CHEB Chebyshev differentiation matrix and nodes on [-1, 1].
3+%
4+% x = cos(pi*(0:N)/N), ordered from 1 down to -1, and D is the (N+1)-by-(N+1)
5+% matrix with (D*v)(i) = p'(x_i) for the polynomial p of degree <= N through
6+% the data v. This is Trefethen's cheb.m (Spectral Methods in MATLAB, 2000).
7+x = cos(pi*(0:N)'/N);
8+c = [2; ones(N-1, 1); 2] .* (-1).^(0:N)';
9+X = repmat(x, 1, N+1);
10+dX = X - X';
11+D = (c * (1 ./ c)') ./ (dX + eye(N+1)); % off-diagonal entries
12+D = D - diag(sum(D, 2)); % diagonal, by the zero-row-sum rule
13+end
src/matlab/lp_curve.madded+25−0View file
@@ -0,0 +1,25 @@
1+function C = lp_curve(a, nlobe, t)
2+% LP_CURVE A smooth closed curve and its first two derivatives.
3+%
4+% r(t) = 1 + a cos(nlobe t), x = r cos t, y = r sin t
5+%
6+% Returned as one matrix with columns [x y x' y' x'' y''] rather than as six
7+% outputs, because the drivers call this from inside timing loops and numbl
8+% compiles a multiple-output call in a loop through its slow path.
9+%
10+% The curve is analytic and traversed counterclockwise, which is what makes the
11+% trapezoid rule spectrally accurate on it and fixes the sign of the outward
12+% normal (y', -x')/|(x', y')|.
13+t = t(:);
14+r = 1 + a * cos(nlobe * t);
15+rp = -a * nlobe * sin(nlobe * t);
16+rpp = -a * nlobe^2 * cos(nlobe * t);
17+c = cos(t);
18+s = sin(t);
19+C = [r .* c, ...
20+ r .* s, ...
21+ rp .* c - r .* s, ...
22+ rp .* s + r .* c, ...
23+ rpp .* c - 2*rp .* s - r .* c, ...
24+ rpp .* s + 2*rp .* c - r .* s];
25+end
src/matlab/lp_dst1.madded+14−0View file
@@ -0,0 +1,14 @@
1+function Y = lp_dst1(X)
2+% LP_DST1 Discrete sine transform of type I along the columns of X.
3+%
4+% Y(k,:) = sqrt(2/(n+1)) * sum_j X(j,:) sin(pi*j*k/(n+1)). The scaling makes
5+% the transform orthogonal, so lp_dst1(lp_dst1(X)) is X, and it is computed
6+% with one complex FFT of the odd extension of X, of length 2(n+1).
7+n = size(X, 1);
8+m = size(X, 2);
9+Z = zeros(2*(n+1), m);
10+Z(2:n+1, :) = X;
11+Z(n+3:2*(n+1), :) = -flipud(X);
12+W = fft(Z);
13+Y = -imag(W(2:n+1, :)) * sqrt(1/(2*(n+1)));
14+end
src/matlab/lp_exact1d.madded+43−0View file
@@ -0,0 +1,43 @@
1+function [u, f] = lp_exact1d(kind, x)
2+% LP_EXACT1D A manufactured solution of -u'' = f on [0, 1], and its source.
3+%
4+% Pick u first, differentiate it by hand, and call the result f. Then u is by
5+% construction the exact solution of the boundary value problem
6+%
7+% -u''(x) = f(x) on (0, 1), u(0), u(1) given by u itself,
8+%
9+% and any difference between a computed solution and u is error, known to full
10+% precision at every point and at no cost. The five choices below differ only
11+% in how hard they are to resolve.
12+switch kind
13+ case 'smooth'
14+ % entire, moderate amplitude: the easy case
15+ s = sin(3*pi*x);
16+ u = exp(s);
17+ f = (3*pi)^2 * (s - cos(3*pi*x).^2) .* u;
18+ case 'oscil'
19+ % eight wavelengths across the interval: still entire, but a mesh must
20+ % resolve the wavelength before any convergence rate appears
21+ u = sin(8*pi*x) + x;
22+ f = (8*pi)^2 * sin(8*pi*x);
23+ case 'peak'
24+ % a Gaussian of width 1/20 at x = 0.4: entire, but with all of its
25+ % structure in a tenth of the interval
26+ a = 200;
27+ u = exp(-a*(x - 0.4).^2);
28+ f = (2*a - 4*a^2*(x - 0.4).^2) .* u;
29+ case 'layer'
30+ % a boundary layer of width eps at x = 1, as in a singular perturbation
31+ ep = 0.02;
32+ u = exp(-(1 - x)/ep);
33+ f = -u/ep^2;
34+ case 'lowreg'
35+ % u is C^2 but not C^3: |x - 1/2|^(5/2) has a square root in u''.
36+ % Every polynomial-based method feels this.
37+ c = 0.5;
38+ u = abs(x - c).^(5/2);
39+ f = -(15/4) * sqrt(abs(x - c));
40+ otherwise
41+ error('lp_exact1d: unknown solution "%s"', kind);
42+end
43+end
src/matlab/lp_exact2d.madded+33−0View file
@@ -0,0 +1,33 @@
1+function U = lp_exact2d(kind, X, Y, par)
2+% LP_EXACT2D A harmonic function on the unit square, used as ground truth.
3+%
4+% Every choice here satisfies Laplace's equation exactly, so taking its trace
5+% on the boundary as Dirichlet data makes it the exact solution of the
6+% Dirichlet problem, and the difference between a computed solution and this
7+% function is the whole of the error. par is a per-kind shape parameter.
8+if nargin < 4
9+ par = 1;
10+end
11+switch kind
12+ case 'expsin'
13+ % exp(x) sin(y): entire, gentle, the standard test
14+ U = exp(X) .* sin(Y);
15+ case 'poly'
16+ % Re (x + iy)^k: a harmonic polynomial of degree k
17+ k = round(par);
18+ U = real((X + 1i*Y).^k);
19+ case 'wave'
20+ % exp(k*y) sin(k*x): oscillatory along one edge, growing across the
21+ % square by a factor exp(k), which is where its conditioning goes
22+ k = par;
23+ U = exp(k*Y) .* sin(k*X);
24+ case 'log'
25+ % log|z - z0| with z0 a distance par outside the right edge: harmonic
26+ % inside, but with a singularity just off the domain, so the solution
27+ % has structure on the length scale par near that edge
28+ z0 = (1 + par) + 0.5i;
29+ U = log(abs((X + 1i*Y) - z0));
30+ otherwise
31+ error('lp_exact2d: unknown solution "%s"', kind);
32+end
33+end
src/matlab/lp_lexact.madded+37−0View file
@@ -0,0 +1,37 @@
1+function U = lp_lexact(kind, X, Y)
2+% LP_LEXACT Exact solutions on the L-shaped domain (-1,1)^2 minus the quadrant
3+% x > 0, y < 0.
4+%
5+% The re-entrant corner sits at the origin, with interior angle 3*pi/2 between
6+% the edge y = 0, x > 0 and the edge x = 0, y < 0. Measuring the angle from the
7+% first of those, counterclockwise, the domain is 0 < theta < 3*pi/2, and
8+%
9+% s(r, theta) = r^(2/3) sin(2 theta / 3)
10+%
11+% is harmonic there and vanishes on both edges of the corner. It is the slowest
12+% decaying member of the family r^(2k/3) sin(2k theta/3) that a general
13+% Dirichlet datum excites, and it is the reason the corner is a problem: s is
14+% continuous, but its gradient is infinite at the origin and it has no second
15+% derivative there, so every convergence rate derived from a Taylor expansion
16+% is void.
17+%
18+% 'sing' the singular function alone
19+% 'smooth' log|z - z0| with z0 outside the domain: harmonic, and smooth right
20+% up to the corner, so the corner costs nothing
21+% 'both' the sum, which is the generic case
22+z0 = 1.6 + 1.3i;
23+th = atan2(Y, X);
24+th(th < -1e-13) = th(th < -1e-13) + 2*pi; % theta in [0, 3*pi/2] on the domain
25+r = sqrt(X.^2 + Y.^2);
26+sing = r.^(2/3) .* sin(2*th/3);
27+smooth = 0.5 * log((X - real(z0)).^2 + (Y - imag(z0)).^2);
28+
29+switch kind
30+ case 'sing'
31+ U = sing;
32+ case 'smooth'
33+ U = smooth;
34+ otherwise
35+ U = sing + smooth;
36+end
37+end
src/matlab/lp_orders.madded+14−0View file
@@ -0,0 +1,14 @@
1+function q = lp_orders(e, n)
2+% LP_ORDERS Observed convergence order between successive refinements.
3+%
4+% q(k) = log(e(k-1)/e(k)) / log(n(k)/n(k-1)), and NaN wherever the ratio is not
5+% a positive finite number (numbl's log of a NaN comes back complex, so the
6+% guard is not optional).
7+q = NaN(size(e));
8+for k = 2:numel(e)
9+ r = e(k-1) / e(k);
10+ if isfinite(r) && r > 0
11+ q(k) = log(r) / log(n(k) / n(k-1));
12+ end
13+end
14+end
src/matlab/lp_params.madded+4−0View file
@@ -0,0 +1,4 @@
1+function p = lp_params()
2+% LP_PARAMS The control settings for this run, as a struct.
3+p = jsondecode(fileread('params.json'));
4+end
src/matlab/lp_rows.madded+12−0View file
@@ -0,0 +1,12 @@
1+function C = lp_rows(M)
2+% LP_ROWS Matrix to a cell array of rows, for jsonencode.
3+%
4+% jsonencode([1 2 3]) is [1,2,3], not [[1,2,3]], so a matrix with a single row
5+% arrives on the JavaScript side with the wrong nesting. Encoding a cell array
6+% of rows instead always gives an array of arrays.
7+n = size(M, 1);
8+C = cell(1, n);
9+for i = 1:n
10+ C{i} = M(i, :);
11+end
12+end
src/matlab/lp_thin.madded+12−0View file
@@ -0,0 +1,12 @@
1+function [Ms, xs, ys] = lp_thin(M, x, y, maxn)
2+% LP_THIN Subsample a field so that a picture of it fits in a JSON message.
3+%
4+% Keeps at most maxn points in each direction, always including both ends.
5+nx = numel(x);
6+ny = numel(y);
7+ix = round(linspace(1, nx, min(nx, maxn)));
8+iy = round(linspace(1, ny, min(ny, maxn)));
9+Ms = M(iy, ix);
10+xs = x(ix);
11+ys = y(iy);
12+end
src/matlab/lp_write.madded+6−0View file
@@ -0,0 +1,6 @@
1+function lp_write(out)
2+% LP_WRITE Hand a struct back to the page as out.json.
3+f = fopen('out.json', 'w');
4+fwrite(f, jsonencode(out));
5+fclose(f);
6+end
src/plot/Curves.tsxadded+188−0View file
@@ -0,0 +1,188 @@
1+// Marks drawn inside a Plot frame: polylines, markers, reference slopes and
2+// direct end-labels. Every mark is thin (2px lines, >=8px markers) and every
3+// marker carries a shape as well as a colour.
4+import type { ReactNode } from 'react'
5+import type { Num } from '../engine/types.ts'
6+import { ink, surface } from './palette.ts'
7+import { linePath, type Frame } from './scales.ts'
8+
9+export type Shape = 'circle' | 'square' | 'triangle' | 'diamond' | 'cross' | 'plus' | 'star'
10+
11+export function Line({
12+ x,
13+ y,
14+ f,
15+ color,
16+ dash,
17+ width = 2,
18+ opacity = 1,
19+}: {
20+ x: readonly number[]
21+ y: readonly Num[]
22+ f: Frame
23+ color: string
24+ dash?: string
25+ width?: number
26+ opacity?: number
27+}) {
28+ return (
29+ <path
30+ d={linePath(x, y, f)}
31+ fill="none"
32+ stroke={color}
33+ strokeWidth={width}
34+ strokeDasharray={dash}
35+ strokeLinejoin="round"
36+ strokeLinecap="round"
37+ opacity={opacity}
38+ />
39+ )
40+}
41+
42+/** One marker, centred on (cx, cy) in pixels, about 9px across. */
43+export function markerPath(shape: Shape, cx: number, cy: number, r = 4.5): string {
44+ switch (shape) {
45+ case 'square':
46+ return `M${cx - r} ${cy - r}h${2 * r}v${2 * r}h${-2 * r}z`
47+ case 'triangle':
48+ return `M${cx} ${cy - r * 1.15}L${cx + r} ${cy + r * 0.8}L${cx - r} ${cy + r * 0.8}z`
49+ case 'diamond':
50+ return `M${cx} ${cy - r * 1.2}L${cx + r * 1.2} ${cy}L${cx} ${cy + r * 1.2}L${cx - r * 1.2} ${cy}z`
51+ case 'cross':
52+ return `M${cx - r} ${cy - r}L${cx + r} ${cy + r}M${cx + r} ${cy - r}L${cx - r} ${cy + r}`
53+ case 'plus':
54+ return `M${cx - r * 1.2} ${cy}h${2.4 * r}M${cx} ${cy - r * 1.2}v${2.4 * r}`
55+ case 'star': {
56+ let d = ''
57+ for (let k = 0; k < 3; k++) {
58+ const a = (Math.PI * k) / 3
59+ d += `M${cx - r * Math.cos(a)} ${cy - r * Math.sin(a)}L${cx + r * Math.cos(a)} ${cy + r * Math.sin(a)}`
60+ }
61+ return d
62+ }
63+ default:
64+ return `M${cx - r} ${cy}a${r} ${r} 0 1 0 ${2 * r} 0a${r} ${r} 0 1 0 ${-2 * r} 0z`
65+ }
66+}
67+
68+export function Markers({
69+ x,
70+ y,
71+ f,
72+ color,
73+ shape = 'circle',
74+ r = 4.5,
75+}: {
76+ x: readonly number[]
77+ y: readonly Num[]
78+ f: Frame
79+ color: string
80+ shape?: Shape
81+ r?: number
82+}) {
83+ const stroked = shape === 'cross' || shape === 'plus' || shape === 'star'
84+ const out: ReactNode[] = []
85+ for (let i = 0; i < x.length && i < y.length; i++) {
86+ const v = y[i]
87+ if (v == null || !isFinite(v) || !isFinite(x[i])) continue
88+ const px = f.sx(x[i])
89+ const py = f.sy(v)
90+ if (py < -20 || py > f.ih + 20) continue
91+ out.push(
92+ <path
93+ key={i}
94+ d={markerPath(shape, px, py, r)}
95+ fill={stroked ? 'none' : color}
96+ stroke={stroked ? color : surface}
97+ strokeWidth={stroked ? 2 : 1.5}
98+ />,
99+ )
100+ }
101+ return <g>{out}</g>
102+}
103+
104+/**
105+ * A dashed reference line of slope p through a chosen point, for log-log
106+ * convergence plots: the eye reads "is the data parallel to this?" much better
107+ * than it reads a printed rate.
108+ */
109+export function RefSlope({
110+ f,
111+ p,
112+ through,
113+ label,
114+ span = 0.6,
115+ labelAt = 'start',
116+}: {
117+ f: Frame
118+ /** the exponent: error ~ n^p */
119+ p: number
120+ /** a point the line passes through, in data coordinates */
121+ through: [number, number]
122+ label: string
123+ /** fraction of the x range the line covers */
124+ span?: number
125+ labelAt?: 'start' | 'end'
126+}) {
127+ const [lx, hx] = f.xDomain
128+ const l = Math.log10(lx)
129+ const h = Math.log10(hx)
130+ const c = Math.log10(through[0])
131+ const x0 = Math.pow(10, Math.max(l, c - (span * (h - l)) / 2))
132+ const x1 = Math.pow(10, Math.min(h, c + (span * (h - l)) / 2))
133+ const y = (xx: number) => through[1] * Math.pow(xx / through[0], p)
134+ return (
135+ <g>
136+ <path
137+ d={`M${f.sx(x0)} ${f.sy(y(x0))}L${f.sx(x1)} ${f.sy(y(x1))}`}
138+ stroke={ink.muted}
139+ strokeWidth={1.5}
140+ strokeDasharray="5 4"
141+ fill="none"
142+ />
143+ {labelAt === 'end' ? (
144+ <text x={f.sx(x1) - 3} y={f.sy(y(x1)) - 6} fill={ink.muted} fontSize={11} textAnchor="end">
145+ {label}
146+ </text>
147+ ) : (
148+ <text x={f.sx(x0) + 3} y={f.sy(y(x0)) + 14} fill={ink.muted} fontSize={11} textAnchor="start">
149+ {label}
150+ </text>
151+ )}
152+ </g>
153+ )
154+}
155+
156+/** A label placed at the end of a curve, so the reader need not use the legend. */
157+export function EndLabel({
158+ f,
159+ x,
160+ y,
161+ text,
162+ color,
163+ dx = 6,
164+ dy = 4,
165+ anchor = 'start',
166+}: {
167+ f: Frame
168+ x: number
169+ y: number
170+ text: string
171+ color: string
172+ dx?: number
173+ dy?: number
174+ anchor?: 'start' | 'end' | 'middle'
175+}) {
176+ return (
177+ <text
178+ x={f.sx(x) + dx}
179+ y={f.sy(y) + dy}
180+ fill={color}
181+ fontSize={11}
182+ textAnchor={anchor}
183+ style={{ paintOrder: 'stroke', stroke: surface, strokeWidth: 3 }}
184+ >
185+ {text}
186+ </text>
187+ )
188+}
src/plot/Field.tsxadded+213−0View file
@@ -0,0 +1,213 @@
1+// A scalar field on a rectangle, drawn as an image with a colour bar.
2+//
3+// The field is painted into a canvas one pixel per grid cell and then scaled to
4+// the frame, which keeps a 512-by-512 grid as cheap to draw as a 16-by-16 one.
5+// Cells outside the domain (an L-shaped region, the outside of a curve) are
6+// passed as null and left as surface.
7+import { useEffect, useRef, useState, useCallback } from 'react'
8+import type { Num } from '../engine/types.ts'
9+import { divColor, ink, seqColor, surface } from './palette.ts'
10+import { formatTick, niceTicks } from './scales.ts'
11+
12+export interface FieldProps {
13+ /** z[i][j] is the value at (x_j, y_i); null is "not in the domain" */
14+ z: Num[][]
15+ xDomain: [number, number]
16+ yDomain: [number, number]
17+ height?: number
18+ /** 'seq' for a magnitude, 'div' for a signed quantity, 'log' for |value| over decades */
19+ kind?: 'seq' | 'div' | 'log'
20+ /** upper end of the colour scale; taken from the data when omitted */
21+ zmax?: number
22+ /** how many decades the 'log' scale spans below zmax */
23+ decades?: number
24+ title?: string
25+ /** polylines in data coordinates, drawn over the field (a boundary, a mesh) */
26+ curves?: { x: number[]; y: number[]; color?: string; dash?: string }[]
27+ points?: { x: number[]; y: number[]; color?: string; r?: number }[]
28+ formatValue?: (v: number) => string
29+}
30+
31+function useWidth(): [(el: HTMLDivElement | null) => void, number] {
32+ const [width, setWidth] = useState(560)
33+ const obs = useRef<ResizeObserver | null>(null)
34+ const ref = useCallback((el: HTMLDivElement | null) => {
35+ obs.current?.disconnect()
36+ if (!el) return
37+ const ro = new ResizeObserver((entries) => {
38+ const w = entries[0]?.contentRect.width
39+ if (w && w > 0) setWidth(w)
40+ })
41+ ro.observe(el)
42+ obs.current = ro
43+ setWidth(el.clientWidth || 560)
44+ }, [])
45+ useEffect(() => () => obs.current?.disconnect(), [])
46+ return [ref, width]
47+}
48+
49+const M = { l: 44, r: 74, t: 18, b: 30 }
50+
51+export default function Field(props: FieldProps) {
52+ const { z, xDomain, yDomain, height = 340, kind = 'seq', decades = 8, title, curves, points } = props
53+ const [hostRef, width] = useWidth()
54+ const canvasRef = useRef<HTMLCanvasElement>(null)
55+
56+ const nrow = z.length
57+ const ncol = nrow > 0 ? z[0].length : 0
58+
59+ // colour scale limits
60+ let lo = Infinity
61+ let hi = -Infinity
62+ for (const row of z) {
63+ for (const v of row) {
64+ if (v == null || !isFinite(v)) continue
65+ if (v < lo) lo = v
66+ if (v > hi) hi = v
67+ }
68+ }
69+ if (!isFinite(lo)) {
70+ lo = 0
71+ hi = 1
72+ }
73+ const amax = props.zmax ?? Math.max(Math.abs(lo), Math.abs(hi))
74+ const logHi = props.zmax ?? (hi > 0 ? Math.log10(hi) : 0)
75+ const logLo = logHi - decades
76+
77+ const colorOf = (v: Num): string | null => {
78+ if (v == null || !isFinite(v)) return null
79+ if (kind === 'div') return divColor(amax > 0 ? v / amax : 0)
80+ if (kind === 'log') {
81+ const l = Math.log10(Math.max(Math.abs(v), Number.MIN_VALUE))
82+ return seqColor((l - logLo) / (logHi - logLo))
83+ }
84+ return seqColor(hi > lo ? (v - lo) / (hi - lo) : 0.5)
85+ }
86+
87+ const iw = Math.max(10, width - M.l - M.r)
88+ const ih = Math.max(10, height - M.t - M.b)
89+
90+ useEffect(() => {
91+ const cv = canvasRef.current
92+ if (!cv || ncol === 0) return
93+ cv.width = ncol
94+ cv.height = nrow
95+ const ctx = cv.getContext('2d')
96+ if (!ctx) return
97+ const img = ctx.createImageData(ncol, nrow)
98+ for (let i = 0; i < nrow; i++) {
99+ // row 0 of z is the bottom of the picture
100+ const dst = (nrow - 1 - i) * ncol
101+ for (let j = 0; j < ncol; j++) {
102+ const c = colorOf(z[i][j])
103+ const o = (dst + j) * 4
104+ if (c == null) {
105+ img.data[o] = 26
106+ img.data[o + 1] = 26
107+ img.data[o + 2] = 25
108+ img.data[o + 3] = 255
109+ } else {
110+ const n = parseInt(c.slice(1), 16)
111+ img.data[o] = (n >> 16) & 255
112+ img.data[o + 1] = (n >> 8) & 255
113+ img.data[o + 2] = n & 255
114+ img.data[o + 3] = 255
115+ }
116+ }
117+ }
118+ ctx.putImageData(img, 0, 0)
119+ })
120+
121+ const sx = (v: number) => ((v - xDomain[0]) / (xDomain[1] - xDomain[0])) * iw
122+ const sy = (v: number) => ih - ((v - yDomain[0]) / (yDomain[1] - yDomain[0])) * ih
123+
124+ const xt = niceTicks(xDomain[0], xDomain[1], 5)
125+ const yt = niceTicks(yDomain[0], yDomain[1], 5)
126+
127+ // colour bar: 24 steps from the same scale
128+ const bar = Array.from({ length: 24 }, (_, k) => k / 23)
129+ const barLabel = (t: number): string => {
130+ if (kind === 'div') return formatTick(Number(((2 * t - 1) * amax).toPrecision(2)))
131+ if (kind === 'log') {
132+ const e = logLo + t * (logHi - logLo)
133+ return `1e${Math.round(e)}`
134+ }
135+ return formatTick(Number((lo + t * (hi - lo)).toPrecision(2)))
136+ }
137+
138+ return (
139+ <div className="plot-host field-host" ref={hostRef}>
140+ <div style={{ position: 'relative', width, height }}>
141+ <canvas
142+ ref={canvasRef}
143+ style={{
144+ position: 'absolute',
145+ left: M.l,
146+ top: M.t,
147+ width: iw,
148+ height: ih,
149+ imageRendering: ncol < 96 ? 'pixelated' : 'auto',
150+ }}
151+ />
152+ <svg width={width} height={height} style={{ position: 'absolute', left: 0, top: 0 }} role="img">
153+ <g transform={`translate(${M.l},${M.t})`}>
154+ {curves?.map((c, k) => (
155+ <path
156+ key={`c${k}`}
157+ d={c.x.map((xv, i) => `${i ? 'L' : 'M'}${sx(xv).toFixed(2)} ${sy(c.y[i]).toFixed(2)}`).join('')}
158+ fill="none"
159+ stroke={c.color ?? ink.node}
160+ strokeWidth={1.5}
161+ strokeDasharray={c.dash}
162+ />
163+ ))}
164+ {points?.map((p, k) => (
165+ <g key={`p${k}`}>
166+ {p.x.map((xv, i) => (
167+ <circle key={i} cx={sx(xv)} cy={sy(p.y[i])} r={p.r ?? 1.8} fill={p.color ?? ink.node} />
168+ ))}
169+ </g>
170+ ))}
171+ <rect x={0} y={0} width={iw} height={ih} fill="none" stroke={ink.axis} strokeWidth={1} />
172+ {xt.map((v) => (
173+ <text key={`x${v}`} x={sx(v)} y={ih + 15} fill={ink.muted} fontSize={11} textAnchor="middle">
174+ {formatTick(v)}
175+ </text>
176+ ))}
177+ {yt.map((v) => (
178+ <text key={`y${v}`} x={-7} y={sy(v) + 4} fill={ink.muted} fontSize={11} textAnchor="end">
179+ {formatTick(v)}
180+ </text>
181+ ))}
182+ {title && (
183+ <text x={0} y={-6} fill={ink.secondary} fontSize={11}>
184+ {title}
185+ </text>
186+ )}
187+ {/* colour bar */}
188+ <g transform={`translate(${iw + 12},0)`}>
189+ {bar.map((t, k) => (
190+ <rect
191+ key={k}
192+ x={0}
193+ y={ih - ((k + 1) * ih) / bar.length}
194+ width={12}
195+ height={ih / bar.length + 0.5}
196+ fill={kind === 'div' ? divColor(2 * t - 1) : seqColor(t)}
197+ />
198+ ))}
199+ <rect x={0} y={0} width={12} height={ih} fill="none" stroke={ink.axis} strokeWidth={1} />
200+ {[0, 0.5, 1].map((t) => (
201+ <text key={t} x={16} y={ih - t * ih + 4} fill={ink.muted} fontSize={10}>
202+ {barLabel(t)}
203+ </text>
204+ ))}
205+ </g>
206+ </g>
207+ </svg>
208+ </div>
209+ </div>
210+ )
211+}
212+
213+export { surface }
src/plot/Legend.tsxadded+47−0View file
@@ -0,0 +1,47 @@
1+import { ink } from './palette.ts'
2+
3+export interface LegendItem {
4+ label: string
5+ color: string
6+ /** stroke-dasharray for a line swatch */
7+ dash?: string
8+ shape?: 'line' | 'dot' | 'cross'
9+ /** e.g. the max error for this series, shown after the label */
10+ value?: string
11+}
12+
13+/**
14+ * Identity is never carried by colour alone: every series that appears in a
15+ * plot appears here too, and the ones the reader most needs are also labelled
16+ * directly on the plot.
17+ */
18+export default function Legend({ items }: { items: LegendItem[] }) {
19+ return (
20+ <div className="legend">
21+ {items.map((it) => (
22+ <span className="legend-item" key={it.label}>
23+ <svg width={18} height={10} aria-hidden="true">
24+ {it.shape === 'dot' ? (
25+ <circle cx={9} cy={5} r={3.5} fill={it.color} stroke={ink.grid} strokeWidth={1} />
26+ ) : it.shape === 'cross' ? (
27+ <path d="M5 1 L13 9 M13 1 L5 9" stroke={it.color} strokeWidth={2} strokeLinecap="round" />
28+ ) : (
29+ <line
30+ x1={1}
31+ x2={17}
32+ y1={5}
33+ y2={5}
34+ stroke={it.color}
35+ strokeWidth={2}
36+ strokeDasharray={it.dash}
37+ strokeLinecap="round"
38+ />
39+ )}
40+ </svg>
41+ <span>{it.label}</span>
42+ {it.value && <span className="legend-value">{it.value}</span>}
43+ </span>
44+ ))}
45+ </div>
46+ )
47+}
src/plot/Plot.tsxadded+199−0View file
@@ -0,0 +1,199 @@
1+import { useCallback, useEffect, useId, useRef, useState, type ReactNode } from 'react'
2+import { ink } from './palette.ts'
3+import { decadeTicks, formatTick, niceTicks, type Frame } from './scales.ts'
4+
5+export interface PlotProps {
6+ height: number
7+ xDomain: [number, number]
8+ yDomain: [number, number]
9+ xLog?: boolean
10+ yLog?: boolean
11+ xLabel?: string
12+ yLabel?: string
13+ xTicks?: number[]
14+ yTicks?: number[]
15+ formatX?: (v: number) => string
16+ formatY?: (v: number) => string
17+ margin?: Partial<{ l: number; r: number; t: number; b: number }>
18+ /** drawn under the data, inside the frame but outside the clip */
19+ under?: (f: Frame) => ReactNode
20+ /** drawn over the data and outside the clip, so labels may reach the margin */
21+ over?: (f: Frame) => ReactNode
22+ children: (f: Frame) => ReactNode
23+ /** enables the crosshair; called with the data x under the pointer, or null */
24+ onHoverX?: (x: number | null) => void
25+ /** data x at which to draw the crosshair (usually what onHoverX last gave) */
26+ hoverX?: number | null
27+ /** rendered as an absolutely positioned box over the plot */
28+ overlay?: ReactNode
29+}
30+
31+const DEFAULT_MARGIN = { l: 54, r: 14, t: 10, b: 30 }
32+
33+/** Measures its own width so plots reflow with the panel. */
34+function useWidth(): [(el: HTMLDivElement | null) => void, number] {
35+ const [width, setWidth] = useState(640)
36+ const obs = useRef<ResizeObserver | null>(null)
37+ const ref = useCallback((el: HTMLDivElement | null) => {
38+ obs.current?.disconnect()
39+ if (!el) return
40+ const ro = new ResizeObserver((entries) => {
41+ const w = entries[0]?.contentRect.width
42+ if (w && w > 0) setWidth(w)
43+ })
44+ ro.observe(el)
45+ obs.current = ro
46+ setWidth(el.clientWidth || 640)
47+ }, [])
48+ useEffect(() => () => obs.current?.disconnect(), [])
49+ return [ref, width]
50+}
51+
52+export default function Plot(props: PlotProps) {
53+ const {
54+ height,
55+ xDomain,
56+ yDomain,
57+ xLog = false,
58+ yLog = false,
59+ xLabel,
60+ yLabel,
61+ formatX = formatTick,
62+ formatY = formatTick,
63+ under,
64+ over,
65+ children,
66+ onHoverX,
67+ hoverX,
68+ overlay,
69+ } = props
70+ const m = { ...DEFAULT_MARGIN, ...props.margin }
71+ const [hostRef, width] = useWidth()
72+ const clipId = useId().replace(/:/g, '')
73+ const svgRef = useRef<SVGSVGElement>(null)
74+
75+ const iw = Math.max(10, width - m.l - m.r)
76+ const ih = Math.max(10, height - m.t - m.b)
77+
78+ const fwd = (v: number, [lo, hi]: [number, number], log: boolean, span: number, flip: boolean) => {
79+ const t = log
80+ ? (Math.log10(Math.max(v, Number.MIN_VALUE)) - Math.log10(lo)) / (Math.log10(hi) - Math.log10(lo))
81+ : (v - lo) / (hi - lo)
82+ return flip ? span * (1 - t) : span * t
83+ }
84+
85+ const f: Frame = {
86+ sx: (v) => fwd(v, xDomain, xLog, iw, false),
87+ sy: (v) => fwd(v, yDomain, yLog, ih, true),
88+ ix: (px) => {
89+ const t = px / iw
90+ return xLog
91+ ? Math.pow(10, Math.log10(xDomain[0]) + t * (Math.log10(xDomain[1]) - Math.log10(xDomain[0])))
92+ : xDomain[0] + t * (xDomain[1] - xDomain[0])
93+ },
94+ iw,
95+ ih,
96+ xDomain,
97+ yDomain,
98+ }
99+
100+ const xt = props.xTicks ?? (xLog ? decadeTicks(xDomain[0], xDomain[1]) : niceTicks(xDomain[0], xDomain[1], 7))
101+ const yt = props.yTicks ?? (yLog ? decadeTicks(yDomain[0], yDomain[1]) : niceTicks(yDomain[0], yDomain[1], 5))
102+
103+ const pointer = (e: React.PointerEvent) => {
104+ if (!onHoverX) return
105+ const rect = svgRef.current?.getBoundingClientRect()
106+ if (!rect) return
107+ 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))))
109+ }
110+
111+ return (
112+ <div className="plot-host" ref={hostRef}>
113+ <svg
114+ ref={svgRef}
115+ width={width}
116+ height={height}
117+ role="img"
118+ onPointerMove={pointer}
119+ onPointerLeave={() => onHoverX?.(null)}
120+ >
121+ <defs>
122+ <clipPath id={clipId}>
123+ <rect x={0} y={0} width={iw} height={ih} />
124+ </clipPath>
125+ </defs>
126+ <g transform={`translate(${m.l},${m.t})`}>
127+ {/* gridlines, recessive */}
128+ {xt.map((v) => (
129+ <line key={`gx${v}`} x1={f.sx(v)} x2={f.sx(v)} y1={0} y2={ih} stroke={ink.grid} strokeWidth={1} />
130+ ))}
131+ {yt.map((v) => (
132+ <line key={`gy${v}`} x1={0} x2={iw} y1={f.sy(v)} y2={f.sy(v)} stroke={ink.grid} strokeWidth={1} />
133+ ))}
134+
135+ {under?.(f)}
136+
137+ <g clipPath={`url(#${clipId})`}>{children(f)}</g>
138+
139+ {over?.(f)}
140+
141+ {hoverX != null && hoverX >= Math.min(...xDomain) && hoverX <= Math.max(...xDomain) && (
142+ <line
143+ x1={f.sx(hoverX)}
144+ x2={f.sx(hoverX)}
145+ y1={0}
146+ y2={ih}
147+ stroke={ink.secondary}
148+ strokeWidth={1}
149+ strokeDasharray="3 3"
150+ opacity={0.7}
151+ pointerEvents="none"
152+ />
153+ )}
154+
155+ {/* axes */}
156+ <line x1={0} x2={iw} y1={ih} y2={ih} stroke={ink.axis} strokeWidth={1} />
157+ <line x1={0} x2={0} y1={0} y2={ih} stroke={ink.axis} strokeWidth={1} />
158+ {xt.map((v) => (
159+ <text
160+ key={`tx${v}`}
161+ x={f.sx(v)}
162+ y={ih + 15}
163+ fill={ink.muted}
164+ fontSize={11}
165+ textAnchor="middle"
166+ style={{ fontVariantNumeric: 'tabular-nums' }}
167+ >
168+ {formatX(v)}
169+ </text>
170+ ))}
171+ {yt.map((v) => (
172+ <text
173+ key={`ty${v}`}
174+ x={-7}
175+ y={f.sy(v) + 4}
176+ fill={ink.muted}
177+ fontSize={11}
178+ textAnchor="end"
179+ style={{ fontVariantNumeric: 'tabular-nums' }}
180+ >
181+ {formatY(v)}
182+ </text>
183+ ))}
184+ {xLabel && (
185+ <text x={iw} y={ih + 27} fill={ink.muted} fontSize={11} textAnchor="end">
186+ {xLabel}
187+ </text>
188+ )}
189+ {yLabel && (
190+ <text x={-m.l + 4} y={-1} fill={ink.muted} fontSize={11} textAnchor="start">
191+ {yLabel}
192+ </text>
193+ )}
194+ </g>
195+ </svg>
196+ {overlay}
197+ </div>
198+ )
199+}
src/plot/palette.tsadded+127−0View file
@@ -0,0 +1,127 @@
1+// Colours, assigned by the job each one does, taken unchanged from the
2+// reference data-viz palette (dark steps, dark surface #1a1a19) and validated
3+// as a set:
4+//
5+// node scripts/validate_palette.js \
6+// "#3987e5,#d95926,#199e70,#c98500,#d55181,#008300,#9085e9,#e66767" --mode dark
7+// lightness band PASS chroma floor PASS CVD separation PASS (worst 8.4)
8+// normal-vision floor PASS (worst 19.3) contrast vs surface PASS
9+//
10+// The worst adjacent pair sits in the 6-8 CVD band's neighbourhood, so every
11+// plot that shows more than two methods also carries a legend, a distinct
12+// marker shape per method, and direct end-labels: identity is never colour
13+// alone. If you add a series, re-run the validator rather than picking a hex
14+// that looks nice.
15+
16+export const surface = '#1a1a19'
17+
18+/**
19+ * Identity: which method produced a curve. Fixed slots in the palette's
20+ * documented order; never cycled, never reassigned when a plot shows a subset.
21+ */
22+export const method = {
23+ /** second-order finite differences */
24+ fd2: '#3987e5',
25+ /** fourth-order finite differences (1-D) and the compact 9-point stencil (2-D) */
26+ fd4: '#d95926',
27+ /** Chebyshev spectral collocation */
28+ spectral: '#199e70',
29+ /** the FFT-based fast Poisson solver */
30+ fast: '#c98500',
31+ /** Krylov and multigrid iterations */
32+ iterative: '#d55181',
33+ /** the method of particular solutions */
34+ mps: '#9085e9',
35+ /** boundary integral equations */
36+ bie: '#e66767',
37+} as const
38+
39+export type MethodKey = keyof typeof method
40+
41+/** Marker shapes, so that identity survives greyscale and colour blindness. */
42+export const marker: Record<MethodKey, 'circle' | 'square' | 'triangle' | 'diamond' | 'cross' | 'plus' | 'star'> = {
43+ fd2: 'circle',
44+ fd4: 'square',
45+ spectral: 'triangle',
46+ fast: 'diamond',
47+ iterative: 'cross',
48+ mps: 'plus',
49+ bie: 'star',
50+}
51+
52+/** State. Reserved: never reused as a series colour. */
53+export const status = {
54+ good: '#008300',
55+ warning: '#c98500',
56+ critical: '#e66767',
57+} as const
58+
59+/** Sequential ramp for magnitude (a field, an error, a residual): one hue. */
60+export const sequential = [
61+ '#cde2fb',
62+ '#b7d3f6',
63+ '#9ec5f4',
64+ '#86b6ef',
65+ '#6da7ec',
66+ '#5598e7',
67+ '#3987e5',
68+ '#2a78d6',
69+ '#256abf',
70+ '#1c5cab',
71+ '#184f95',
72+ '#104281',
73+ '#0d366b',
74+] as const
75+
76+/** Diverging pair for a signed field: two poles, neutral grey midpoint. */
77+export const diverging = {
78+ low: '#3987e5',
79+ mid: '#383835',
80+ high: '#e66767',
81+} as const
82+
83+/** Chart chrome and ink. Text always wears these, never a series colour. */
84+export const ink = {
85+ primary: '#e8e6df',
86+ secondary: '#c3c2b7',
87+ muted: '#8a8880',
88+ grid: '#272724',
89+ axis: '#3c3c38',
90+ /** the exact solution: the reference everything is measured against */
91+ reference: '#8a8880',
92+ /** grid nodes, boundary points, collocation points */
93+ node: '#e8e6df',
94+ /** a family of curves drawn as a mass */
95+ family: '#4a4a45',
96+} as const
97+
98+function hexToRgb(hex: string): [number, number, number] {
99+ const v = parseInt(hex.slice(1), 16)
100+ return [(v >> 16) & 255, (v >> 8) & 255, v & 255]
101+}
102+
103+function rgbToHex(r: number, g: number, b: number): string {
104+ const c = (x: number) => Math.max(0, Math.min(255, Math.round(x))).toString(16).padStart(2, '0')
105+ return `#${c(r)}${c(g)}${c(b)}`
106+}
107+
108+function mix(a: string, b: string, t: number): string {
109+ const [r1, g1, b1] = hexToRgb(a)
110+ const [r2, g2, b2] = hexToRgb(b)
111+ return rgbToHex(r1 + (r2 - r1) * t, g1 + (g2 - g1) * t, b1 + (b2 - b1) * t)
112+}
113+
114+/** Sequential colour for t in [0, 1]: light at 0 (near zero), dark at 1. */
115+export function seqColor(t: number): string {
116+ if (!isFinite(t)) return ink.grid
117+ const u = Math.max(0, Math.min(1, t)) * (sequential.length - 1)
118+ const i = Math.floor(u)
119+ return i >= sequential.length - 1 ? sequential[sequential.length - 1] : mix(sequential[i], sequential[i + 1], u - i)
120+}
121+
122+/** Diverging colour for t in [-1, 1]: blue below, grey at zero, red above. */
123+export function divColor(t: number): string {
124+ if (!isFinite(t)) return ink.grid
125+ const u = Math.max(-1, Math.min(1, t))
126+ return u < 0 ? mix(diverging.mid, diverging.low, -u) : mix(diverging.mid, diverging.high, u)
127+}
src/plot/scales.tsadded+114−0View file
@@ -0,0 +1,114 @@
1+import type { Num } from '../engine/types.ts'
2+
3+export interface Frame {
4+ /** data x to pixel x, within the plotting area */
5+ sx: (v: number) => number
6+ /** data y to pixel y */
7+ sy: (v: number) => number
8+ /** pixel x back to data x, for hover */
9+ ix: (px: number) => number
10+ iw: number
11+ ih: number
12+ xDomain: [number, number]
13+ yDomain: [number, number]
14+}
15+
16+export function niceTicks(lo: number, hi: number, target = 6): number[] {
17+ if (!isFinite(lo) || !isFinite(hi) || hi <= lo) return []
18+ const raw = (hi - lo) / target
19+ const mag = Math.pow(10, Math.floor(Math.log10(raw)))
20+ const norm = raw / mag
21+ const step = (norm >= 5 ? 10 : norm >= 2 ? 5 : norm >= 1 ? 2 : 1) * mag
22+ const out: number[] = []
23+ for (let v = Math.ceil(lo / step) * step; v <= hi + step * 1e-9; v += step) {
24+ out.push(Math.abs(v) < step * 1e-9 ? 0 : v)
25+ }
26+ return out
27+}
28+
29+/** Decade ticks, thinned so the labels do not collide over a wide range. */
30+export function decadeTicks(lo: number, hi: number, maxCount = 9): number[] {
31+ if (!(lo > 0) || !(hi > 0)) return []
32+ const a = Math.floor(Math.log10(lo))
33+ const b = Math.ceil(Math.log10(hi))
34+ const every = Math.max(1, Math.ceil((b - a + 1) / maxCount))
35+ const out: number[] = []
36+ for (let e = a; e <= b; e += every) {
37+ const v = Math.pow(10, e)
38+ if (v >= lo * 0.999 && v <= hi * 1.001) out.push(v)
39+ }
40+ return out
41+}
42+
43+export function formatTick(v: number): string {
44+ if (v === 0) return '0'
45+ const a = Math.abs(v)
46+ if (a >= 1e5 || a < 1e-3) {
47+ const e = Math.round(Math.log10(a))
48+ if (Math.abs(a - Math.pow(10, e)) < Math.pow(10, e) * 1e-6) return `1e${e}`
49+ return v.toExponential(0)
50+ }
51+ if (Number.isInteger(v)) return String(v)
52+ return String(Number(v.toPrecision(3)))
53+}
54+
55+export function padDomain([lo, hi]: [number, number], frac = 0.06): [number, number] {
56+ if (!isFinite(lo) || !isFinite(hi)) return [0, 1]
57+ if (hi === lo) return [lo - 0.5, hi + 0.5]
58+ const p = (hi - lo) * frac
59+ return [lo - p, hi + p]
60+}
61+
62+/** Range of the finite values in one or more series. */
63+export function extent(...series: (readonly Num[] | undefined)[]): [number, number] {
64+ let lo = Infinity
65+ let hi = -Infinity
66+ for (const s of series) {
67+ if (!s) continue
68+ for (const v of s) {
69+ if (v == null || !isFinite(v)) continue
70+ if (v < lo) lo = v
71+ if (v > hi) hi = v
72+ }
73+ }
74+ return isFinite(lo) ? [lo, hi] : [0, 1]
75+}
76+
77+/**
78+ * A polyline through (xs, ys). Runs of missing values break the path rather
79+ * than being bridged, and values far outside the frame are clamped to a band
80+ * just off-screen so that a curve running off to infinity still leaves the
81+ * frame in the right direction instead of producing unusable path data.
82+ */
83+export function linePath(xs: readonly number[], ys: readonly Num[], f: Frame): string {
84+ const bound = f.ih * 12
85+ let out = ''
86+ let pen = false
87+ for (let i = 0; i < xs.length && i < ys.length; i++) {
88+ const y = ys[i]
89+ if (y == null || !isFinite(y) || !isFinite(xs[i])) {
90+ pen = false
91+ continue
92+ }
93+ const px = f.sx(xs[i])
94+ const py = Math.min(bound, Math.max(-bound, f.sy(y)))
95+ out += `${pen ? 'L' : 'M'}${px.toFixed(2)} ${py.toFixed(2)}`
96+ pen = true
97+ }
98+ return out
99+}
100+
101+/** Index of the sample nearest a data-space x. */
102+export function nearestIndex(xs: readonly number[], v: number): number {
103+ if (xs.length === 0) return -1
104+ let lo = 0
105+ let hi = xs.length - 1
106+ if (v <= xs[lo]) return lo
107+ if (v >= xs[hi]) return hi
108+ while (hi - lo > 1) {
109+ const mid = (lo + hi) >> 1
110+ if (xs[mid] <= v) lo = mid
111+ else hi = mid
112+ }
113+ return v - xs[lo] <= xs[hi] - v ? lo : hi
114+}
src/ui/ChapterView.tsxadded+46−0View file
@@ -0,0 +1,46 @@
1+// A chapter: prose with its labs dropped in where the text calls for them.
2+//
3+// The prose is one markdown file per chapter, with a line
4+//
5+// <!--lab:convergence-->
6+//
7+// wherever a lab belongs. Splitting on those markers keeps the writing in one
8+// readable file instead of scattered through JSX.
9+import { useMemo } from 'react'
10+import Markdown from './Markdown.tsx'
11+import Lab from './Lab.tsx'
12+import type { Chapter } from './types.ts'
13+
14+const MARKER = /<!--\s*lab:([a-z0-9_-]+)\s*-->/gi
15+
16+export default function ChapterView({ chapter }: { chapter: Chapter }) {
17+ const pieces = useMemo(() => {
18+ const out: ({ md: string } | { labId: string })[] = []
19+ let last = 0
20+ MARKER.lastIndex = 0
21+ let m: RegExpExecArray | null
22+ while ((m = MARKER.exec(chapter.text)) !== null) {
23+ out.push({ md: chapter.text.slice(last, m.index) })
24+ out.push({ labId: m[1] })
25+ last = m.index + m[0].length
26+ }
27+ out.push({ md: chapter.text.slice(last) })
28+ return out
29+ }, [chapter.text])
30+
31+ return (
32+ <article className="chapter">
33+ <header className="chapter-head">
34+ <div className="chapter-part">{chapter.part}</div>
35+ <h2>{chapter.title}</h2>
36+ <p className="chapter-blurb">{chapter.blurb}</p>
37+ </header>
38+ {pieces.map((p, i) => {
39+ if ('md' in p) return p.md.trim() === '' ? null : <Markdown key={i} src={p.md} />
40+ const lab = chapter.labs.find((l) => l.id === p.labId)
41+ if (!lab) return <p key={i} className="lab-error">missing lab: {p.labId}</p>
42+ return <Lab key={`${chapter.id}-${lab.id}`} spec={lab} />
43+ })}
44+ </article>
45+ )
46+}
src/ui/Controls.tsxadded+93−0View file
@@ -0,0 +1,93 @@
1+// The controls for one lab, declared by the lab and laid out in a single row
2+// above its figures.
3+import type { Params, ParamValue } from '../engine/types.ts'
4+import type { Control } from './types.ts'
5+
6+export default function Controls({
7+ controls,
8+ params,
9+ onChange,
10+}: {
11+ controls: Control[]
12+ params: Params
13+ onChange: (key: string, value: ParamValue) => void
14+}) {
15+ if (controls.length === 0) return null
16+ return (
17+ <div className="controls">
18+ {controls.map((c) => {
19+ if (c.kind === 'select') {
20+ return (
21+ <label className="control" key={c.key}>
22+ <span className="control-label">{c.label}</span>
23+ <select
24+ value={String(params[c.key])}
25+ onChange={(e) => {
26+ const raw = e.target.value
27+ const opt = c.options.find((o) => String(o.value) === raw)
28+ onChange(c.key, opt ? opt.value : raw)
29+ }}
30+ >
31+ {c.options.map((o) => (
32+ <option key={String(o.value)} value={String(o.value)}>
33+ {o.label}
34+ </option>
35+ ))}
36+ </select>
37+ </label>
38+ )
39+ }
40+ if (c.kind === 'toggle') {
41+ return (
42+ <label className="control control-toggle" key={c.key}>
43+ <input
44+ type="checkbox"
45+ checked={Boolean(params[c.key])}
46+ onChange={(e) => onChange(c.key, e.target.checked)}
47+ />
48+ <span className="control-label">{c.label}</span>
49+ </label>
50+ )
51+ }
52+ const value = Number(params[c.key])
53+ const fmt = c.format ?? ((v: number) => String(v))
54+ if (c.values) {
55+ const idx = Math.max(
56+ 0,
57+ c.values.findIndex((v) => v === value),
58+ )
59+ return (
60+ <label className="control" key={c.key}>
61+ <span className="control-label">
62+ {c.label} <span className="control-value">{fmt(value)}</span>
63+ </span>
64+ <input
65+ type="range"
66+ min={0}
67+ max={c.values.length - 1}
68+ step={1}
69+ value={idx}
70+ onChange={(e) => onChange(c.key, c.values![Number(e.target.value)])}
71+ />
72+ </label>
73+ )
74+ }
75+ return (
76+ <label className="control" key={c.key}>
77+ <span className="control-label">
78+ {c.label} <span className="control-value">{fmt(value)}</span>
79+ </span>
80+ <input
81+ type="range"
82+ min={c.min}
83+ max={c.max}
84+ step={c.step ?? 1}
85+ value={value}
86+ onChange={(e) => onChange(c.key, Number(e.target.value))}
87+ />
88+ </label>
89+ )
90+ })}
91+ </div>
92+ )
93+}
src/ui/Lab.tsxadded+144−0View file
@@ -0,0 +1,144 @@
1+// One runnable section: controls, figures, and the MATLAB behind them.
2+//
3+// Each lab owns a numbl session for as long as it is on screen. Moving a
4+// control rewrites params.json and re-runs; editing the script reboots the
5+// session, which is why an edited lab waits for the Run button.
6+import { useCallback, useEffect, useMemo, useRef, useState } from 'react'
7+import ScriptEditor from '../editor/ScriptEditor.tsx'
8+import Controls from './Controls.tsx'
9+import { Engine } from '../engine/runner.ts'
10+import type { Params, ParamValue, RunResult } from '../engine/types.ts'
11+import type { LabSpec } from './types.ts'
12+
13+export default function Lab({ spec }: { spec: LabSpec }) {
14+ const [script, setScript] = useState(spec.script)
15+ const [pending, setPending] = useState(spec.script)
16+ const [params, setParams] = useState<Params>(spec.defaults)
17+ const [result, setResult] = useState<RunResult<unknown> | null>(null)
18+ const [running, setRunning] = useState(false)
19+ const [ranOnce, setRanOnce] = useState(false)
20+ const [showEditor, setShowEditor] = useState(Boolean(spec.editorOpen))
21+ const [showLog, setShowLog] = useState(false)
22+ const engine = useMemo(() => new Engine(), [])
23+ const seq = useRef(0)
24+
25+ useEffect(() => () => engine.dispose(), [engine])
26+
27+ const run = useCallback(
28+ async (src: string, p: Params) => {
29+ const my = ++seq.current
30+ setRunning(true)
31+ const res = await engine.run(spec.driver, src, spec.derive ? spec.derive(p) : p)
32+ if (seq.current !== my) return
33+ setResult(res)
34+ setRunning(false)
35+ setRanOnce(true)
36+ },
37+ [engine, spec.driver, spec.derive],
38+ )
39+
40+ // first run, unless the lab is heavy enough to wait to be asked
41+ useEffect(() => {
42+ if (!spec.manual) void run(script, params)
43+ // eslint-disable-next-line react-hooks/exhaustive-deps
44+ }, [])
45+
46+ // a control moved: re-run, debounced, but only once the lab is live
47+ const paramsRef = useRef(params)
48+ paramsRef.current = params
49+ const onControl = useCallback(
50+ (key: string, value: ParamValue) => {
51+ const next = { ...paramsRef.current, [key]: value }
52+ paramsRef.current = next
53+ setParams(next)
54+ },
55+ [],
56+ )
57+ const firstParams = useRef(true)
58+ useEffect(() => {
59+ if (firstParams.current) {
60+ firstParams.current = false
61+ return
62+ }
63+ if (!ranOnce) return
64+ const t = setTimeout(() => void run(script, params), 120)
65+ return () => clearTimeout(t)
66+ // eslint-disable-next-line react-hooks/exhaustive-deps
67+ }, [params])
68+
69+ const dirty = pending !== script
70+ const onRun = useCallback(() => {
71+ setScript(pending)
72+ void run(pending, paramsRef.current)
73+ }, [pending, run])
74+
75+ const data = result && result.ok ? result.data : null
76+ const summary = data && spec.summary ? spec.summary(data, params) : []
77+
78+ return (
79+ <section className="lab" id={`lab-${spec.id}`}>
80+ <header className="lab-head">
81+ <h3>{spec.title}</h3>
82+ {spec.caption && <p className="lab-caption">{spec.caption}</p>}
83+ </header>
84+
85+ <Controls controls={spec.controls} params={params} onChange={onControl} />
86+
87+ <div className="lab-bar">
88+ <button className={dirty || !ranOnce ? 'primary' : ''} onClick={onRun} disabled={running}>
89+ {running ? 'running...' : dirty ? 'run edited script' : 'run'}
90+ </button>
91+ {result && (
92+ <span className="stat">
93+ <span className="stat-label">wall clock</span>
94+ <span className="stat-value">{result.ms < 1000 ? `${result.ms.toFixed(0)} ms` : `${(result.ms / 1000).toFixed(2)} s`}</span>
95+ </span>
96+ )}
97+ {summary.map((s) => (
98+ <span className="stat" key={s.label}>
99+ <span className="stat-label">{s.label}</span>
100+ <span className="stat-value">{s.value}</span>
101+ </span>
102+ ))}
103+ <span className="lab-bar-gap" />
104+ <button className="ghost" onClick={() => setShowEditor((v) => !v)}>
105+ {showEditor ? 'hide script' : 'show script'}
106+ </button>
107+ {result && result.output.trim() !== '' && (
108+ <button className="ghost" onClick={() => setShowLog((v) => !v)}>
109+ {showLog ? 'hide output' : 'output'}
110+ </button>
111+ )}
112+ </div>
113+
114+ {result && !result.ok && <pre className="lab-error">{result.error}</pre>}
115+
116+ <div className="lab-figures">
117+ {data ? spec.render(data, params) : <div className="lab-placeholder">{running ? 'running' : 'press run'}</div>}
118+ </div>
119+
120+ {showLog && result && <pre className="lab-log">{result.output}</pre>}
121+
122+ {showEditor && (
123+ <div className="lab-editor">
124+ <div className="lab-editor-bar">
125+ <span>the method, in MATLAB syntax</span>
126+ <span className="lab-bar-gap" />
127+ {dirty && <span className="warn">edited; press run</span>}
128+ <button
129+ className="ghost"
130+ onClick={() => {
131+ setPending(spec.script)
132+ setScript(spec.script)
133+ void run(spec.script, paramsRef.current)
134+ }}
135+ >
136+ reset
137+ </button>
138+ </div>
139+ <ScriptEditor value={pending} onChange={setPending} onRun={onRun} />
140+ </div>
141+ )}
142+ </section>
143+ )
144+}
src/ui/LibraryView.tsxadded+39−0View file
@@ -0,0 +1,39 @@
1+// The shared .m files, shown read-only. The editable script in each lab is
2+// the method itself; these are the pieces every lab is allowed to assume.
3+import { useState } from 'react'
4+
5+const raw = import.meta.glob('../matlab/*.m', { query: '?raw', import: 'default', eager: true }) as Record<
6+ string,
7+ string
8+>
9+
10+const FILES = Object.entries(raw)
11+ .map(([path, content]) => ({ name: path.slice(path.lastIndexOf('/') + 1), content }))
12+ .sort((a, b) => a.name.localeCompare(b.name))
13+
14+export default function LibraryView() {
15+ const [open, setOpen] = useState<string | null>(FILES[0]?.name ?? null)
16+ return (
17+ <article className="chapter">
18+ <header className="chapter-head">
19+ <div className="chapter-part">Reference</div>
20+ <h2>The shared MATLAB library</h2>
21+ <p className="chapter-blurb">
22+ Every lab boots with these function files in its session, alongside the editable script and a driver that
23+ reads the controls and writes the figures. They are the small pieces that would otherwise be repeated in
24+ every chapter.
25+ </p>
26+ </header>
27+ <div className="lib-list">
28+ {FILES.map((f) => (
29+ <div className="lib-file" key={f.name}>
30+ <button className="lib-name" onClick={() => setOpen(open === f.name ? null : f.name)}>
31+ {open === f.name ? '▾' : '▸'} {f.name}
32+ </button>
33+ {open === f.name && <pre className="lib-src">{f.content}</pre>}
34+ </div>
35+ ))}
36+ </div>
37+ </article>
38+ )
39+}
src/ui/Markdown.tsxadded+35−0View file
@@ -0,0 +1,35 @@
1+// Prose is written as markdown with TeX between dollar signs. KaTeX renders
2+// the maths first, into placeholders that marked will not touch, and the
3+// result is put in with dangerouslySetInnerHTML: the source is ours, written
4+// at build time, not anything a reader supplies.
5+import { useMemo } from 'react'
6+import { marked } from 'marked'
7+import katex from 'katex'
8+import 'katex/dist/katex.min.css'
9+
10+function renderMath(src: string): { text: string; slots: string[] } {
11+ const slots: string[] = []
12+ const keep = (html: string) => {
13+ slots.push(html)
14+ return `@@MATH${slots.length - 1}@@`
15+ }
16+ // display maths first, so that $$...$$ is not eaten by the inline rule
17+ let text = src.replace(/\$\$([\s\S]+?)\$\$/g, (_m, tex: string) =>
18+ keep(katex.renderToString(tex.trim(), { displayMode: true, throwOnError: false })),
19+ )
20+ text = text.replace(/(^|[^\\$])\$([^$\n]+?)\$/g, (_m, pre: string, tex: string) =>
21+ pre + keep(katex.renderToString(tex.trim(), { displayMode: false, throwOnError: false })),
22+ )
23+ return { text, slots }
24+}
25+
26+export function markdownToHtml(src: string): string {
27+ const { text, slots } = renderMath(src)
28+ const html = marked.parse(text, { async: false }) as string
29+ return html.replace(/@@MATH(\d+)@@/g, (_m, i: string) => slots[Number(i)])
30+}
31+
32+export default function Markdown({ src }: { src: string }) {
33+ const html = useMemo(() => markdownToHtml(src), [src])
34+ return <div className="prose" dangerouslySetInnerHTML={{ __html: html }} />
35+}
src/ui/types.tsadded+59−0View file
@@ -0,0 +1,59 @@
1+import type { ReactNode } from 'react'
2+import type { Params } from '../engine/types.ts'
3+
4+export type Control =
5+ | {
6+ kind: 'slider'
7+ key: string
8+ label: string
9+ min: number
10+ max: number
11+ step?: number
12+ /** discrete choices; the slider then moves over their indices */
13+ values?: number[]
14+ format?: (v: number) => string
15+ }
16+ | {
17+ kind: 'select'
18+ key: string
19+ label: string
20+ options: { value: string | number; label: string }[]
21+ }
22+ | { kind: 'toggle'; key: string; label: string }
23+
24+export interface LabSpec {
25+ id: string
26+ title: string
27+ /** one line telling the reader what to watch */
28+ caption?: string
29+ /** the fixed part of the MATLAB layer, hidden from the reader */
30+ driver: string
31+ /** the editable part: the method itself */
32+ script: string
33+ controls: Control[]
34+ defaults: Params
35+ /** turns the control settings into what params.json should contain */
36+ derive?: (p: Params) => Params
37+ /** heavy labs wait for the Run button instead of running on sight */
38+ manual?: boolean
39+ /** the editor starts open when the script is the point of the section */
40+ editorOpen?: boolean
41+ // eslint-disable-next-line @typescript-eslint/no-explicit-any
42+ render: (data: any, params: Params) => ReactNode
43+ /** numbers for the status strip beside the timing */
44+ // eslint-disable-next-line @typescript-eslint/no-explicit-any
45+ summary?: (data: any, params: Params) => { label: string; value: string }[]
46+}
47+
48+export interface Chapter {
49+ id: string
50+ /** short name for the sidebar */
51+ nav: string
52+ title: string
53+ /** one sentence, shown under the title and in the contents */
54+ blurb: string
55+ part: 'One dimension' | 'Two dimensions' | 'Beyond the square'
56+ /** the prose, with <!--lab:id--> lines marking where each lab goes */
57+ text: string
58+ labs: LabSpec[]
59+}
src/vite-env.d.tsadded+11−0View file
@@ -0,0 +1,11 @@
1+/// <reference types="vite/client" />
2+
3+declare module '*.m?raw' {
4+ const src: string
5+ export default src
6+}
7+
8+declare module '*.md?raw' {
9+ const src: string
10+ export default src
11+}
tsconfig.app.jsonadded+25−0View file
@@ -0,0 +1,25 @@
1+{
2+ "compilerOptions": {
3+ "tsBuildInfoFile": "./node_modules/.tmp/tsconfig.app.tsbuildinfo",
4+ "target": "es2023",
5+ "lib": ["ES2023", "DOM", "WebWorker"],
6+ "module": "esnext",
7+ "types": ["vite/client"],
8+ "allowArbitraryExtensions": true,
9+ "skipLibCheck": true,
10+
11+ "moduleResolution": "bundler",
12+ "allowImportingTsExtensions": true,
13+ "verbatimModuleSyntax": true,
14+ "moduleDetection": "force",
15+ "noEmit": true,
16+ "jsx": "react-jsx",
17+
18+ "strict": true,
19+ "noUnusedLocals": true,
20+ "noUnusedParameters": true,
21+ "erasableSyntaxOnly": true,
22+ "noFallthroughCasesInSwitch": true
23+ },
24+ "include": ["src"]
25+}
tsconfig.jsonadded+4−0View file
@@ -0,0 +1,4 @@
1+{
2+ "files": [],
3+ "references": [{ "path": "./tsconfig.app.json" }, { "path": "./tsconfig.node.json" }]
4+}
tsconfig.node.jsonadded+21−0View file
@@ -0,0 +1,21 @@
1+{
2+ "compilerOptions": {
3+ "tsBuildInfoFile": "./node_modules/.tmp/tsconfig.node.tsbuildinfo",
4+ "target": "es2023",
5+ "lib": ["ES2023"],
6+ "module": "esnext",
7+ "types": ["node"],
8+ "skipLibCheck": true,
9+
10+ "moduleResolution": "bundler",
11+ "allowImportingTsExtensions": true,
12+ "verbatimModuleSyntax": true,
13+ "moduleDetection": "force",
14+ "noEmit": true,
15+
16+ "strict": true,
17+ "noUnusedLocals": true,
18+ "noUnusedParameters": true
19+ },
20+ "include": ["vite.config.ts"]
21+}
vite.config.tsadded+11−0View file
@@ -0,0 +1,11 @@
1+import { defineConfig } from 'vite'
2+import react from '@vitejs/plugin-react'
3+
4+export default defineConfig({
5+ plugins: [react()],
6+ // relative base so the same build works at the repository root and under a
7+ // sub-path on a static host
8+ base: './',
9+ // .m and .md files are pulled in with ?raw; tell vite they are assets
10+ assetsInclude: ['**/*.m'],
11+})