e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 1// Static figures for Chapters 1 to 4.
2import { Arrowheads, Figure, Label, Node, arrow } from './svg.tsx'
3import { ink, method } from '../plot/palette.ts'
5/** The mesh, the three-point stencil, and the band it puts in the matrix. */
6export 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]} />
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} />
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>
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 ))}
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}
98/** Where manufactured data comes from, and what is then compared with what. */
99export 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}
151/** The V: discretisation down, roundoff up, and the best attainable error. */
152export 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>
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" />
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>
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}
230/** Numerov changes the right-hand side and leaves the matrix alone. */
231export 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 ))}
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')}
274 <line x1={12} y1={110} x2={704} y2={110} stroke={ink.grid} strokeWidth={1} />
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}
285/** Equispaced against Chebyshev points, and where they come from. */
286export 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>
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>
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>
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}