088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 1// Figures for the iterative-solver chapter: how the residual falls per
2// iteration, how the true error falls, and how both compare with a direct
3// solve of the same system.
4import Plot from '../../plot/Plot.tsx'
5import Legend from '../../plot/Legend.tsx'
6import { EndLabel, Line, Markers } from '../../plot/Curves.tsx'
7import { ink, marker, method, type MethodKey } from '../../plot/palette.ts'
8import type { Num } from '../../engine/types.ts'
9import { fmtErr } from '../onedim/panels.tsx'
11export interface IterOut {
12 names: string[]
13 N: number
14 its: number[]
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 15 R: Num[][]
16 E: Num[][]
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 17 T: Num[][]
18 discErr: number
19 directMs: number
20 unknowns: number
21}
23// one palette slot per method, in the order the script lists them
24const SLOTS: MethodKey[] = ['fd2', 'fd4', 'fast', 'spectral', 'iterative']
26function lowHigh(rows: Num[][], floor = 1e-14): [number, number] {
27 let lo = Infinity
28 let hi = -Infinity
29 for (const r of rows) {
30 for (const v of r) {
31 if (v == null || !isFinite(v) || v <= 0) continue
32 if (v < lo) lo = v
33 if (v > hi) hi = v
34 }
35 }
36 if (!isFinite(lo)) return [floor, 1]
37 return [Math.max(lo / 3, floor), hi * 3]
38}
40function lastIndex(row: Num[]): number {
41 for (let i = row.length - 1; i >= 0; i--) if (row[i] != null && isFinite(row[i] as number)) return i
42 return 0
43}
45export function IterPanels({ data }: { data: IterOut }) {
46 const its = data.its
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 47 const R = data.R
48 const E = data.E
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 49 const xd: [number, number] = [0, its[its.length - 1]]
50 const rd = lowHigh(R)
51 const ed = lowHigh([...E, [data.discErr]])
52 const td = lowHigh(data.T.map((r) => r.map((v) => (v != null && v > 0 ? v : null))))
53 td[0] = Math.max(td[0], 0.05)
55 const legend = data.names.map((n, k) => ({ label: n, color: method[SLOTS[k]] }))
57 return (
58 <>
59 <div className="plot-title">relative residual against iteration count</div>
60 <Plot height={300} xDomain={xd} yDomain={rd} yLog xLabel="iteration" yLabel="‖r‖ / ‖r₀‖">
61 {(f) => (
62 <>
63 {R.map((row, k) => (
64 <g key={k}>
65 <Line x={its} y={row} f={f} color={method[SLOTS[k]]} />
66 {row.length <= 40 && (
67 <Markers x={its} y={row} f={f} color={method[SLOTS[k]]} shape={marker[SLOTS[k]]} r={3} />
68 )}
69 </g>
70 ))}
71 </>
72 )}
73 </Plot>
74 <Legend items={legend} />
76 <div className="plot-title">
77 true error against iteration count; the dashed line is the discretisation error, which is where the exact
78 solution of this linear system sits
79 </div>
80 <Plot
81 height={300}
82 xDomain={xd}
83 yDomain={ed}
84 yLog
85 xLabel="iteration"
86 yLabel="max |error|"
87 margin={{ r: 30 }}
88 over={(f) => (
89 <>
90 {E.map((row, k) => {
91 const i = lastIndex(row)
92 return (
93 <EndLabel
94 key={k}
95 f={f}
96 x={its[i]}
97 y={row[i] as number}
98 text={data.names[k]}
99 color={method[SLOTS[k]]}
100 anchor="end"
101 dx={-6}
102 dy={-7}
103 />
104 )
105 })}
106 </>
107 )}
108 >
109 {(f) => (
110 <>
111 <Line x={xd} y={[data.discErr, data.discErr]} f={f} color={ink.muted} dash="6 4" width={1.5} />
112 {E.map((row, k) => (
113 <Line key={k} x={its} y={row} f={f} color={method[SLOTS[k]]} />
114 ))}
115 </>
116 )}
117 </Plot>
118 <Legend
119 items={[...legend, { label: 'discretisation error', color: ink.muted, dash: '6 4' }]}
120 />
122 <div className="plot-title">
123 true error against milliseconds of solver time; the marked point is the sine-transform direct solve of the same
124 system
125 </div>
126 <Plot height={300} xDomain={td} yDomain={ed} xLog yLog xLabel="milliseconds" yLabel="max |error|">
127 {(f) => (
128 <>
129 {E.map((row, k) => (
130 <Line
131 key={k}
132 x={(data.T[k] ?? []).map((v) => (v == null || v <= 0 ? NaN : v))}
133 y={row}
134 f={f}
135 color={method[SLOTS[k]]}
136 />
137 ))}
138 <Markers
139 x={[data.directMs]}
140 y={[data.discErr]}
141 f={f}
142 color={ink.node}
143 shape="star"
144 r={7}
145 />
146 <EndLabel
147 f={f}
148 x={data.directMs}
149 y={data.discErr}
150 text="direct solve"
151 color={ink.secondary}
152 dx={-12}
153 dy={-10}
154 anchor="end"
155 />
156 </>
157 )}
158 </Plot>
159 <Legend items={[...legend, { label: 'direct solve', color: ink.node, shape: 'cross' }]} />
160 </>
161 )
162}
164export function iterSummary(d: IterOut) {
e5fc8c6An introduction and static figures, and the MATLAB written plainly againJeremy Magland 165 const mg = d.E[d.E.length - 1]
088be4cAn interactive tutorial on the Dirichlet problem, accuracy against costJeremy Magland 166 const i = lastIndex(mg)
167 return [
168 { label: 'unknowns', value: String(d.unknowns) },
169 { label: 'discretisation error', value: fmtErr(d.discErr) },
170 { label: 'multigrid cycles used', value: String(d.its[i]) },
171 ]
172}