e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 1// Static figures for Chapters 5 to 8.
2import type { ReactNode } from 'react'
3import { Arrowheads, Figure, Label, Node, arrow } from './svg.tsx'
4import { ink, method, seqColor } from '../plot/palette.ts'
6/** The grid, the five-point stencil, and the sparsity it produces. */
7export 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} />
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 ))}
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>
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}
116/** The sine-transform solve, as three steps. */
117export 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)} />
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}
163/** Five-point against the compact nine-point. */
164export 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>
206 <line x1={358} y1={40} x2={358} y2={190} stroke={ink.grid} strokeWidth={1} />
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}
219/** The multigrid V-cycle. */
220export 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 ))}
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 ))}
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>
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}
302/** What smoothing does to the error, as a picture rather than a claim. */
303export 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} />
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>
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}