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