An introduction and static figures, and the MATLAB written plainly again
Two things.
An introduction before Chapter 1 defines Laplace's and Poisson's equations and
the Dirichlet, Neumann and Robin conditions, says what the tutorial measures,
and lays out the four domains the chapters work on. Every chapter gains one or
two drawn diagrams, addressed from the prose by `<!--fig:id-->` next to the
existing `<!--lab:id-->`: several chapters, 8 to 10 in particular, showed
nothing at all until the reader pressed Run.
And the MATLAB layer goes back to ordinary MATLAB. The drivers were full of
shapes chosen to keep numbl's compiler on its fast path: sweeps written as
recursions because a call inside a loop was compiled generically, histories
handed back as squares, a copy routed through a function call to avoid a slice
assignment. Those limitations were reported and fixed upstream in numbl 0.4.19,
so the sweeps are plain `for` loops again, roughly a hundred and fifty lines
lighter, and the comments no longer explain a compiler instead of a method.
Timings and every asserted number are unchanged; the requirement is recorded in
the README, and the one rule that still earns its keep, keeping the *timed*
loop type-stable, is in CLAUDE.md.
39 changed files+1731−407
CLAUDE.mdmodified+82−82View file
@@ -16,8 +16,10 @@ src/
1616 ui/ Lab (one runnable section), ChapterView, Controls, Markdown
1717 engine/ the numbl session wrapper
1818 plot/ a small hand-rolled SVG layer, plus a canvas field plot
19+ figures/ static diagrams, one .tsx per part, keyed in figures/index.ts
1920 matlab/ the shared .m library, shown to the reader on its own page
2021 chapters/
22+ intro/ the introduction: front matter, no labs
2123 onedim/driver.m the driver for chapters 1-4, plus their shared panels
2224 twodim/driver.m the driver for chapters 5-7, plus their shared panels
2325 chNN/ per-chapter text.md, editable .m script, chapter.ts
@@ -35,105 +37,82 @@ become **local functions of the driver** and the driver can call them. That is t
3537 whole of the contract; each driver's header comment states which functions its
3638 script must define.
3739
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.
40+## Working with numbl
41+
42+**This tutorial needs numbl 0.4.19 or later.** It runs correctly on 0.4.18 and
43+is drastically slower there: the L-shaped sweep in Chapter 9 takes about three
44+minutes instead of six seconds, and Chapter 8 about six seconds instead of a
45+fifth of one.
46+
47+The reason is worth recording. Most of what this section used to say was a list
48+of ways to fall off numbl's JIT fast path, found by profiling this tutorial: a
49+`switch` anywhere in a call chain, a struct parameter, `A(k,:) = row` with a
50+loop index, a bare `tic;`, and the expensive one: the rule that nothing nested
51+under a loop the JIT had declined was compiled either, so a single unsupported
52+line could put a whole solver on the interpreter. They were fixed
53+upstream in 0.4.19, and the drivers here were rewritten back into ordinary
54+MATLAB once they were: the sweeps are plain `for` loops again. If you find a
55+comment describing a recursion-instead-of-a-loop workaround, it is stale; the
56+git history has the old shape.
57+
58+What is left is worth knowing.
59+
60+### Keep the *timed* loop type-stable
61+
62+numbl specialises a loop on the types flowing through it. A loop whose iterate
63+is assigned from several different functions, or whose body indexes a cell
64+array, is not something it can specialise, and it runs interpreted. That is
65+usually harmless, since the expensive callee inside is still compiled, but when
66+the loop being measured *is* the thing under test the overhead lands in the
67+number on screen. Two consequences visible in `ch08/iterative.m`:
68+
69+- each iterative method has its own loop rather than one loop with a switch
70+ inside it, and multigrid keeps its grid hierarchy in one flat vector rather
71+ than a cell array;
72+- residual and error histories accumulate as squares and are square-rooted
73+ once at the end. A `sqrt` inside the loop costs about 60 per cent on the
74+ sweep, because sqrt may return a complex value and the history array is then
75+ no longer a plain real array.
76+
77+Neither is a correctness matter. Check with `%!numbl:assert_jit` inside a loop
78+body (it throws if that unit ran interpreted) and `NUMBL_JIT_LOG=1`, which
79+prints every JIT decline with its reason.
10480
10581 ### `jsondecode` of a JSON array of arrays gives a cell of column vectors
10682
107-MATLAB gives a matrix when the rows are the same length. The L-shaped driver unpacks
108-`p.nsm` by hand.
83+MATLAB gives a matrix when the rows are the same length. The L-shaped driver
84+unpacks `p.nsm` by hand.
10985
11086 ### `jsonencode` refuses complex values, and flattens single-row matrices
11187
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.
88+`jsonencode([1 2 3])` is `[1,2,3]`, not `[[1,2,3]]`, so anything matrix-shaped
89+goes through `lp_rows`, which converts it to a cell array of rows first. NaN and
90+Inf cross as `null`, which is why the TypeScript side uses `Num = number | null`
91+everywhere and every plot helper skips nulls. A value whose imaginary part is
92+entirely zero counts as real and encodes fine.
11693
11794 ### Backslash densifies sparse matrices
11895
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.
96+`A \ b` converts a sparse `A` to dense before calling LAPACK, so there is no
97+sparse direct solver here. That is why Chapter 5 stops at $N = 48$ and why
98+Chapters 6 to 10 are about methods that do not need one. The text says so where
99+it matters.
122100
123101 ### Scalar loops are often faster than vectorised slices
124102
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.
103+A hand-written double loop over a stencil beats the equivalent
104+`A(2:end-1,...)` expression by about a factor of three, because the vectorised
105+form allocates several temporaries per sweep. The kernels in
106+`ch08/iterative.m` are written as loops for that reason, which also happens to
107+make them readable.
129108
130109 ### The browser is slower than the CLI
131110
132111 The CLI uses a native LAPACK addon; the browser falls back to a TypeScript
133112 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.
113+noticeably slower in the page than in `npm run test:matlab`. Nothing in the
114+tutorial depends on the absolute numbers, but do not quote a CLI timing as if it
115+were what a reader will see.
137116
138117 ## Measuring time honestly
139118
@@ -160,6 +139,27 @@ and checks that nothing throws and that the figures contain data. It cannot judg
160139 whether a picture is any good; `--shots` writes full-page screenshots into `tmp/` so
161140 that a human, or an agent with vision, can look.
162141
142+## Static figures
143+
144+`src/figures/` holds the diagrams that are drawn rather than computed: the problem
145+statement, the stencils, the roundoff-floor V, the V-cycle, the re-entrant corner,
146+the close-evaluation spike. They are addressed from the prose by `<!--fig:id-->`, the
147+same way labs are addressed by `<!--lab:id-->`, and resolved through the registry in
148+`figures/index.ts`. They exist because several chapters, Chapters 8 to 10 in
149+particular, show nothing at all until the reader presses Run.
150+
151+They are plain inline SVG with no data behind them, written against a 716-unit-wide
152+viewBox and `width="100%"`, using `plot/palette.ts` for colour so that a curve
153+labelled `fd2` is the same blue everywhere. Two habits are worth keeping. Text is
154+placed by hand, so a caption or a legend that fits at one size may collide at
155+another: render and look. And the introduction's accuracy-against-cost picture is a
156+schematic, which its caption says out loud; do not let it drift into looking like
157+measured data.
158+
159+To look at them, point a headless browser at each chapter and screenshot
160+`figure.static-figure`; `npm run test:browser --shots` writes whole pages, which is
161+coarser but needs no extra script.
162+
163163 ## Plotting
164164
165165 `src/plot/` is a small SVG layer with no chart library. `palette.ts` records which
README.mdmodified+19−1View file
@@ -12,6 +12,12 @@ point to full precision. What the tutorial then measures is not whether a method
1212 works but what accuracy it delivers for what cost, and under which conditions that
1313 answer changes.
1414
15+An introduction, before Chapter 1, defines Laplace's and Poisson's equations and the
16+Dirichlet, Neumann and Robin boundary conditions, states what the tutorial measures,
17+and lays out the four domains the chapters work on. Every chapter also carries one or
18+two static diagrams, drawn as inline SVG, so that a page says something before its
19+labs are run.
20+
1521 ## The problems
1622
1723 In one dimension, Laplace's equation is $u'' = 0$ and its solution is a straight
@@ -28,6 +34,10 @@ error of the method.
2834
2935 ## The chapters
3036
37+**Introduction** — Laplace's and Poisson's equations, the Dirichlet, Neumann and
38+Robin boundary conditions, what accuracy-against-cost means here, and the four
39+domains the chapters work on.
40+
3141 **One dimension**
3242
3343 1. **Ground truth by construction** — manufactured solutions, the three-point
@@ -74,6 +84,12 @@ runs from inside the editor, and "reset" restores the original.
7484 Some labs run as soon as the chapter opens; the heavier ones wait for the run
7585 button, which the caption says.
7686
87+## Requirements
88+
89+numbl **0.4.19 or later**. The tutorial runs on 0.4.18 but is far slower there:
90+several JIT limitations that this tutorial's profiling turned up were fixed in
91+0.4.19, and the MATLAB here is written in the plain form those fixes allow.
92+
7793 ## Development
7894
7995 ```
@@ -84,7 +100,9 @@ npm run test:matlab # the .m layer, through the numbl CLI (needs a numbl chec
84100 npm run test:browser # the built app in headless Chrome
85101 ```
86102
87-`test:matlab` needs a numbl clone; set `NUMBL_DIR` if it is not at `~/src/numbl`.
103+`test:matlab` needs a numbl clone; set `NUMBL_DIR` if it is not at `~/src/numbl`,
104+and build its CLI (`npm run build:cli` in that checkout) after updating it — the
105+tests run `dist-cli/cli.js`, not the TypeScript sources.
88106 Every assertion in it is a number claimed somewhere in the prose, so a failure means
89107 either the code or the text is wrong. `test:browser` needs `npm run build` first, and
90108 takes `--shots` to write a full-page screenshot of each chapter into `tmp/`.
package-lock.jsonmodified+4−4View file
@@ -17,7 +17,7 @@
1717 "codemirror": "^6.0.2",
1818 "katex": "^0.16.11",
1919 "marked": "^15.0.0",
20- "numbl": "^0.4.18",
20+ "numbl": "^0.4.19",
2121 "react": "^19.2.0",
2222 "react-dom": "^19.2.0"
2323 },
@@ -3733,9 +3733,9 @@
37333733 }
37343734 },
37353735 "node_modules/numbl": {
3736- "version": "0.4.18",
3737- "resolved": "https://registry.npmjs.org/numbl/-/numbl-0.4.18.tgz",
3738- "integrity": "sha512-ZFmkwP0jZvXJSNUQbIDSB4/QeKtkXDWlP36oFlyH1F+DiHrovjBcRQBoYb1Qw8e6z/XSaPxtq183hWm6Rd4+WQ==",
3736+ "version": "0.4.19",
3737+ "resolved": "https://registry.npmjs.org/numbl/-/numbl-0.4.19.tgz",
3738+ "integrity": "sha512-wfTAhrdOSuzUoiglglTwykwEguU0Vl5qBJNhj2hyP8FzHACd/w1qX69OhEC7sbyE+gGEH9gjPjFsBkoi9UTGqw==",
37393739 "hasInstallScript": true,
37403740 "license": "Apache-2.0",
37413741 "dependencies": {
package.jsonmodified+1−1View file
@@ -22,7 +22,7 @@
2222 "codemirror": "^6.0.2",
2323 "katex": "^0.16.11",
2424 "marked": "^15.0.0",
25- "numbl": "^0.4.18",
25+ "numbl": "^0.4.19",
2626 "react": "^19.2.0",
2727 "react-dom": "^19.2.0"
2828 },
scripts/matlab-test.mjsmodified+3−3View file
@@ -166,15 +166,15 @@ console.log('\nChapter 8: iterative solvers')
166166 nit: 60,
167167 tol: 1e-11,
168168 })
169- const mg = out.R2[4].filter((v) => v != null)
169+ const mg = out.R[4].filter((v) => v != null)
170170 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)
171+ const mgErr = out.E[4].filter((v) => v != null)
172172 check(
173173 'and it lands on the discretisation error',
174174 near(mgErr[mgErr.length - 1] / out.discErr, 1, 0.02),
175175 `${mgErr[mgErr.length - 1].toExponential(3)} vs ${out.discErr.toExponential(3)}`,
176176 )
177- const jac = out.R2[0].filter((v) => v != null).map(Math.sqrt)
177+ const jac = out.R[0].filter((v) => v != null)
178178 check('Jacobi is nowhere near after 60 sweeps', jac[jac.length - 1] > 1e-2, jac[jac.length - 1].toExponential(2))
179179 }
180180
src/App.tsxmodified+5−2View file
@@ -3,7 +3,10 @@ import ChapterView from './ui/ChapterView.tsx'
33 import LibraryView from './ui/LibraryView.tsx'
44 import { CHAPTERS } from './chapters/index.ts'
55
6-const PARTS = ['One dimension', 'Two dimensions', 'Beyond the square'] as const
6+const PARTS = ['Introduction', 'One dimension', 'Two dimensions', 'Beyond the square'] as const
7+
8+/** The introduction is front matter, so the chapter numbers start after it. */
9+const NUMBERED = CHAPTERS.filter((c) => c.part !== 'Introduction')
710
811 function readHash(): string {
912 const h = window.location.hash.replace(/^#\/?/, '')
@@ -40,7 +43,7 @@ export default function App() {
4043 {CHAPTERS.filter((c) => c.part === part).map((c) => (
4144 <li key={c.id}>
4245 <button className={c.id === route ? 'nav-item active' : 'nav-item'} onClick={() => go(c.id)}>
43- <span className="nav-num">{CHAPTERS.indexOf(c) + 1}</span>
46+ <span className="nav-num">{NUMBERED.indexOf(c) + 1 || ''}</span>
4447 {c.nav}
4548 </button>
4649 </li>
src/chapters/ch01/text.mdmodified+8−2View file
@@ -1,8 +1,10 @@
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,
1+Part I works in one dimension, on Poisson's equation
22
33 $$-u''(x) = f(x), \qquad 0 < x < 1, \qquad u(0) = \alpha, \quad u(1) = \beta.$$
44
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.
5+As the introduction noted, Laplace's equation in one dimension is $u'' = 0$ and has nothing left to compute, so Poisson's is the model problem here. 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. From Chapter 5 on we set $f = 0$ and let the second dimension supply the difficulty.
6+
7+Before any of that can be measured, we need something to measure against, and that is where this chapter starts.
68
79 ### Knowing the answer in advance
810
@@ -14,6 +16,8 @@ $$f := -u'', \qquad \alpha := u(0), \qquad \beta := u(1).$$
1416
1517 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.
1618
19+<!--fig:manufactured-->
20+
1721 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.
1822
1923 ### The three-point scheme
@@ -26,6 +30,8 @@ A Taylor expansion of the right-hand side gives $u''(x_i) + \frac{h^2}{12}u''''(
2630
2731 $$\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}.$$
2832
33+<!--fig:three-point-->
34+
2935 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.
3036
3137 <!--lab:first-->
src/chapters/ch02/text.mdmodified+2−0View file
@@ -16,6 +16,8 @@ $$N_{\text{best}} \;\sim\; \varepsilon^{-1/4} \approx 10^4, \qquad \text{error}_
1616
1717 About half the digits, and no more, whatever the machine and however long one is willing to wait.
1818
19+<!--fig:roundoff-floor-->
20+
1921 <!--lab:floor-->
2022
2123 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.
src/chapters/ch03/text.mdmodified+2−0View file
@@ -12,6 +12,8 @@ $$\frac{\delta^2 u_i}{h^2} \;=\; -f_i - \frac{h^2}{12}\,\frac{\delta^2 f_i}{h^2}
1212
1313 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.
1414
15+<!--fig:compact-->
16+
1517 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.
1618
1719 <!--lab:orders-->
src/chapters/ch04/text.mdmodified+2−0View file
@@ -14,6 +14,8 @@ What is bought is the approximation power of polynomials of high degree. For $u$
1414
1515 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.
1616
17+<!--fig:chebyshev-points-->
18+
1719 <!--lab:single-->
1820
1921 ### Where the curve stops
src/chapters/ch05/text.mdmodified+2−0View file
@@ -25,6 +25,8 @@ $$\frac{u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1} - 4u_{i,j}}{h^2} \;=\; 0,$
2525
2626 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.
2727
28+<!--fig:five-point-->
29+
2830 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.
2931
3032 <!--lab:solve-->
src/chapters/ch06/text.mdmodified+2−0View file
@@ -12,6 +12,8 @@ $$A\,\big(v_k \otimes v_l\big) \;=\; (\lambda_k + \lambda_l)\,\big(v_k \otimes v
1212
1313 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.
1414
15+<!--fig:fast-poisson-->
16+
1517 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.
1618
1719 <!--lab:compare-->
src/chapters/ch07/text.mdmodified+2−0View file
@@ -14,6 +14,8 @@ $$20 - 8\cos a - 8\cos b - 4\cos a \cos b,$$
1414
1515 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.
1616
17+<!--fig:nine-point-->
18+
1719 ### Chebyshev in two dimensions
1820
1921 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,
src/chapters/ch08/driver.mmodified+9−21View file
@@ -5,11 +5,9 @@
55 % names = method_names()
66 % H = run_iter(mid, N, U0, UE, nit, tol)
77 %
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.
8+% where H has, in its three rows and one column per iteration, the relative
9+% residual, the maximum error, and the elapsed milliseconds. Columns past the
10+% iteration the method stopped at are NaN.
1311 %
1412 % and the driver supplies the starting array U0, which holds the exact
1513 % boundary data on the boundary and zeros inside, and the exact solution UE for
@@ -40,39 +38,29 @@ ms_direct = toc(t0) * 1000;
4038 disc_err = max(max(abs(Ud - UE)));
4139
4240 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
41+tol = p.tol;
4542 R = NaN(nm, nit+1);
4643 E = NaN(nm, nit+1);
4744 T = NaN(nm, nit+1);
4845 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
5246 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);
47+ R(k, :) = H(1, :);
48+ E(k, :) = H(2, :);
49+ T(k, :) = H(3, :);
5650 end
5751
5852 out = struct();
5953 out.names = names;
6054 out.N = N;
6155 out.its = 0:nit;
62-out.R2 = lp_rows(R);
63-out.E2 = lp_rows(E);
56+out.R = lp_rows(R);
57+out.E = lp_rows(E);
6458 out.T = lp_rows(T);
6559 out.discErr = disc_err;
6660 out.directMs = ms_direct;
6761 out.unknowns = (N-1)^2;
6862 lp_write(out);
6963
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-
7664 function U = drv_fast(U, N)
7765 n = N - 1;
7866 B = zeros(n, n);
src/chapters/ch08/iterative.mmodified+35−52View file
@@ -15,10 +15,10 @@
1515 %
1616 % A note on the shape of this file. Each method has its own loop rather than
1717 % 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.
18+% vector rather than a cell array of grids. Both keep the *timed* loop
19+% type-stable for numbl's compiler: an iterate assigned from several different
20+% functions, or a body that indexes a cell, is not something it can specialise,
21+% and here that loop is the measurement.
2222
2323 function names = method_names()
2424 names = {'Jacobi', 'Gauss-Seidel', 'SOR', 'conjugate gradients', 'multigrid V-cycles'};
@@ -47,6 +47,7 @@ function H = loop_jacobi(N, U, UE, nit, tol)
4747 res = NaN(1, nit+1);
4848 err = NaN(1, nit+1);
4949 tms = NaN(1, nit+1);
50+nf = 1;
5051 ss0 = res_sumsq(U, N);
5152 res(1) = ss0;
5253 err(1) = maxsq(U - UE);
@@ -60,15 +61,12 @@ for it = 1:nit
6061 res(it+1) = res_sumsq(U, N);
6162 err(it+1) = maxsq(U - UE);
6263 tms(it+1) = tacc * 1000;
64+ nf = it + 1;
6365 if res(it+1) < stop
6466 break;
6567 end
6668 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];
69+H = hist_pack(res, err, tms, ss0, nf);
7270 end
7371
7472 function H = loop_sor(N, U, UE, nit, tol, omega)
@@ -80,6 +78,7 @@ function H = loop_sor(N, U, UE, nit, tol, omega)
8078 res = NaN(1, nit+1);
8179 err = NaN(1, nit+1);
8280 tms = NaN(1, nit+1);
81+nf = 1;
8382 ss0 = res_sumsq(U, N);
8483 res(1) = ss0;
8584 err(1) = maxsq(U - UE);
@@ -93,15 +92,12 @@ for it = 1:nit
9392 res(it+1) = res_sumsq(U, N);
9493 err(it+1) = maxsq(U - UE);
9594 tms(it+1) = tacc * 1000;
95+ nf = it + 1;
9696 if res(it+1) < stop
9797 break;
9898 end
9999 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];
100+H = hist_pack(res, err, tms, ss0, nf);
105101 end
106102
107103 function H = loop_cg(N, U, UE, nit, tol)
@@ -113,6 +109,7 @@ function H = loop_cg(N, U, UE, nit, tol)
113109 res = NaN(1, nit+1);
114110 err = NaN(1, nit+1);
115111 tms = NaN(1, nit+1);
112+nf = 1;
116113 R = neg_residual(U, N); % the right-hand side is zero at interior points
117114 P = R;
118115 rs = dot2(R, R);
@@ -135,15 +132,12 @@ for it = 1:nit
135132 res(it+1) = rsn;
136133 err(it+1) = maxsq(U - UE);
137134 tms(it+1) = tacc * 1000;
135+ nf = it + 1;
138136 if res(it+1) < stop
139137 break;
140138 end
141139 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];
140+H = hist_pack(res, err, tms, ss0, nf);
147141 end
148142
149143 function H = loop_mg(N, U, UE, nit, tol)
@@ -177,6 +171,7 @@ WE = flatten_into(zeros((N+1)^2, 1), UE, 0, N);
177171 res = NaN(1, nit+1);
178172 err = NaN(1, nit+1);
179173 tms = NaN(1, nit+1);
174+nf = 1;
180175 ss0 = res_sumsq_flat(W, 0, tot, N);
181176 res(1) = ss0;
182177 err(1) = maxsq_flat(W, WE, N);
@@ -190,15 +185,12 @@ for it = 1:nit
190185 res(it+1) = res_sumsq_flat(W, 0, tot, N);
191186 err(it+1) = maxsq_flat(W, WE, N);
192187 tms(it+1) = tacc * 1000;
188+ nf = it + 1;
193189 if res(it+1) < stop
194190 break;
195191 end
196192 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];
193+H = hist_pack(res, err, tms, ss0, nf);
202194 end
203195
204196 % ---------------------------------------------------------------- kernels ---
@@ -258,6 +250,19 @@ for j = 2:N
258250 end
259251 end
260252
253+function H = hist_pack(res, err, tms, ss0, nf)
254+% Both histories are accumulated as squares and rooted here, over the first nf
255+% columns — the part the iteration actually filled; the rest stays NaN, which
256+% is what the page draws as a gap. Rooting inside the loop instead costs about
257+% 60 per cent on the sweep itself, because sqrt may return a complex value and
258+% the history array is then no longer a plain real array to numbl.
259+r = res;
260+e = err;
261+r(1:nf) = sqrt(res(1:nf) / ss0);
262+e(1:nf) = sqrt(err(1:nf));
263+H = [r; e; tms];
264+end
265+
261266 function m = maxsq(E)
262267 m = max(max(E .* E));
263268 end
@@ -310,24 +315,10 @@ end
310315
311316 function W = mg_down(W, off, offc, tot, n, nu)
312317 % 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.
318+% iterate at zero. The residual is formed inline, one 9-point full-weighting
319+% stencil per coarse point, rather than materialised on the fine grid first.
320+W = gs_flat(W, off, tot, n, nu);
321321 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-
331322 nc = n / 2;
332323 mc = nc + 1;
333324 for jc = 2:nc
@@ -387,20 +378,12 @@ for j = 1:m
387378 W(b+i) = W(b+i) + Z(i,j);
388379 end
389380 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
381+W = gs_flat(W, off, tot, n, nu);
398382 end
399383
400384 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.
385+% Gauss-Seidel on one level of the flat hierarchy: nu sweeps over the interior
386+% of the level at `off`, reading its right-hand side from the second half of W.
404387 m = n + 1;
405388 for s = 1:nu
406389 for j = 2:n
src/chapters/ch08/panels.tsxmodified+5−10View file
@@ -12,19 +12,14 @@ export interface IterOut {
1212 names: string[]
1313 N: number
1414 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[][]
15+ R: Num[][]
16+ E: Num[][]
1917 T: Num[][]
2018 discErr: number
2119 directMs: number
2220 unknowns: number
2321 }
2422
25-const rootOf = (rows: Num[][]): Num[][] =>
26- rows.map((r) => r.map((v) => (v == null || !isFinite(v) || v < 0 ? null : Math.sqrt(v))))
27-
2823 // one palette slot per method, in the order the script lists them
2924 const SLOTS: MethodKey[] = ['fd2', 'fd4', 'fast', 'spectral', 'iterative']
3025
@@ -49,8 +44,8 @@ function lastIndex(row: Num[]): number {
4944
5045 export function IterPanels({ data }: { data: IterOut }) {
5146 const its = data.its
52- const R = rootOf(data.R2)
53- const E = rootOf(data.E2)
47+ const R = data.R
48+ const E = data.E
5449 const xd: [number, number] = [0, its[its.length - 1]]
5550 const rd = lowHigh(R)
5651 const ed = lowHigh([...E, [data.discErr]])
@@ -167,7 +162,7 @@ export function IterPanels({ data }: { data: IterOut }) {
167162 }
168163
169164 export function iterSummary(d: IterOut) {
170- const mg = d.E2[d.E2.length - 1]
165+ const mg = d.E[d.E.length - 1]
171166 const i = lastIndex(mg)
172167 return [
173168 { label: 'unknowns', value: String(d.unknowns) },
src/chapters/ch08/text.mdmodified+4−0View file
@@ -16,8 +16,12 @@ which turns $O(N^2)$ sweeps into $O(N)$. That is a genuine change of order, and
1616
1717 **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.
1818
19+<!--fig:smoothing-->
20+
1921 **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.
2022
23+<!--fig:vcycle-->
24+
2125 <!--lab:history-->
2226
2327 ### Reading the three figures
src/chapters/ch09/driver.mmodified+18−44View file
@@ -55,10 +55,10 @@ else
5555 py = Y(inside);
5656 ve = lp_lexact(kind, px, py);
5757
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
58+ % Each method gets its own list of resolutions: intervals per unit length
59+ % for the difference scheme, basis terms for the expansion. numbl's
60+ % jsondecode turns a JSON array of arrays into a cell of column vectors
61+ % where MATLAB would give a matrix, so unpack it by hand.
6262 NSc = p.nsm;
6363 nn = numel(NSc{1});
6464 NS = zeros(nm, nn);
@@ -66,15 +66,21 @@ else
6666 row = NSc{k};
6767 NS(k, :) = row(:)';
6868 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
69+ E = NaN(nm, nn);
70+ T = NaN(nm, nn);
71+ Q = NaN(nm, nn);
7772 for k = 1:nm
73+ vw = solve_l(k, NS(k, 1), g, px(1), py(1)); % warm-up
74+ for j = 1:nn
75+ n = NS(k, j);
76+ if ~isfinite(n) || n <= 0
77+ continue;
78+ end
79+ t0 = tic;
80+ v = solve_l(k, n, g, px, py);
81+ T(k, j) = toc(t0) * 1000;
82+ E(k, j) = max(abs(v - ve));
83+ end
7884 Q(k, :) = lp_orders(E(k, :), NS(k, :));
7985 end
8086 out.ns = NS(1, :);
@@ -85,35 +91,3 @@ else
8591 end
8692
8793 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.mmodified+1−5View file
@@ -52,11 +52,7 @@ U = zeros(m, m);
5252 U(~unknown) = g(X(~unknown), Y(~unknown)); % boundary values
5353 U(~inL) = 0; % outside: never read
5454
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);
55+S = cg_solve(U, mask, m, 1e-9, 4000); % U holds the boundary data, S the solution
6056
6157 % the requested points are grid points, so sampling is a lookup
6258 v = zeros(numel(px), 1);
src/chapters/ch09/text.mdmodified+2−0View file
@@ -14,6 +14,8 @@ $$s_k(r,\theta) = r^{2k/3}\sin\frac{2k\theta}{3}, \qquad k = 1, 2, 3, \dots$$
1414
1515 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.
1616
17+<!--fig:reentrant-corner-->
18+
1719 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.
1820
1921 Three ground truths are available in the labs, so that the effect can be isolated:
src/chapters/ch10/driver.mmodified+22−54View file
@@ -74,15 +74,29 @@ else
7474
7575 ns = p.ns(:)';
7676 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);
77+ E = NaN(nm, nn);
78+ T = NaN(nm, nn);
79+ Q = NaN(nm, nn);
80+ EN = NaN(nm, nn);
8281 for k = 1:nm
83- E(k, :) = ET(3*k-2, :);
84- EN(k, :) = ET(3*k-1, :);
85- T(k, :) = ET(3*k, :);
82+ vw = solve_bie(k, 16, a, nlobe, p.par, px(1), py(1)); % warm-up
83+ for j = 1:nn
84+ n = ns(j);
85+ t0 = tic;
86+ v = solve_bie(k, n, a, nlobe, p.par, px, py);
87+ el = toc(t0);
88+ rep = 0;
89+ while el < 0.005 && rep < 4
90+ rep = rep + 1;
91+ t0 = tic;
92+ v = solve_bie(k, n, a, nlobe, p.par, px, py);
93+ el = min(el, toc(t0));
94+ end
95+ q = solve_bie(k, n, a, nlobe, p.par, qx, qy);
96+ E(k, j) = max(abs(v - ve));
97+ EN(k, j) = max(abs(q - qe));
98+ T(k, j) = el * 1000;
99+ end
86100 Q(k, :) = lp_orders(E(k, :), ns);
87101 end
88102 out.Enear = lp_rows(EN);
@@ -93,49 +107,3 @@ else
93107 end
94108
95109 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.mdmodified+4−0View file
@@ -4,6 +4,8 @@ $$r(t) = 1 + a\cos(n_{\text{lobe}}\, t),$$
44
55 which no tensor-product grid fits, and the unknown lives on the boundary alone.
66
7+<!--fig:area-boundary-->
8+
79 ### Reformulating on the boundary
810
911 Look for the solution as the field of a dipole layer spread over the curve with density $\sigma$:
@@ -44,6 +46,8 @@ The first row is the point of the whole tutorial in one line: full double precis
4446
4547 ### Three things this does not do
4648
49+<!--fig:close-evaluation-->
50+
4751 **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.
4852
4953 **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.
src/chapters/index.tsmodified+2−1View file
@@ -1,4 +1,5 @@
11 import type { Chapter } from '../ui/types.ts'
2+import { intro } from './intro/chapter.ts'
23 import { ch01 } from './ch01/chapter.ts'
34 import { ch02 } from './ch02/chapter.ts'
45 import { ch03 } from './ch03/chapter.ts'
@@ -10,4 +11,4 @@ import { ch08 } from './ch08/chapter.ts'
1011 import { ch09 } from './ch09/chapter.ts'
1112 import { ch10 } from './ch10/chapter.ts'
1213
13-export const CHAPTERS: Chapter[] = [ch01, ch02, ch03, ch04, ch05, ch06, ch07, ch08, ch09, ch10]
14+export const CHAPTERS: Chapter[] = [intro, ch01, ch02, ch03, ch04, ch05, ch06, ch07, ch08, ch09, ch10]
src/chapters/intro/chapter.tsadded+12−0View file
@@ -0,0 +1,12 @@
1+import type { Chapter } from '../../ui/types.ts'
2+import text from './text.md?raw'
3+
4+export const intro: Chapter = {
5+ id: 'introduction',
6+ nav: 'Introduction',
7+ part: 'Introduction',
8+ title: 'The Dirichlet problem',
9+ blurb: 'What the problem is, what the tutorial measures, and how the ten chapters are arranged.',
10+ text,
11+ labs: [],
12+}
src/chapters/intro/text.mdadded+43−0View file
@@ -0,0 +1,43 @@
1+Many of the steady states of classical physics satisfy one equation. The temperature in a conducting plate that has stopped changing, the electrostatic potential in a region free of charge, the velocity potential of an incompressible irrotational flow, and the displacement of a stretched membrane at rest are all solutions of Laplace's equation,
2+
3+$$\Delta u \;=\; \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \;=\; 0 .$$
4+
5+A function satisfying it is called *harmonic*. If the interior carries a source, a heat source, a charge density, a load on the membrane, the equation gains a right-hand side and becomes Poisson's equation, $-\Delta u = f$. The two differ only in that right-hand side, and every method in this tutorial applies to both, though not always at the same order of accuracy.
6+
7+The equation by itself does not determine $u$. Adding any harmonic function to a solution gives another solution, so a condition on the boundary $\partial\Omega$ of the region $\Omega$ is needed. Three kinds are standard. A **Dirichlet** condition prescribes $u = g$ on the boundary: the plate is held at a known temperature along its edge. A **Neumann** condition prescribes the outward normal derivative $\partial u/\partial n = g$: the flux through the edge is known instead. A **Robin**, or mixed, condition prescribes a combination of the two. We say the *Dirichlet problem* is the first of these:
8+
9+$$\Delta u = 0 \ \ \text{in } \Omega, \qquad u = g \ \ \text{on } \partial\Omega,$$
10+
11+or the same with $-\Delta u = f$ inside. It is well posed: for a reasonable domain and continuous data a solution exists, it is unique, and it depends continuously on $g$. Neumann and mixed conditions are not treated anywhere in these chapters, and neither are variable coefficients, nonlinearity, or time dependence. What is left is the simplest elliptic boundary value problem there is, which is exactly why it is worth ten chapters: nearly every technique used on harder elliptic problems can be introduced here, where the answer is known and nothing else is in the way.
12+
13+<!--fig:dirichlet-problem-->
14+
15+### What is being measured
16+
17+A tutorial on numerical methods has to decide what counts as success. Showing that a method converges is not difficult and not very informative, since almost anything consistent and stable converges. The question a user actually faces is narrower: for the accuracy I need, on the problem I have, which method gets there in the least time, and under what conditions does that answer change.
18+
19+To ask that question honestly one needs an error, and an error needs a truth to measure against. Reference solutions computed at higher resolution by the same family of methods are the usual substitute, and they are least reliable exactly where the interesting behaviour is. Here we avoid them entirely. Every problem in this tutorial has an exact solution known in closed form, either manufactured by choosing $u$ first and defining the data from it, or taken from the supply of harmonic functions that can be written down. The error of every computation is therefore known at every point, to full precision, at every resolution.
20+
21+What the chapters then compare is accuracy against cost, with the cost measured as time for one solve inside the numerical layer. The spread turns out to be large. Across the ten chapters, methods that are reasonable choices for the same problem differ by roughly fourteen orders of magnitude in the accuracy they deliver for comparable work.
22+
23+<!--fig:accuracy-cost-->
24+
25+Two things drive that spread, and they recur in every chapter. The first is the order of the method, which sets the slope of the curve rather than its position, and therefore decides who wins only once the accuracy demand is serious. The second is structure: a method that exploits a special feature of the problem, a rectangle, a separable operator, an analytic solution, a smooth boundary, can be faster than a general method by factors that no amount of tuning recovers. The conditions attached to each such method are as much the subject of this tutorial as the speedups are.
26+
27+### The domains
28+
29+The geometry is what removes those special features, one at a time, so the chapters are organised by domain rather than by method.
30+
31+<!--fig:domains-->
32+
33+**Part I, one dimension.** Poisson's equation on $[0,1]$. In one dimension Laplace's equation is $u'' = 0$, whose solution is the straight line through the two boundary values, so there is nothing to compute and the model problem has to be Poisson's. It has the same conditioning and the same convergence theory as the two-dimensional problem at a hundredth of the cost, which makes it the right place to develop the tools: manufactured solutions and the three-point scheme (Chapter 1), the roundoff floor and Richardson extrapolation (Chapter 2), a compact fourth-order scheme (Chapter 3), and Chebyshev collocation (Chapter 4).
34+
35+**Part II, the unit square.** Laplace's equation, with harmonic ground truth. The five-point Laplacian and the cost of a dense solve (Chapter 5), the same system diagonalised by the sine transform (Chapter 6), a sixth-order compact stencil for the same price (Chapter 7), and what to do when no direct solver is available (Chapter 8).
36+
37+**Part III, beyond the square.** An L-shaped domain, where a re-entrant corner puts a singularity into the solution and second order is lost (Chapter 9), and a smooth closed curve, where the unknown moves onto the boundary and the answer comes back at full double precision (Chapter 10).
38+
39+### How the pages work
40+
41+Each chapter is prose with runnable sections dropped into it. Every method is a short MATLAB-syntax script that you can read, edit, and re-run, with sliders for the resolution and for which exact solution is being asked for. The code runs in your browser through [numbl](https://numbl.org), a MATLAB-to-JavaScript compiler: there is no server, nothing to install, and nothing is sent anywhere. Heavier sections wait for a Run button rather than starting on sight.
42+
43+The timings you see are your own machine's, measured around the solve and nothing else. Any absolute number quoted in the prose was measured at the command line, where a native linear algebra library is available; in the browser the same code is a few times slower. Nothing in the tutorial depends on the absolute numbers, only on the ratios.
src/chapters/onedim/driver.mmodified+27−63View file
@@ -11,13 +11,10 @@
1111 % function handle giving the source term, so a method is free to evaluate it
1212 % anywhere.
1313 %
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.
14+% Timings come from tic/toc around the solve and nothing else. Each method gets
15+% a warm-up call before its first timed one, because the first call also pays
16+% for compiling the method, and anything under 5 ms is repeated with the best
17+% kept: a sub-millisecond solve otherwise measures the clock.
2118 p = lp_params();
2219 kind = p.kind;
2320 names = method_names();
@@ -47,15 +44,32 @@ if strcmp(p.mode, 'single')
4744 else
4845 ns = p.ns(:)';
4946 nn = numel(ns);
50- E = zeros(nm, nn);
51- T = zeros(nm, nn);
52- Q = zeros(nm, nn);
47+ E = NaN(nm, nn);
48+ T = NaN(nm, nn);
49+ Q = NaN(nm, nn);
5350 nmax = p.nmaxs(:)';
5451 timed = p.timed;
55- ET = drv_methods(1, NaN(2*nm, nn), ns, f, ub(1), ub(2), kind, timed, nmax, nm);
5652 for k = 1:nm
57- E(k, :) = ET(2*k-1, :);
58- T(k, :) = ET(2*k, :);
53+ [xw, uw] = solve1d(k, 8, f, ub(1), ub(2)); % warm-up
54+ for j = 1:nn
55+ N = ns(j);
56+ if N > nmax(k)
57+ continue;
58+ end
59+ t0 = tic;
60+ [x, u] = solve1d(k, N, f, ub(1), ub(2));
61+ el = toc(t0);
62+ rep = 0;
63+ while timed > 0 && el < 0.005 && rep < 4
64+ rep = rep + 1;
65+ t0 = tic;
66+ [x, u] = solve1d(k, N, f, ub(1), ub(2));
67+ el = min(el, toc(t0));
68+ end
69+ uex = lp_exact1d(kind, x(:));
70+ E(k, j) = max(abs(u(:) - uex(:)));
71+ T(k, j) = el * 1000;
72+ end
5973 Q(k, :) = lp_orders(E(k, :), ns);
6074 end
6175 out.ns = ns;
@@ -66,56 +80,6 @@ end
6680
6781 lp_write(out);
6882
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-
11983 function y = drv_rhs(kind, x)
12084 [~, y] = lp_exact1d(kind, x);
12185 end
src/chapters/twodim/driver.mmodified+23−51View file
@@ -42,15 +42,32 @@ if strcmp(p.mode, 'single')
4242 else
4343 ns = p.ns(:)';
4444 nn = numel(ns);
45- E = zeros(nm, nn);
46- T = zeros(nm, nn);
47- Q = zeros(nm, nn);
45+ E = NaN(nm, nn);
46+ T = NaN(nm, nn);
47+ Q = NaN(nm, nn);
4848 nmax = p.nmaxs(:)';
4949 timed = p.timed;
50- ET = drv_methods2(1, NaN(2*nm, nn), ns, g, kind, par, nmax, timed, nm);
5150 for k = 1:nm
52- E(k, :) = ET(2*k-1, :);
53- T(k, :) = ET(2*k, :);
51+ [xw, yw, Uw] = solve2d(k, 8, g); % warm-up
52+ for j = 1:nn
53+ N = ns(j);
54+ if N > nmax(k)
55+ continue;
56+ end
57+ t0 = tic;
58+ [x, y, U] = solve2d(k, N, g);
59+ el = toc(t0);
60+ rep = 0;
61+ while timed > 0 && el < 0.005 && rep < 4
62+ rep = rep + 1;
63+ t0 = tic;
64+ [x, y, U] = solve2d(k, N, g);
65+ el = min(el, toc(t0));
66+ end
67+ [X, Y] = meshgrid(x(:), y(:));
68+ E(k, j) = max(max(abs(U - lp_exact2d(kind, X, Y, par))));
69+ T(k, j) = el * 1000;
70+ end
5471 Q(k, :) = lp_orders(E(k, :), ns);
5572 end
5673 out.ns = ns;
@@ -60,48 +77,3 @@ else
6077 end
6178
6279 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/engine/runner.tsmodified+0−0View file
Binary file not shown.
src/figures/beyond.tsxadded+227−0View file
@@ -0,0 +1,227 @@
1+// Static figures for Chapters 9 and 10.
2+import { Arrowheads, Figure, Label, Node, arrow } from './svg.tsx'
3+import { ink, method, seqColor } from '../plot/palette.ts'
4+
5+/** The L-shaped domain, its re-entrant corner, and the profile of r^(2/3). */
6+export function ReentrantCorner() {
7+ const S = 62 // half-width of the square, in px
8+ const CX = 150
9+ const CY = 112
10+ return (
11+ <Figure
12+ width={716}
13+ height={236}
14+ caption="The re-entrant corner and what it does to the solution. A convex corner of angle π/2 gives r², which is harmless; the 3π/2 corner gives r²ᐟ³, whose gradient is infinite at the origin. No difference scheme recovers second order from that."
15+ >
16+ <Arrowheads id="rc" colors={[method.mps, ink.muted]} />
17+
18+ {/* the domain */}
19+ <path
20+ d={`M ${CX - S} ${CY - S} H ${CX + S} V ${CY} H ${CX} V ${CY + S} H ${CX - S} Z`}
21+ fill="#20242c"
22+ stroke={method.iterative}
23+ strokeWidth={2.2}
24+ strokeLinejoin="round"
25+ />
26+ <Label x={CX - 30} y={CY - 34} size={13} fill={ink.secondary}>
27+ Ω
28+ </Label>
29+ <Label x={CX} y={18} size={12.5} weight={600} fill={ink.primary}>
30+ the L-shaped domain
31+ </Label>
32+
33+ {/* the corner angle */}
34+ <path
35+ d={`M ${CX + 26} ${CY} A 26 26 0 1 0 ${CX} ${CY + 26}`}
36+ fill="none"
37+ stroke={method.mps}
38+ strokeWidth={1.6}
39+ markerEnd={arrow('rc', 0)}
40+ />
41+ <Label x={CX - 32} y={CY + 54} size={12} anchor="middle" fill={method.mps}>
42+ α = 3π/2
43+ </Label>
44+ <Node x={CX} y={CY} r={3.6} fill={method.mps} />
45+
46+ <Label x={CX} y={CY + S + 26} size={11.5} fill={ink.muted}>
47+ π/2 at every other corner
48+ </Label>
49+
50+ {/* the profile */}
51+ <g transform="translate(354 0)">
52+ <Label x={160} y={18} size={12.5} weight={600} fill={ink.primary}>
53+ radial profile at the corner
54+ </Label>
55+ <line x1={20} y1={186} x2={300} y2={186} stroke={ink.axis} strokeWidth={1.2} />
56+ <line x1={20} y1={186} x2={20} y2={38} stroke={ink.axis} strokeWidth={1.2} />
57+ <Label x={160} y={206} size={11.5} fill={ink.muted} italic>
58+ distance r from the corner
59+ </Label>
60+ {(
61+ [
62+ { p: 2 / 3, c: method.mps, l: 'r²ᐟ³ at the re-entrant corner' },
63+ { p: 2, c: ink.muted, l: 'r² at a convex corner' },
64+ ] as const
65+ ).map((s, k) => {
66+ const pts: string[] = []
67+ for (let i = 0; i <= 100; i++) {
68+ const t = i / 100
69+ pts.push(`${(20 + t * 280).toFixed(1)} ${(186 - 140 * Math.pow(t, s.p)).toFixed(1)}`)
70+ }
71+ return (
72+ <g key={k}>
73+ <path d={'M ' + pts.join(' L ')} fill="none" stroke={s.c} strokeWidth={2.2} />
74+ <line x1={34} y1={52 + k * 19} x2={52} y2={52 + k * 19} stroke={s.c} strokeWidth={2.2} />
75+ <Label x={60} y={52 + k * 19} size={11.5} anchor="start" fill={s.c}>
76+ {s.l}
77+ </Label>
78+ </g>
79+ )
80+ })}
81+
82+ </g>
83+ </Figure>
84+ )
85+}
86+
87+/** Discretising the area against discretising the boundary. */
88+export function AreaAgainstBoundary() {
89+ const blobPts = (cx: number, cy: number, rr: number, m: number) =>
90+ Array.from({ length: m }, (_, i) => {
91+ const t = (2 * Math.PI * i) / m
92+ const r = rr * (1 + 0.18 * Math.cos(3 * t))
93+ return { x: cx + r * Math.cos(t), y: cy + r * Math.sin(t) * 0.92, t }
94+ })
95+ const path = (cx: number, cy: number, rr: number) =>
96+ 'M ' + blobPts(cx, cy, rr, 160).map((p) => `${p.x.toFixed(1)} ${p.y.toFixed(1)}`).join(' L ') + ' Z'
97+
98+ const cx1 = 168
99+ const cx2 = 548
100+ const cy = 118
101+ const rr = 78
102+
103+ return (
104+ <Figure
105+ width={716}
106+ height={258}
107+ caption="Two ways to spend unknowns on the same problem. Filling the area needs a grid the curve does not fit; putting the unknown on the curve itself needs a few hundred points, and the equation is then satisfied exactly in the interior by construction."
108+ >
109+ {/* area discretisation */}
110+ <Label x={cx1} y={18} size={12.5} weight={600} fill={ink.primary}>
111+ unknowns in the area
112+ </Label>
113+ <clipPath id="ab-clip">
114+ <path d={path(cx1, cy, rr)} />
115+ </clipPath>
116+ <path d={path(cx1, cy, rr)} fill="#20242c" stroke={ink.axis} strokeWidth={1.6} />
117+ <g clipPath="url(#ab-clip)">
118+ {Array.from({ length: 19 }, (_, i) => (
119+ <line key={`v${i}`} x1={cx1 - 96 + i * 11} y1={cy - 92} x2={cx1 - 96 + i * 11} y2={cy + 96} stroke={ink.axis} strokeWidth={1} />
120+ ))}
121+ {Array.from({ length: 18 }, (_, j) => (
122+ <line key={`h${j}`} x1={cx1 - 100} y1={cy - 88 + j * 11} x2={cx1 + 100} y2={cy - 88 + j * 11} stroke={ink.axis} strokeWidth={1} />
123+ ))}
124+ {Array.from({ length: 19 }, (_, i) =>
125+ Array.from({ length: 17 }, (_, j) => (
126+ <Node key={`${i}-${j}`} x={cx1 - 96 + i * 11} y={cy - 88 + j * 11} r={1.9} fill={seqColor(0.55)} />
127+ )),
128+ )}
129+ </g>
130+ <Label x={cx1} y={cy + rr + 44} size={11.5} fill={ink.muted}>
131+ O(N²) unknowns, and the boundary falls between nodes
132+ </Label>
133+
134+ {/* boundary discretisation */}
135+ <Label x={cx2} y={18} size={12.5} weight={600} fill={ink.primary}>
136+ unknowns on the boundary
137+ </Label>
138+ <path d={path(cx2, cy, rr)} fill="#1e1e1c" stroke={method.bie} strokeWidth={2.2} />
139+ {blobPts(cx2, cy, rr, 36).map((p, i) => (
140+ <Node key={i} x={p.x} y={p.y} r={3} fill={method.bie} />
141+ ))}
142+ <Label x={cx2} y={cy} size={12} fill={ink.secondary}>
143+ Δu = 0 exactly
144+ </Label>
145+ <Label x={cx2} y={cy + rr + 44} size={11.5} fill={ink.muted}>
146+ O(n) unknowns, placed where the data is
147+ </Label>
148+ </Figure>
149+ )
150+}
151+
152+/** Why evaluation near the curve is the hard part. */
153+export function CloseEvaluation() {
154+ const cx = 200
155+ const cy = 116
156+ const rr = 86
157+ const pts = (m: number) =>
158+ Array.from({ length: m }, (_, i) => {
159+ const t = (2 * Math.PI * i) / m
160+ const r = rr * (1 + 0.16 * Math.cos(3 * t))
161+ return { x: cx + r * Math.cos(t), y: cy + r * Math.sin(t) * 0.94 }
162+ })
163+ const d = 'M ' + pts(180).map((p) => `${p.x.toFixed(1)} ${p.y.toFixed(1)}`).join(' L ') + ' Z'
164+ const far = { x: cx - 12, y: cy + 6 }
165+ const near = { x: cx + 74, y: cy - 26 }
166+
167+ // the integrand along the curve, seen from each target
168+ const K = (spread: number, amp: number) => {
169+ const out: string[] = []
170+ for (let i = 0; i <= 200; i++) {
171+ const t = -1 + (2 * i) / 200
172+ const v = amp / (1 + (t / spread) ** 2)
173+ out.push(`${(420 + (t + 1) * 140).toFixed(1)} ${(186 - v).toFixed(1)}`)
174+ }
175+ return 'M ' + out.join(' L ')
176+ }
177+
178+ return (
179+ <Figure
180+ width={716}
181+ height={256}
182+ caption="Close evaluation. The density on the curve is accurate in both cases; what fails is the quadrature used to evaluate the potential. Seen from a target well inside, the integrand is a gentle periodic function; seen from a target one grid spacing away, it is a spike the trapezoid rule cannot resolve."
183+ >
184+ <Arrowheads id="ce" colors={[ink.muted]} />
185+ <path d={d} fill="#20242c" stroke={method.bie} strokeWidth={2.2} />
186+ {pts(40).map((p, i) => (
187+ <Node key={i} x={p.x} y={p.y} r={2.6} fill={method.bie} />
188+ ))}
189+ <Node x={far.x} y={far.y} r={4.4} fill={method.spectral} />
190+ <Label x={far.x} y={far.y + 20} size={11.5} fill={method.spectral}>
191+ target well inside
192+ </Label>
193+ <Node x={near.x} y={near.y} r={4.4} fill={method.iterative} />
194+ <line x1={near.x + 6} y1={near.y - 4} x2={near.x + 40} y2={near.y - 28} stroke={method.iterative} strokeWidth={1} />
195+ <Label x={near.x + 44} y={near.y - 32} size={11.5} anchor="start" fill={method.iterative}>
196+ target one grid
197+ </Label>
198+ <Label x={near.x + 44} y={near.y - 17} size={11.5} anchor="start" fill={method.iterative}>
199+ spacing away
200+ </Label>
201+
202+ <line x1={332} y1={130} x2={382} y2={130} stroke={ink.muted} strokeWidth={1.3} markerEnd={arrow('ce', 0)} />
203+ <Label x={357} y={114} size={11} fill={ink.muted}>
204+ integrand
205+ </Label>
206+
207+ <line x1={420} y1={186} x2={700} y2={186} stroke={ink.axis} strokeWidth={1.2} />
208+ <Label x={560} y={206} size={11.5} fill={ink.muted} italic>
209+ arclength around the curve
210+ </Label>
211+ <path d={K(0.55, 42)} fill="none" stroke={method.spectral} strokeWidth={2.2} />
212+ <path d={K(0.035, 138)} fill="none" stroke={method.iterative} strokeWidth={2.2} />
213+ {Array.from({ length: 21 }, (_, i) => (
214+ <Node key={i} x={420 + i * 14} y={186} r={2.2} fill={ink.family} />
215+ ))}
216+ <line x1={424} y1={228} x2={442} y2={228} stroke={method.spectral} strokeWidth={2.2} />
217+ <Label x={450} y={228} size={11.5} anchor="start" fill={method.spectral}>
218+ from the far target: smooth
219+ </Label>
220+ <line x1={424} y1={246} x2={442} y2={246} stroke={method.iterative} strokeWidth={2.2} />
221+ <Label x={450} y={246} size={11.5} anchor="start" fill={method.iterative}>
222+ from the near target: a spike
223+ </Label>
224+
225+ </Figure>
226+ )
227+}
src/figures/index.tsadded+26−0View file
@@ -0,0 +1,26 @@
1+// Static figures, addressed from the prose by <!--fig:id--> the same way labs
2+// are addressed by <!--lab:id-->.
3+import type { ComponentType } from 'react'
4+import { AccuracyCost, DirichletProblem, Domains } from './intro.tsx'
5+import { ChebyshevPoints, CompactScheme, ManufacturedLoop, RoundoffFloor, ThreePointStencil } from './onedim.tsx'
6+import { FastPoissonSteps, FivePointStencil, NinePointStencil, SmoothingIdea, VCycle } from './twodim.tsx'
7+import { AreaAgainstBoundary, CloseEvaluation, ReentrantCorner } from './beyond.tsx'
8+
9+export const FIGURES: Record<string, ComponentType> = {
10+ 'dirichlet-problem': DirichletProblem,
11+ domains: Domains,
12+ 'accuracy-cost': AccuracyCost,
13+ 'three-point': ThreePointStencil,
14+ manufactured: ManufacturedLoop,
15+ 'roundoff-floor': RoundoffFloor,
16+ compact: CompactScheme,
17+ 'chebyshev-points': ChebyshevPoints,
18+ 'five-point': FivePointStencil,
19+ 'fast-poisson': FastPoissonSteps,
20+ 'nine-point': NinePointStencil,
21+ smoothing: SmoothingIdea,
22+ vcycle: VCycle,
23+ 'reentrant-corner': ReentrantCorner,
24+ 'area-boundary': AreaAgainstBoundary,
25+ 'close-evaluation': CloseEvaluation,
26+}
src/figures/intro.tsxadded+280−0View file
@@ -0,0 +1,280 @@
1+// Figures for the introduction: what the problem is, on what domains, and
2+// what the spread in accuracy-against-cost looks like before any of it is
3+// measured.
4+import type { ReactNode } from 'react'
5+import { Arrowheads, Figure, Label, Node, arrow } from './svg.tsx'
6+import { ink, method, seqColor } from '../plot/palette.ts'
7+
8+/** A closed wobbly curve, so that "a region" does not read as "a disc". */
9+function blob(cx: number, cy: number, rx: number, ry: number, phase: number): string {
10+ const pts: string[] = []
11+ for (let i = 0; i < 72; i++) {
12+ const t = (2 * Math.PI * i) / 72
13+ const r = 1 + 0.1 * Math.cos(3 * t + phase) + 0.045 * Math.sin(5 * t + 1.3)
14+ pts.push(`${(cx + rx * r * Math.cos(t)).toFixed(1)} ${(cy + ry * r * Math.sin(t)).toFixed(1)}`)
15+ }
16+ return `M ${pts.join(' L ')} Z`
17+}
18+
19+/** Outward normal of the blob at parameter t, for the Neumann panel. */
20+function blobPoint(cx: number, cy: number, rx: number, ry: number, phase: number, t: number) {
21+ const r = 1 + 0.1 * Math.cos(3 * t + phase) + 0.045 * Math.sin(5 * t + 1.3)
22+ const dr = -0.3 * Math.sin(3 * t + phase) + 0.225 * Math.cos(5 * t + 1.3)
23+ const x = cx + rx * r * Math.cos(t)
24+ const y = cy + ry * r * Math.sin(t)
25+ const dx = rx * (dr * Math.cos(t) - r * Math.sin(t))
26+ const dy = ry * (dr * Math.sin(t) + r * Math.cos(t))
27+ const n = Math.hypot(dx, dy)
28+ return { x, y, nx: dy / n, ny: -dx / n }
29+}
30+
31+const PANEL_W = 232
32+const CY = 108
33+
34+function Panel({
35+ x,
36+ heading,
37+ interior,
38+ sub,
39+ muted,
40+ children,
41+}: {
42+ x: number
43+ heading: string
44+ interior: string
45+ sub: string
46+ muted?: boolean
47+ children?: ReactNode
48+}) {
49+ return (
50+ <g transform={`translate(${x} 0)`}>
51+ <Label x={PANEL_W / 2} y={16} size={13} weight={600} fill={muted ? ink.muted : ink.primary}>
52+ {heading}
53+ </Label>
54+ {children}
55+ <Label x={PANEL_W / 2} y={CY} size={14} fill={muted ? ink.muted : ink.primary}>
56+ {interior}
57+ </Label>
58+ <Label x={PANEL_W / 2} y={214} size={12} fill={ink.muted}>
59+ {sub}
60+ </Label>
61+ </g>
62+ )
63+}
64+
65+/** Laplace against Poisson, Dirichlet against Neumann, in one row. */
66+export function DirichletProblem() {
67+ const cx = PANEL_W / 2
68+ const rx = 74
69+ const ry = 54
70+ const normals = [0.4, 1.5, 2.6, 3.7, 4.8, 5.9].map((t) => blobPoint(cx, CY, rx, ry, 0.7, t))
71+ return (
72+ <Figure
73+ width={716}
74+ height={232}
75+ caption="The interior condition and the boundary condition are chosen independently. This tutorial takes the first two panels: Laplace and Poisson with data prescribed on the whole boundary."
76+ >
77+ <Arrowheads id="dp" colors={[ink.muted]} />
78+
79+ <Panel x={0} heading="Laplace, Dirichlet" interior="Δu = 0" sub="Chapters 5–10">
80+ <path d={blob(cx, CY, rx, ry, 0.7)} fill="#20242c" stroke={method.fd2} strokeWidth={2.4} />
81+ <Label x={cx} y={CY + 22} size={11.5} fill={ink.muted}>
82+ no source
83+ </Label>
84+ <Label x={cx} y={CY + ry + 26} size={12} fill={method.fd2}>
85+ u = g on ∂Ω
86+ </Label>
87+ </Panel>
88+
89+ <Panel x={242} heading="Poisson, Dirichlet" interior="−Δu = f" sub="Chapters 1–4">
90+ <path d={blob(cx, CY, rx, ry, 0.7)} fill="#20242c" stroke={method.fd2} strokeWidth={2.4} />
91+ <ellipse cx={cx - 44} cy={CY - 24} rx={17} ry={13} fill={seqColor(0.75)} opacity={0.5} />
92+ <ellipse cx={cx + 46} cy={CY + 26} rx={12} ry={10} fill={seqColor(0.5)} opacity={0.45} />
93+ <Label x={cx} y={CY + 22} size={11.5} fill={ink.muted}>
94+ a source f inside
95+ </Label>
96+ <Label x={cx} y={CY + ry + 26} size={12} fill={method.fd2}>
97+ u = g on ∂Ω
98+ </Label>
99+ </Panel>
100+
101+ <Panel x={484} heading="Neumann" interior="Δu = 0" sub="not treated here" muted>
102+ <path d={blob(cx, CY, rx, ry, 0.7)} fill="#1e1e1c" stroke={ink.axis} strokeWidth={2.4} />
103+ {normals.map((p, i) => (
104+ <line
105+ key={i}
106+ x1={p.x}
107+ y1={p.y}
108+ x2={p.x + 15 * p.nx}
109+ y2={p.y + 15 * p.ny}
110+ stroke={ink.muted}
111+ strokeWidth={1.3}
112+ markerEnd={arrow('dp', 0)}
113+ />
114+ ))}
115+ <Label x={cx} y={CY + ry + 26} size={12} fill={ink.muted}>
116+ ∂u/∂n = g on ∂Ω
117+ </Label>
118+ </Panel>
119+ </Figure>
120+ )
121+}
122+
123+/** The four domains, in the order the tutorial meets them. */
124+export function Domains() {
125+ const c = ink.node
126+ const box = (i: number) => 12 + i * 176
127+ return (
128+ <Figure
129+ width={716}
130+ height={196}
131+ caption="The four domains. Each step to the right removes a structural assumption: the second dimension, then the convexity of the corner, then the straight edge itself."
132+ >
133+ {/* 1-D interval */}
134+ <g transform={`translate(${box(0)} 0)`}>
135+ <Label x={76} y={16} size={12.5} weight={600} fill={ink.primary}>
136+ interval
137+ </Label>
138+ <line x1={10} y1={92} x2={142} y2={92} stroke={ink.axis} strokeWidth={1.5} />
139+ {Array.from({ length: 9 }, (_, i) => (
140+ <Node key={i} x={10 + (i * 132) / 8} y={92} r={i === 0 || i === 8 ? 4.5 : 3} fill={i === 0 || i === 8 ? method.fd2 : c} />
141+ ))}
142+ <Label x={10} y={112} size={11} fill={method.fd2}>
143+ u(0)=α
144+ </Label>
145+ <Label x={142} y={112} size={11} fill={method.fd2}>
146+ u(1)=β
147+ </Label>
148+ <Label x={76} y={168} size={11.5} fill={ink.muted}>
149+ Chapters 1–4
150+ </Label>
151+ </g>
152+
153+ {/* unit square */}
154+ <g transform={`translate(${box(1)} 0)`}>
155+ <Label x={76} y={16} size={12.5} weight={600} fill={ink.primary}>
156+ square
157+ </Label>
158+ <rect x={22} y={38} width={108} height={108} fill="#20242c" stroke={method.fd2} strokeWidth={2.2} />
159+ {Array.from({ length: 5 }, (_, i) => (
160+ <g key={i}>
161+ <line x1={22 + (i + 1) * 18} y1={38} x2={22 + (i + 1) * 18} y2={146} stroke={ink.axis} strokeWidth={1} />
162+ <line x1={22} y1={38 + (i + 1) * 18} x2={130} y2={38 + (i + 1) * 18} stroke={ink.axis} strokeWidth={1} />
163+ </g>
164+ ))}
165+ <Label x={76} y={168} size={11.5} fill={ink.muted}>
166+ Chapters 5–8
167+ </Label>
168+ </g>
169+
170+ {/* L-shape */}
171+ <g transform={`translate(${box(2)} 0)`}>
172+ <Label x={76} y={16} size={12.5} weight={600} fill={ink.primary}>
173+ L-shape
174+ </Label>
175+ <path
176+ d="M 22 38 H 130 V 92 H 76 V 146 H 22 Z"
177+ fill="#20242c"
178+ stroke={method.iterative}
179+ strokeWidth={2.2}
180+ strokeLinejoin="round"
181+ />
182+ <circle cx={76} cy={92} r={5} fill="none" stroke={method.mps} strokeWidth={2} />
183+ <Label x={95} y={106} size={11} anchor="start" fill={method.mps}>
184+ 3π/2
185+ </Label>
186+ <Label x={76} y={168} size={11.5} fill={ink.muted}>
187+ Chapter 9
188+ </Label>
189+ </g>
190+
191+ {/* smooth curve */}
192+ <g transform={`translate(${box(3)} 0)`}>
193+ <Label x={76} y={16} size={12.5} weight={600} fill={ink.primary}>
194+ smooth curve
195+ </Label>
196+ <path
197+ d={blob(76, 92, 56, 54, 2.1)}
198+ fill="#20242c"
199+ stroke={method.bie}
200+ strokeWidth={2.2}
201+ />
202+ {Array.from({ length: 18 }, (_, i) => {
203+ const p = blobPoint(76, 92, 56, 54, 2.1, (2 * Math.PI * i) / 18)
204+ return <Node key={i} x={p.x} y={p.y} r={2.4} fill={ink.node} />
205+ })}
206+ <Label x={76} y={168} size={11.5} fill={ink.muted}>
207+ Chapter 10
208+ </Label>
209+ </g>
210+ </Figure>
211+ )
212+}
213+
214+/** Schematic of the spread the tutorial measures. Not data. */
215+export function AccuracyCost() {
216+ const L = 66
217+ const R = 470
218+ const T = 26
219+ const B = 176
220+ // log10 time from -1 to 3.5 ms across; log10 error from 0 at the top to -16
221+ const sx = (v: number) => L + ((v + 1) / 4.5) * (R - L)
222+ const sy = (v: number) => T + (-v / 16) * (B - T)
223+ const curve = (pts: [number, number][]) =>
224+ 'M ' + pts.map(([t, e]) => `${sx(t).toFixed(1)} ${sy(e).toFixed(1)}`).join(' L ')
225+
226+ const series: { key: keyof typeof method; label: string; pts: [number, number][] }[] = [
227+ { key: 'fd2', label: 'second order', pts: [[-0.9, -1.0], [0.0, -2.6], [1.2, -4.4], [2.4, -6.0], [3.4, -7.3]] },
228+ { key: 'fd4', label: 'sixth order', pts: [[-0.5, -1.4], [0.1, -5.0], [0.7, -9.0], [1.3, -13.0], [2.6, -13.4]] },
229+ { key: 'spectral', label: 'spectral', pts: [[-0.2, -2.0], [0.15, -6.5], [0.5, -11.5], [0.9, -14.0], [2.0, -14.0]] },
230+ { key: 'bie', label: 'boundary integral', pts: [[0.0, -3.0], [0.35, -9.0], [0.8, -15.4], [1.9, -15.6]] },
231+ ]
232+
233+ return (
234+ <Figure
235+ width={716}
236+ height={214}
237+ caption="Schematic, not measured data: the shape of the comparison every chapter returns to. Higher order changes the slope, and structure changes which curve is available at all. The labs produce the real numbers."
238+ >
239+ <line x1={L} y1={T} x2={L} y2={B} stroke={ink.axis} strokeWidth={1.2} />
240+ <line x1={L} y1={B} x2={R} y2={B} stroke={ink.axis} strokeWidth={1.2} />
241+ {[0, -4, -8, -12, -16].map((e) => (
242+ <g key={e}>
243+ <line x1={L} y1={sy(e)} x2={R} y2={sy(e)} stroke={ink.grid} strokeWidth={1} />
244+ <Label x={L - 8} y={sy(e)} size={10.5} anchor="end" fill={ink.muted} mono>
245+ {e === 0 ? '1' : `1e${e}`}
246+ </Label>
247+ </g>
248+ ))}
249+ <g transform={`translate(${L - 46} ${(T + B) / 2}) rotate(-90)`}>
250+ <Label x={0} y={0} size={11.5} fill={ink.muted}>
251+ error
252+ </Label>
253+ </g>
254+ <Label x={(L + R) / 2} y={B + 22} size={11.5} fill={ink.muted}>
255+ time for one solve →
256+ </Label>
257+
258+ {series.map((s, i) => {
259+ const last = s.pts[s.pts.length - 1]
260+ const ly = T + 16 + i * 22
261+ return (
262+ <g key={s.key}>
263+ <path d={curve(s.pts)} fill="none" stroke={method[s.key]} strokeWidth={2.2} strokeLinecap="round" />
264+ <Node x={sx(last[0])} y={sy(last[1])} r={3.4} fill={method[s.key]} />
265+ <line x1={R + 26} y1={ly} x2={R + 40} y2={ly} stroke={method[s.key]} strokeWidth={2.4} />
266+ <Label x={R + 48} y={ly} size={12} anchor="start" fill={method[s.key]}>
267+ {s.label}
268+ </Label>
269+ </g>
270+ )
271+ })}
272+ <Label x={R + 26} y={T + 16 + 4 * 22 + 6} size={11} anchor="start" fill={ink.muted}>
273+ each curve stops at
274+ </Label>
275+ <Label x={R + 26} y={T + 16 + 4 * 22 + 22} size={11} anchor="start" fill={ink.muted}>
276+ its own roundoff floor
277+ </Label>
278+ </Figure>
279+ )
280+}
src/figures/onedim.tsxadded+348−0View file
@@ -0,0 +1,348 @@
1+// Static figures for Chapters 1 to 4.
2+import { Arrowheads, Figure, Label, Node, arrow } from './svg.tsx'
3+import { ink, method } from '../plot/palette.ts'
4+
5+/** The mesh, the three-point stencil, and the band it puts in the matrix. */
6+export function ThreePointStencil() {
7+ const L = 70
8+ const R = 470
9+ const Y = 96
10+ const N = 8
11+ const xs = Array.from({ length: N + 1 }, (_, i) => L + (i * (R - L)) / N)
12+ const hit = [3, 4, 5]
13+ return (
14+ <Figure
15+ width={716}
16+ height={200}
17+ caption="The three-point scheme. Each interior node is coupled to its two neighbours and to nothing else, so the matrix has three non-zero diagonals and the solve is two sweeps."
18+ >
19+ <Arrowheads id="tps" colors={[ink.muted]} />
20+
21+ {/* the mesh */}
22+ <line x1={L} y1={Y} x2={R} y2={Y} stroke={ink.axis} strokeWidth={1.4} />
23+ {xs.map((x, i) => (
24+ <g key={i}>
25+ <Node x={x} y={Y} r={i === 0 || i === N ? 5 : 3.6} fill={i === 0 || i === N ? method.fd2 : ink.node} />
26+ <Label x={x} y={Y + 20} size={11} fill={ink.muted} mono>
27+ {i === N ? (
28+ <>
29+ u<tspan fontSize={9} dy={3}>N</tspan>
30+ </>
31+ ) : (
32+ `u${'₀₁₂₃₄₅₆₇₈'[i]}`
33+ )}
34+ </Label>
35+ </g>
36+ ))}
37+ <Label x={L} y={Y - 46} size={11.5} fill={method.fd2}>
38+ u(0) = α
39+ </Label>
40+ <Label x={R} y={Y - 46} size={11.5} fill={method.fd2}>
41+ u(1) = β
42+ </Label>
43+ <line x1={L} y1={Y - 36} x2={L} y2={Y - 10} stroke={method.fd2} strokeWidth={1} />
44+ <line x1={R} y1={Y - 36} x2={R} y2={Y - 10} stroke={method.fd2} strokeWidth={1} />
45+
46+ {/* h */}
47+ <line x1={xs[1]} y1={Y + 44} x2={xs[2]} y2={Y + 44} stroke={ink.muted} strokeWidth={1} markerStart={arrow('tps', 0)} markerEnd={arrow('tps', 0)} />
48+ <Label x={(xs[1] + xs[2]) / 2} y={Y + 60} size={11.5} fill={ink.muted} italic>
49+ h = 1/N
50+ </Label>
51+
52+ {/* the stencil, sitting on one interior node */}
53+ <path
54+ d={`M ${xs[hit[0]] - 12} ${Y - 26} H ${xs[hit[2]] + 12} V ${Y + 8} H ${xs[hit[0]] - 12} Z`}
55+ fill="none"
56+ stroke={method.fd2}
57+ strokeWidth={1.4}
58+ strokeDasharray="4 3"
59+ />
60+ {hit.map((i, k) => (
61+ <Label key={i} x={xs[i]} y={Y - 38} size={12} fill={method.fd2} mono>
62+ {[' −1', ' 2', ' −1'][k]}
63+ </Label>
64+ ))}
65+
66+ {/* the matrix band */}
67+ <g transform="translate(524 22)">
68+ <Label x={72} y={-4} size={11.5} fill={ink.muted}>
69+ the matrix, to scale
70+ </Label>
71+ {Array.from({ length: 7 }, (_, r) =>
72+ Array.from({ length: 7 }, (_, c) => {
73+ const d = Math.abs(r - c)
74+ if (d > 1) return null
75+ return (
76+ <rect
77+ key={`${r}-${c}`}
78+ x={16 + c * 19}
79+ y={8 + r * 19}
80+ width={17}
81+ height={17}
82+ fill={d === 0 ? method.fd2 : '#1c5cab'}
83+ opacity={d === 0 ? 0.95 : 0.6}
84+ rx={2}
85+ />
86+ )
87+ }),
88+ )}
89+ <rect x={14} y={6} width={137} height={137} fill="none" stroke={ink.axis} strokeWidth={1} rx={3} />
90+ <Label x={82} y={160} size={11} fill={ink.muted}>
91+ three diagonals, O(N) to solve
92+ </Label>
93+ </g>
94+ </Figure>
95+ )
96+}
97+
98+/** Where manufactured data comes from, and what is then compared with what. */
99+export function ManufacturedLoop() {
100+ const boxes = [
101+ { x: 6, w: 126, t: 'pick u(x)', s: 'any function you like' },
102+ { x: 148, w: 206, t: 'f := −u″, α := u(0), β := u(1)', s: 'differentiate by hand' },
103+ { x: 370, w: 146, t: 'solve for uₕ', s: 'the method under test' },
104+ { x: 532, w: 178, t: 'error = uₕ − u', s: 'known exactly, everywhere' },
105+ ]
106+ const Y = 42
107+ return (
108+ <Figure
109+ width={716}
110+ height={104}
111+ caption="The method of manufactured solutions. The data is derived from the answer rather than the other way round, so no reference computation stands between the method and its error."
112+ >
113+ <Arrowheads id="mfg" colors={[ink.muted]} />
114+ {boxes.map((b, i) => (
115+ <g key={i}>
116+ <rect x={b.x} y={Y} width={b.w} height={46} rx={6} fill="#201f1e" stroke={i === 3 ? method.fd2 : ink.axis} strokeWidth={1.4} />
117+ <Label x={b.x + b.w / 2} y={Y + 18} size={12.5} fill={ink.primary}>
118+ {b.t}
119+ </Label>
120+ <Label x={b.x + b.w / 2} y={Y + 35} size={10.5} fill={ink.muted}>
121+ {b.s}
122+ </Label>
123+ {i < boxes.length - 1 ? (
124+ <line
125+ x1={b.x + b.w + 4}
126+ y1={Y + 23}
127+ x2={boxes[i + 1].x - 4}
128+ y2={Y + 23}
129+ stroke={ink.muted}
130+ strokeWidth={1.3}
131+ markerEnd={arrow('mfg', 0)}
132+ />
133+ ) : null}
134+ </g>
135+ ))}
136+ <path
137+ d={`M 69 ${Y} C 69 14, 621 14, 621 ${Y}`}
138+ fill="none"
139+ stroke={ink.muted}
140+ strokeWidth={1.1}
141+ strokeDasharray="4 3"
142+ markerEnd={arrow('mfg', 0)}
143+ />
144+ <Label x={345} y={14} size={11} fill={ink.muted}>
145+ the same u, kept for comparison
146+ </Label>
147+ </Figure>
148+ )
149+}
150+
151+/** The V: discretisation down, roundoff up, and the best attainable error. */
152+export function RoundoffFloor() {
153+ const L = 66
154+ const R = 546
155+ const T = 24
156+ const B = 168
157+ // x is log10 N from 1 to 7; y is log10 error, 1 at the top and 1e-10 at the foot
158+ const sx = (v: number) => L + ((v - 1) / 6) * (R - L)
159+ const sy = (v: number) => T + (-v / 10) * (B - T)
160+ const disc = (n: number) => -2.3 - 2 * (n - 1)
161+ const round = (n: number) => -16 + 1.925 * n
162+ const total = (n: number) => Math.log10(10 ** disc(n) + 10 ** round(n))
163+ const path = (f: (n: number) => number) => {
164+ const pts: string[] = []
165+ for (let i = 0; i <= 240; i++) {
166+ const n = 1 + (6 * i) / 240
167+ const y = f(n)
168+ if (y > 0 || y < -10) continue
169+ pts.push(`${sx(n).toFixed(1)} ${sy(y).toFixed(1)}`)
170+ }
171+ return 'M ' + pts.join(' L ')
172+ }
173+ const nStar = 4 // where the two contributions are equal
174+ return (
175+ <Figure
176+ width={716}
177+ height={216}
178+ caption="Two error sources with opposite trends. Refining reduces the discretisation error and amplifies roundoff, so the total has a minimum: for a second-order scheme in double precision, near N = 10⁴ and an error of 10⁻⁸."
179+ >
180+ <line x1={L} y1={T} x2={L} y2={B} stroke={ink.axis} strokeWidth={1.2} />
181+ <line x1={L} y1={B} x2={R} y2={B} stroke={ink.axis} strokeWidth={1.2} />
182+ {[0, -2, -4, -6, -8, -10].map((e) => (
183+ <g key={e}>
184+ <line x1={L} y1={sy(e)} x2={R} y2={sy(e)} stroke={ink.grid} strokeWidth={1} />
185+ <Label x={L - 8} y={sy(e)} size={10.5} anchor="end" fill={ink.muted} mono>
186+ {e === 0 ? '1' : `1e${e}`}
187+ </Label>
188+ </g>
189+ ))}
190+ {[1, 2, 3, 4, 5, 6, 7].map((n) => (
191+ <Label key={n} x={sx(n)} y={B + 14} size={10.5} fill={ink.muted} mono baseline="hanging">
192+ {`1e${n}`}
193+ </Label>
194+ ))}
195+ <Label x={(L + R) / 2} y={B + 34} size={11.5} fill={ink.muted} baseline="hanging">
196+ intervals N
197+ </Label>
198+ <g transform={`translate(${L - 48} ${(T + B) / 2}) rotate(-90)`}>
199+ <Label x={0} y={0} size={11.5} fill={ink.muted}>
200+ error
201+ </Label>
202+ </g>
203+
204+ <path d={path(disc)} fill="none" stroke={ink.family} strokeWidth={1.8} strokeDasharray="5 4" />
205+ <path d={path(round)} fill="none" stroke={ink.family} strokeWidth={1.8} strokeDasharray="5 4" />
206+ <path d={path(total)} fill="none" stroke={method.fd2} strokeWidth={2.6} strokeLinecap="round" />
207+
208+ <Label x={sx(2.1)} y={sy(disc(2.1)) - 14} size={11.5} anchor="start" fill={ink.secondary}>
209+ discretisation, C h²
210+ </Label>
211+ <Label x={sx(6.4)} y={sy(round(6.4)) - 14} size={11.5} anchor="end" fill={ink.secondary}>
212+ roundoff, ε κ(A)
213+ </Label>
214+
215+ <line x1={sx(nStar)} y1={sy(total(nStar)) + 8} x2={sx(nStar)} y2={B} stroke={method.fd2} strokeWidth={1} strokeDasharray="3 3" />
216+ <Node x={sx(nStar)} y={sy(total(nStar))} r={4} fill={method.fd2} />
217+ <Label x={R + 12} y={sy(total(nStar)) - 8} size={12} anchor="start" fill={method.fd2}>
218+ best ≈ √ε
219+ </Label>
220+ <Label x={R + 12} y={sy(total(nStar)) + 9} size={11} anchor="start" fill={ink.muted}>
221+ about half the
222+ </Label>
223+ <Label x={R + 12} y={sy(total(nStar)) + 24} size={11} anchor="start" fill={ink.muted}>
224+ digits, and no more
225+ </Label>
226+ </Figure>
227+ )
228+}
229+
230+/** Numerov changes the right-hand side and leaves the matrix alone. */
231+export function CompactScheme() {
232+ const CX = [268, 336, 404]
233+ const cell = (x: number, y: number, v: string | null, colour: string) =>
234+ v === null ? (
235+ <rect x={x - 30} y={y} width={60} height={34} rx={5} fill="none" stroke={ink.grid} strokeWidth={1.2} strokeDasharray="3 3" />
236+ ) : (
237+ <>
238+ <rect x={x - 30} y={y} width={60} height={34} rx={5} fill="#201f1e" stroke={colour} strokeWidth={1.4} />
239+ <Label x={x} y={y + 17} size={12.5} fill={colour} mono>
240+ {v}
241+ </Label>
242+ </>
243+ )
244+ const row = (y: number, name: string, vals: (string | null)[], colour: string, note: string) => (
245+ <g>
246+ <Label x={12} y={y + 17} size={12.5} anchor="start" fill={ink.primary}>
247+ {name}
248+ </Label>
249+ {vals.map((v, i) => (
250+ <g key={i}>{cell(CX[i], y, v, colour)}</g>
251+ ))}
252+ <Label x={452} y={y + 17} size={11.5} anchor="start" fill={ink.muted}>
253+ {note}
254+ </Label>
255+ </g>
256+ )
257+ return (
258+ <Figure
259+ width={716}
260+ height={234}
261+ caption="Numerov's scheme against the three-point scheme. The left-hand side, and therefore the matrix and the cost of the solve, is identical; only the weighting of the data changes."
262+ >
263+ {['i−1', 'i', 'i+1'].map((t, i) => (
264+ <Label key={t} x={CX[i]} y={14} size={11.5} fill={ink.muted}>
265+ {t}
266+ </Label>
267+ ))}
268+
269+ <Label x={12} y={40} size={11} anchor="start" fill={ink.muted}>
270+ left-hand side, on the unknowns u
271+ </Label>
272+ {row(52, 'both schemes', ['−1', '2', '−1'], method.fd2, '÷ h², tridiagonal')}
273+
274+ <line x1={12} y1={110} x2={704} y2={110} stroke={ink.grid} strokeWidth={1} />
275+
276+ <Label x={12} y={130} size={11} anchor="start" fill={ink.muted}>
277+ right-hand side, on the data f
278+ </Label>
279+ {row(142, 'three-point', [null, 'fᵢ', null], ink.secondary, 'order 2')}
280+ {row(184, 'Numerov', ['f/12', '10f/12', 'f/12'], method.fd4, 'order 4, for one pass over an array')}
281+ </Figure>
282+ )
283+}
284+
285+/** Equispaced against Chebyshev points, and where they come from. */
286+export function ChebyshevPoints() {
287+ const N = 16
288+ const L = 172
289+ const R = 562
290+ const cx = (L + R) / 2
291+ const rad = (R - L) / 2
292+ const YARC = 216
293+ const cheb = Array.from({ length: N + 1 }, (_, j) => 0.5 * (1 - Math.cos((Math.PI * j) / N)))
294+ const equi = Array.from({ length: N + 1 }, (_, j) => j / N)
295+ const X = (t: number) => L + t * (R - L)
296+ return (
297+ <Figure
298+ width={716}
299+ height={300}
300+ caption="Chebyshev points are equally spaced angles on the semicircle, projected down. The clustering near the ends is what removes the Runge phenomenon, and it is also why the boundary-layer solution turns out to be the easy case for this method."
301+ >
302+ <path d={`M ${L} ${YARC} A ${rad} ${rad} 0 0 1 ${R} ${YARC}`} fill="none" stroke={ink.axis} strokeWidth={1.3} />
303+ {cheb.map((t, j) => {
304+ const th = Math.PI - (Math.PI * j) / N
305+ const ax = cx + rad * Math.cos(th)
306+ const ay = YARC - rad * Math.sin(th)
307+ return (
308+ <g key={j}>
309+ <line x1={ax} y1={ay} x2={X(t)} y2={YARC} stroke={ink.grid} strokeWidth={1} />
310+ <Node x={ax} y={ay} r={2.6} fill={method.spectral} />
311+ </g>
312+ )
313+ })}
314+ <Label x={20} y={64} size={11.5} anchor="start" fill={ink.muted}>
315+ equal angles
316+ </Label>
317+ <Label x={20} y={80} size={11.5} anchor="start" fill={ink.muted}>
318+ on the semicircle
319+ </Label>
320+
321+ <line x1={L} y1={YARC} x2={R} y2={YARC} stroke={ink.axis} strokeWidth={1.4} />
322+ {cheb.map((t, j) => (
323+ <Node key={j} x={X(t)} y={YARC} r={3.4} fill={method.spectral} />
324+ ))}
325+ <Label x={L - 14} y={YARC} size={12} anchor="end" fill={method.spectral}>
326+ Chebyshev
327+ </Label>
328+
329+ <line x1={L} y1={YARC + 44} x2={R} y2={YARC + 44} stroke={ink.axis} strokeWidth={1.4} />
330+ {equi.map((t, j) => (
331+ <Node key={j} x={X(t)} y={YARC + 44} r={3.4} fill={ink.family} />
332+ ))}
333+ <Label x={L - 14} y={YARC + 44} size={12} anchor="end" fill={ink.muted}>
334+ equispaced
335+ </Label>
336+ <Label x={R + 14} y={YARC + 44} size={11} anchor="start" fill={ink.muted}>
337+ unusable
338+ </Label>
339+ <Label x={R + 14} y={YARC + 59} size={11} anchor="start" fill={ink.muted}>
340+ past N ≈ 20
341+ </Label>
342+
343+ <Label x={cx} y={YARC + 78} size={11.5} fill={ink.muted}>
344+ N = 16 on [0, 1]
345+ </Label>
346+ </Figure>
347+ )
348+}
src/figures/svg.tsxadded+113−0View file
@@ -0,0 +1,113 @@
1+// Static figures: hand-drawn SVG, no data and no session behind them.
2+//
3+// A lab shows what a method does once the reader presses run; these show what
4+// it *is*, on first sight and in print. They are deliberately plain: the same
5+// palette as the plots, unicode rather than KaTeX for the few symbols that
6+// appear, and one idea per picture.
7+import type { ReactNode } from 'react'
8+import { ink } from '../plot/palette.ts'
9+
10+export const SANS = 'ui-sans-serif, system-ui, -apple-system, Segoe UI, Roboto, sans-serif'
11+export const MONO = 'ui-monospace, SFMono-Regular, Menlo, Consolas, monospace'
12+
13+/** Wrapper: a captioned, width-filling SVG with a fixed aspect ratio. */
14+export function Figure({
15+ width,
16+ height,
17+ caption,
18+ maxWidth,
19+ children,
20+}: {
21+ width: number
22+ height: number
23+ caption?: string
24+ /** in px; the figure never grows past this even in a wide column */
25+ maxWidth?: number
26+ children: ReactNode
27+}) {
28+ return (
29+ <figure className="static-figure" style={maxWidth ? { maxWidth } : undefined}>
30+ <svg
31+ viewBox={`0 0 ${width} ${height}`}
32+ width="100%"
33+ preserveAspectRatio="xMidYMid meet"
34+ role="img"
35+ aria-label={caption ?? 'diagram'}
36+ style={{ display: 'block', overflow: 'visible' }}
37+ >
38+ {children}
39+ </svg>
40+ {caption ? <figcaption>{caption}</figcaption> : null}
41+ </figure>
42+ )
43+}
44+
45+/** Arrowhead markers, one set per figure so the ids never collide. */
46+export function Arrowheads({ id, colors }: { id: string; colors: string[] }) {
47+ return (
48+ <defs>
49+ {colors.map((c, i) => (
50+ <marker
51+ key={i}
52+ id={`${id}-ah-${i}`}
53+ viewBox="0 0 8 8"
54+ refX="7"
55+ refY="4"
56+ markerWidth="6"
57+ markerHeight="6"
58+ orient="auto-start-reverse"
59+ >
60+ <path d="M 0 0.7 L 8 4 L 0 7.3 z" fill={c} />
61+ </marker>
62+ ))}
63+ </defs>
64+ )
65+}
66+
67+export const arrow = (id: string, i: number) => `url(#${id}-ah-${i})`
68+
69+/** A text label. `anchor` and `baseline` are the usual SVG values. */
70+export function Label({
71+ x,
72+ y,
73+ children,
74+ fill = ink.secondary,
75+ size = 12,
76+ anchor = 'middle',
77+ baseline = 'middle',
78+ italic = false,
79+ mono = false,
80+ weight,
81+}: {
82+ x: number
83+ y: number
84+ children: ReactNode
85+ fill?: string
86+ size?: number
87+ anchor?: 'start' | 'middle' | 'end'
88+ baseline?: 'middle' | 'hanging' | 'auto' | 'central'
89+ italic?: boolean
90+ mono?: boolean
91+ weight?: number
92+}) {
93+ return (
94+ <text
95+ x={x}
96+ y={y}
97+ fill={fill}
98+ fontSize={size}
99+ fontFamily={mono ? MONO : SANS}
100+ fontStyle={italic ? 'italic' : undefined}
101+ fontWeight={weight}
102+ textAnchor={anchor}
103+ dominantBaseline={baseline}
104+ >
105+ {children}
106+ </text>
107+ )
108+}
109+
110+/** A grid node. */
111+export function Node({ x, y, r = 3.2, fill = ink.node, stroke }: { x: number; y: number; r?: number; fill?: string; stroke?: string }) {
112+ return <circle cx={x} cy={y} r={r} fill={fill} stroke={stroke} strokeWidth={stroke ? 1.2 : undefined} />
113+}
src/figures/twodim.tsxadded+363−0View file
@@ -0,0 +1,363 @@
1+// Static figures for Chapters 5 to 8.
2+import type { ReactNode } from 'react'
3+import { Arrowheads, Figure, Label, Node, arrow } from './svg.tsx'
4+import { ink, method, seqColor } from '../plot/palette.ts'
5+
6+/** The grid, the five-point stencil, and the sparsity it produces. */
7+export function FivePointStencil() {
8+ const n = 8
9+ const s = 20
10+ const X0 = 42
11+ const Y0 = 34
12+ const gx = (i: number) => X0 + i * s
13+ const gy = (j: number) => Y0 + (n - j) * s
14+ const ci = 3
15+ const cj = 5
16+ const nb: [number, number][] = [
17+ [ci - 1, cj],
18+ [ci + 1, cj],
19+ [ci, cj - 1],
20+ [ci, cj + 1],
21+ ]
22+ return (
23+ <Figure
24+ width={716}
25+ height={230}
26+ caption="The five-point Laplacian on a uniform grid, N = 8. Boundary nodes (blue) carry known data; each of the (N−1)² interior nodes gives one equation coupling it to its four neighbours."
27+ >
28+ {/* grid lines */}
29+ {Array.from({ length: n + 1 }, (_, i) => (
30+ <g key={i}>
31+ <line x1={gx(i)} y1={gy(0)} x2={gx(i)} y2={gy(n)} stroke={ink.grid} strokeWidth={1} />
32+ <line x1={gx(0)} y1={gy(i)} x2={gx(n)} y2={gy(i)} stroke={ink.grid} strokeWidth={1} />
33+ </g>
34+ ))}
35+ <rect x={gx(0)} y={gy(n)} width={n * s} height={n * s} fill="none" stroke={method.fd2} strokeWidth={2} />
36+
37+ {/* the stencil arms, drawn under the nodes */}
38+ {nb.map(([i, j], k) => (
39+ <line key={k} x1={gx(ci)} y1={gy(cj)} x2={gx(i)} y2={gy(j)} stroke={method.fd2} strokeWidth={2.4} />
40+ ))}
41+
42+ {/* nodes */}
43+ {Array.from({ length: n + 1 }, (_, i) =>
44+ Array.from({ length: n + 1 }, (_, j) => {
45+ const bnd = i === 0 || j === 0 || i === n || j === n
46+ const inStencil = (i === ci && j === cj) || nb.some(([a, b]) => a === i && b === j)
47+ return (
48+ <Node
49+ key={`${i}-${j}`}
50+ x={gx(i)}
51+ y={gy(j)}
52+ r={inStencil ? 4.6 : bnd ? 3.4 : 2.4}
53+ fill={inStencil ? method.fd2 : bnd ? '#256abf' : '#5a5a53'}
54+ />
55+ )
56+ }),
57+ )}
58+ <Label x={gx(ci) - 15} y={gy(cj) - 13} size={11.5} fill={ink.primary} mono>
59+ −4
60+ </Label>
61+ {nb.map(([i, j], k) => (
62+ <Label
63+ key={k}
64+ x={gx(i) + (j === cj ? (i < ci ? -13 : 13) : 0)}
65+ y={gy(j) + (j === cj ? 0 : j > cj ? -13 : 13)}
66+ size={11}
67+ fill={ink.secondary}
68+ mono
69+ >
70+ 1
71+ </Label>
72+ ))}
73+ <Label x={gx(n) + 8} y={gy(Math.round(n / 2))} size={11.5} anchor="start" fill={method.fd2}>
74+ u = g
75+ </Label>
76+
77+ {/* the sparsity pattern of A = I⊗T + T⊗I */}
78+ <g transform="translate(300 20)">
79+ <Label x={92} y={-4} size={11.5} fill={ink.muted}>
80+ A = I ⊗ T + T ⊗ I, for N = 6
81+ </Label>
82+ {(() => {
83+ const m = 5
84+ const K = m * m
85+ const px = 168 / K
86+ const cells: ReactNode[] = []
87+ for (let r = 0; r < K; r++) {
88+ for (let c = 0; c < K; c++) {
89+ const d = Math.abs(r - c)
90+ if (d !== 0 && d !== 1 && d !== m) continue
91+ if (d === 1 && Math.floor(r / m) !== Math.floor(c / m)) continue
92+ cells.push(
93+ <rect
94+ key={`${r}-${c}`}
95+ x={8 + c * px}
96+ y={8 + r * px}
97+ width={px}
98+ height={px}
99+ fill={d === 0 ? method.fd2 : '#1c5cab'}
100+ opacity={d === 0 ? 0.95 : 0.65}
101+ />,
102+ )
103+ }
104+ }
105+ return cells
106+ })()}
107+ <rect x={8} y={8} width={168} height={168} fill="none" stroke={ink.axis} strokeWidth={1} />
108+ <Label x={92} y={192} size={11} fill={ink.muted}>
109+ five non-zeros per row, and dense elimination ignores that
110+ </Label>
111+ </g>
112+ </Figure>
113+ )
114+}
115+
116+/** The sine-transform solve, as three steps. */
117+export function FastPoissonSteps() {
118+ const boxes = [
119+ { x: 10, w: 160, t: 'F', s: 'right-hand side on the grid', c: ink.axis },
120+ { x: 216, w: 176, t: 'F̃ = S F S', s: '2-D sine transform, O(N² log N)', c: method.fast },
121+ { x: 438, w: 176, t: 'Ũ = F̃ / (λₖ + λₗ)', s: 'a division, one per unknown', c: method.fast },
122+ ]
123+ return (
124+ <Figure
125+ width={716}
126+ height={182}
127+ caption="The fast Poisson solver. The sine vectors diagonalise the five-point matrix exactly, so the solve is a transform, a division and a transform back. The answer is the same one dense elimination gives, to the last bit."
128+ >
129+ <Arrowheads id="fps" colors={[ink.muted]} />
130+ {boxes.map((b, i) => (
131+ <g key={i}>
132+ <rect x={b.x} y={30} width={b.w} height={50} rx={6} fill="#201f1e" stroke={b.c} strokeWidth={1.4} />
133+ <Label x={b.x + b.w / 2} y={49} size={13} fill={ink.primary}>
134+ {b.t}
135+ </Label>
136+ <Label x={b.x + b.w / 2} y={67} size={10.5} fill={ink.muted}>
137+ {b.s}
138+ </Label>
139+ {i < boxes.length - 1 ? (
140+ <line x1={b.x + b.w + 8} y1={55} x2={boxes[i + 1].x - 8} y2={55} stroke={ink.muted} strokeWidth={1.3} markerEnd={arrow('fps', 0)} />
141+ ) : null}
142+ </g>
143+ ))}
144+ <rect x={438} y={106} width={176} height={50} rx={6} fill="#201f1e" stroke={method.fast} strokeWidth={1.4} />
145+ <Label x={526} y={125} size={13} fill={ink.primary}>
146+ U = S Ũ S
147+ </Label>
148+ <Label x={526} y={143} size={10.5} fill={ink.muted}>
149+ transform back
150+ </Label>
151+ <line x1={526} y1={84} x2={526} y2={102} stroke={ink.muted} strokeWidth={1.3} markerEnd={arrow('fps', 0)} />
152+
153+ <Label x={216} y={131} size={12} anchor="start" fill={ink.secondary}>
154+ λₖ = 4 sin²(πk / 2N)
155+ </Label>
156+ <Label x={216} y={149} size={11} anchor="start" fill={ink.muted}>
157+ known in closed form
158+ </Label>
159+ </Figure>
160+ )
161+}
162+
163+/** Five-point against the compact nine-point. */
164+export function NinePointStencil() {
165+ const s = 44
166+ const draw = (ox: number, oy: number, w: (di: number, dj: number) => string | null, colour: string) => {
167+ const cells: { x: number; y: number; v: string }[] = []
168+ for (const dj of [-1, 0, 1]) {
169+ for (const di of [-1, 0, 1]) {
170+ const v = w(di, dj)
171+ if (v !== null) cells.push({ x: ox + di * s, y: oy - dj * s, v })
172+ }
173+ }
174+ return (
175+ <g>
176+ {cells.map((c, i) =>
177+ c.x === ox && c.y === oy ? null : (
178+ <line key={`l${i}`} x1={ox} y1={oy} x2={c.x} y2={c.y} stroke={colour} strokeWidth={1.6} opacity={0.55} />
179+ ),
180+ )}
181+ {cells.map((c, i) => (
182+ <g key={`c${i}`}>
183+ <circle cx={c.x} cy={c.y} r={15} fill="#201f1e" stroke={colour} strokeWidth={1.6} />
184+ <Label x={c.x} y={c.y} size={12} fill={colour} mono>
185+ {c.v}
186+ </Label>
187+ </g>
188+ ))}
189+ </g>
190+ )
191+ }
192+ return (
193+ <Figure
194+ width={716}
195+ height={220}
196+ caption="The compact nine-point (Mehrstellen) stencil beside the five-point one. Both couple a node only to its immediate neighbours, so both are diagonalised by the same sine vectors and cost the same to solve; on Laplace's equation the nine-point scheme is sixth-order accurate."
197+ >
198+ <Label x={168} y={18} size={12.5} weight={600} fill={ink.primary}>
199+ five-point, order 2
200+ </Label>
201+ {draw(168, 116, (di, dj) => (di === 0 && dj === 0 ? '−4' : di === 0 || dj === 0 ? '1' : null), method.fd2)}
202+ <Label x={168} y={202} size={11.5} fill={ink.muted}>
203+ ÷ h²
204+ </Label>
205+
206+ <line x1={358} y1={40} x2={358} y2={190} stroke={ink.grid} strokeWidth={1} />
207+
208+ <Label x={548} y={18} size={12.5} weight={600} fill={ink.primary}>
209+ nine-point, order 6 on Δu = 0
210+ </Label>
211+ {draw(548, 116, (di, dj) => (di === 0 && dj === 0 ? '−20' : di === 0 || dj === 0 ? '4' : '1'), method.fd4)}
212+ <Label x={548} y={202} size={11.5} fill={ink.muted}>
213+ ÷ 6h²
214+ </Label>
215+ </Figure>
216+ )
217+}
218+
219+/** The multigrid V-cycle. */
220+export function VCycle() {
221+ const levels = [
222+ { label: 'h', w: 300 },
223+ { label: '2h', w: 210 },
224+ { label: '4h', w: 140 },
225+ { label: '8h', w: 84 },
226+ ]
227+ const y = (k: number) => 34 + k * 42
228+ const XC = 358
229+ const xl = (k: number) => XC - levels[k].w / 2
230+ const xr = (k: number) => XC + levels[k].w / 2
231+ return (
232+ <Figure
233+ width={716}
234+ height={218}
235+ caption="One multigrid V-cycle. A few Gauss–Seidel sweeps make the error smooth; a smooth error is represented on a grid twice as coarse at a quarter of the cost, and the argument repeats. The whole cycle costs O(N²) and reduces the error by a factor that does not depend on N."
236+ >
237+ <Arrowheads id="vc" colors={[method.iterative, ink.muted]} />
238+ {levels.map((L, k) => (
239+ <g key={k}>
240+ <line x1={xl(k)} y1={y(k)} x2={xr(k)} y2={y(k)} stroke={ink.grid} strokeWidth={1} />
241+ <Label x={xl(k) - 16} y={y(k)} size={11.5} anchor="end" fill={ink.muted} mono>
242+ {L.label}
243+ </Label>
244+ {(() => {
245+ const m = Math.max(3, 9 - 2 * k)
246+ return Array.from({ length: m }, (_, i) => (
247+ <Node key={i} x={xl(k) + (i * L.w) / (m - 1)} y={y(k)} r={2.6} fill={ink.family} />
248+ ))
249+ })()}
250+ </g>
251+ ))}
252+
253+ {/* descending and ascending legs */}
254+ {[0, 1, 2].map((k) => (
255+ <line
256+ key={`d${k}`}
257+ x1={xl(k) + 26}
258+ y1={y(k) + 6}
259+ x2={xl(k + 1) + 20}
260+ y2={y(k + 1) - 6}
261+ stroke={method.iterative}
262+ strokeWidth={2}
263+ markerEnd={arrow('vc', 0)}
264+ />
265+ ))}
266+ {[2, 1, 0].map((k) => (
267+ <line
268+ key={`u${k}`}
269+ x1={xr(k + 1) - 20}
270+ y1={y(k + 1) - 6}
271+ x2={xr(k) - 26}
272+ y2={y(k) + 6}
273+ stroke={method.iterative}
274+ strokeWidth={2}
275+ markerEnd={arrow('vc', 0)}
276+ />
277+ ))}
278+
279+ <Label x={152} y={98} size={11.5} anchor="end" fill={method.iterative}>
280+ smooth,
281+ </Label>
282+ <Label x={152} y={114} size={11.5} anchor="end" fill={method.iterative}>
283+ restrict residual
284+ </Label>
285+ <Label x={564} y={98} size={11.5} anchor="start" fill={method.iterative}>
286+ interpolate correction,
287+ </Label>
288+ <Label x={564} y={114} size={11.5} anchor="start" fill={method.iterative}>
289+ smooth again
290+ </Label>
291+ <Label x={XC} y={y(3) + 22} size={11.5} fill={ink.secondary}>
292+ solve exactly: the coarse grid is tiny
293+ </Label>
294+
295+ <Label x={XC} y={198} size={11.5} fill={ink.muted}>
296+ cost of the whole cycle ≈ 4/3 of one fine-grid sweep
297+ </Label>
298+ </Figure>
299+ )
300+}
301+
302+/** What smoothing does to the error, as a picture rather than a claim. */
303+export function SmoothingIdea() {
304+ const L = 40
305+ const R = 320
306+ const mid = 62
307+ const wig = (amp: number, hf: number, phase: number) => {
308+ const pts: string[] = []
309+ for (let i = 0; i <= 140; i++) {
310+ const t = i / 140
311+ const v =
312+ amp * Math.sin(Math.PI * t) * 0.75 +
313+ hf * (Math.sin(19 * Math.PI * t + phase) * 0.5 + Math.sin(31 * Math.PI * t) * 0.35)
314+ pts.push(`${(L + t * (R - L)).toFixed(1)} ${(mid - v).toFixed(1)}`)
315+ }
316+ return 'M ' + pts.join(' L ')
317+ }
318+ return (
319+ <Figure
320+ width={716}
321+ height={140}
322+ caption="Two Gauss–Seidel sweeps applied to a rough error. What is left is smooth, and a smooth function loses nothing by being carried to a grid twice as coarse. Gauss–Seidel is a poor solver and a good smoother, and multigrid uses only the second property."
323+ >
324+ <Arrowheads id="sm" colors={[ink.muted]} />
325+ <Label x={(L + R) / 2} y={16} size={12} fill={ink.muted}>
326+ error before
327+ </Label>
328+ <line x1={L} y1={mid} x2={R} y2={mid} stroke={ink.grid} strokeWidth={1} />
329+ <path d={wig(26, 11, 0.4)} fill="none" stroke={method.iterative} strokeWidth={1.7} />
330+
331+ <line x1={R + 26} y1={mid} x2={R + 66} y2={mid} stroke={ink.muted} strokeWidth={1.4} markerEnd={arrow('sm', 0)} />
332+ <Label x={R + 46} y={mid - 12} size={11} fill={ink.muted}>
333+ 2 sweeps
334+ </Label>
335+
336+ <g transform={`translate(${R + 92} 0)`}>
337+ <Label x={(L + R) / 2 - L + 8} y={16} size={12} fill={ink.muted}>
338+ error after
339+ </Label>
340+ <line x1={8} y1={mid} x2={8 + (R - L)} y2={mid} stroke={ink.grid} strokeWidth={1} />
341+ <path
342+ d={(() => {
343+ const pts: string[] = []
344+ for (let i = 0; i <= 140; i++) {
345+ const t = i / 140
346+ pts.push(`${(8 + t * (R - L)).toFixed(1)} ${(mid - 26 * Math.sin(Math.PI * t) * 0.75).toFixed(1)}`)
347+ }
348+ return 'M ' + pts.join(' L ')
349+ })()}
350+ fill="none"
351+ stroke={method.iterative}
352+ strokeWidth={1.7}
353+ />
354+ {Array.from({ length: 8 }, (_, i) => (
355+ <Node key={i} x={8 + (i * (R - L)) / 7} y={mid} r={3} fill={seqColor(0.8)} />
356+ ))}
357+ <Label x={8 + (R - L) / 2} y={mid + 40} size={11} fill={ink.muted}>
358+ eight coarse points now suffice
359+ </Label>
360+ </g>
361+ </Figure>
362+ )
363+}
src/index.cssmodified+14−0View file
@@ -419,6 +419,20 @@ button.ghost:hover {
419419 font-size: 12.5px;
420420 }
421421
422+/* ---------- static figures ---------- */
423+
424+.static-figure {
425+ margin: 22px 0;
426+ max-width: 720px;
427+}
428+
429+.static-figure figcaption {
430+ margin-top: 8px;
431+ font-size: 12.5px;
432+ line-height: 1.5;
433+ color: var(--muted);
434+}
435+
422436 /* ---------- plots ---------- */
423437
424438 .plot-host {
src/matlab/lp_curve.mmodified+2−2View file
@@ -4,8 +4,8 @@ function C = lp_curve(a, nlobe, t)
44 % r(t) = 1 + a cos(nlobe t), x = r cos t, y = r sin t
55 %
66 % 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.
7+% outputs: the callers want all six together and a single value keeps the
8+% call sites short.
99 %
1010 % The curve is analytic and traversed counterclockwise, which is what makes the
1111 % trapezoid rule spectrally accurate on it and fixes the sign of the outward
src/matlab/lp_orders.mmodified+2−2View file
@@ -2,8 +2,8 @@ function q = lp_orders(e, n)
22 % LP_ORDERS Observed convergence order between successive refinements.
33 %
44 % 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).
5+% a positive finite number — an error of exactly zero, or a resolution a method
6+% skipped, would otherwise put an Inf or a complex value in the table.
77 q = NaN(size(e));
88 for k = 2:numel(e)
99 r = e(k-1) / e(k);
src/ui/ChapterView.tsxmodified+14−6View file
@@ -1,27 +1,30 @@
1-// A chapter: prose with its labs dropped in where the text calls for them.
1+// A chapter: prose with its labs and figures dropped in where the text calls
2+// for them.
23 //
34 // The prose is one markdown file per chapter, with a line
45 //
56 // <!--lab:convergence-->
7+// <!--fig:five-point-->
68 //
7-// wherever a lab belongs. Splitting on those markers keeps the writing in one
8-// readable file instead of scattered through JSX.
9+// wherever a lab or a static figure belongs. Splitting on those markers keeps
10+// the writing in one readable file instead of scattered through JSX.
911 import { useMemo } from 'react'
1012 import Markdown from './Markdown.tsx'
1113 import Lab from './Lab.tsx'
1214 import type { Chapter } from './types.ts'
15+import { FIGURES } from '../figures/index.ts'
1316
14-const MARKER = /<!--\s*lab:([a-z0-9_-]+)\s*-->/gi
17+const MARKER = /<!--\s*(lab|fig):([a-z0-9_-]+)\s*-->/gi
1518
1619 export default function ChapterView({ chapter }: { chapter: Chapter }) {
1720 const pieces = useMemo(() => {
18- const out: ({ md: string } | { labId: string })[] = []
21+ const out: ({ md: string } | { labId: string } | { figId: string })[] = []
1922 let last = 0
2023 MARKER.lastIndex = 0
2124 let m: RegExpExecArray | null
2225 while ((m = MARKER.exec(chapter.text)) !== null) {
2326 out.push({ md: chapter.text.slice(last, m.index) })
24- out.push({ labId: m[1] })
27+ out.push(m[1] === 'fig' ? { figId: m[2] } : { labId: m[2] })
2528 last = m.index + m[0].length
2629 }
2730 out.push({ md: chapter.text.slice(last) })
@@ -37,6 +40,11 @@ export default function ChapterView({ chapter }: { chapter: Chapter }) {
3740 </header>
3841 {pieces.map((p, i) => {
3942 if ('md' in p) return p.md.trim() === '' ? null : <Markdown key={i} src={p.md} />
43+ if ('figId' in p) {
44+ const Fig = FIGURES[p.figId]
45+ if (!Fig) return <p key={i} className="lab-error">missing figure: {p.figId}</p>
46+ return <Fig key={`${chapter.id}-fig-${p.figId}`} />
47+ }
4048 const lab = chapter.labs.find((l) => l.id === p.labId)
4149 if (!lab) return <p key={i} className="lab-error">missing lab: {p.labId}</p>
4250 return <Lab key={`${chapter.id}-${lab.id}`} spec={lab} />
src/ui/types.tsmodified+1−1View file
@@ -52,7 +52,7 @@ export interface Chapter {
5252 title: string
5353 /** one sentence, shown under the title and in the contents */
5454 blurb: string
55- part: 'One dimension' | 'Two dimensions' | 'Beyond the square'
55+ part: 'Introduction' | 'One dimension' | 'Two dimensions' | 'Beyond the square'
5656 /** the prose, with <!--lab:id--> lines marking where each lab goes */
5757 text: string
5858 labs: LabSpec[]