2// Checks the MATLAB layer against the numbers the chapters claim, through the
3// numbl CLI, with no browser involved.
4//
5// npm run test:matlab (NUMBL_DIR defaults to ~/src/numbl)
6//
7// Every assertion here is a statement made in the prose, so if one of these
8// fails either the code or the text is wrong.
9import { join, dirname } from 'node:path'
10import { fileURLToPath } from 'node:url'
11import { runLab } from './runlab.mjs'
13const root = join(dirname(fileURLToPath(import.meta.url)), '..')
14const src = (p) => join(root, 'src', p)
16let failures = 0
17const check = (name, ok, detail = '') => {
18 console.log(`${ok ? ' ok ' : ' FAIL '} ${name}${detail ? ` ${detail}` : ''}`)
19 if (!ok) failures++
20}
21const near = (a, b, tol) => Math.abs(a - b) <= tol
22const last = (row) => {
23 for (let i = row.length - 1; i >= 0; i--) if (row[i] != null) return row[i]
24 return null
25}
26const best = (row) => Math.min(...row.filter((v) => v != null))
28const ONE = src('chapters/onedim/driver.m')
29const TWO = src('chapters/twodim/driver.m')
31console.log('\nChapter 1: three-point differences')
32{
33 const { out } = runLab(ONE, src('chapters/ch01/fd2.m'), {
34 kind: 'smooth',
35 mode: 'sweep',
36 ns: [16, 32, 64, 128, 256, 512],
37 nmaxs: [1e9],
38 timed: 0,
39 })
40 const q = last(out.orders[0])
41 check('observed order is 2', near(q, 2, 0.02), `q = ${q.toFixed(3)}`)
42 const { out: low } = runLab(ONE, src('chapters/ch01/fd2.m'), {
43 kind: 'lowreg',
44 mode: 'sweep',
45 ns: [64, 128, 256, 512],
46 nmaxs: [1e9],
47 timed: 0,
48 })
49 const ql = last(low.orders[0])
50 check('the low-regularity solution converges more slowly', ql < 1.7, `q = ${ql.toFixed(3)}`)
51}
53console.log('\nChapter 2: the roundoff floor and Richardson')
54{
55 const { out } = runLab(ONE, src('chapters/ch01/fd2.m'), {
56 kind: 'smooth',
57 mode: 'sweep',
58 ns: [1024, 4096, 16384, 65536, 262144],
59 nmaxs: [1e9],
60 timed: 0,
61 })
62 const e = out.E[0]
63 const bi = e.indexOf(best(e))
64 check('the error stops improving before N = 262144', bi < e.length - 1, `best at N = ${out.ns[bi]}`)
65 check('the best error is a few times 1e-9', best(e) < 3e-8 && best(e) > 1e-10, `${best(e).toExponential(2)}`)
67 const { out: r } = runLab(ONE, src('chapters/ch02/richardson.m'), {
68 kind: 'smooth',
69 mode: 'sweep',
70 ns: [32, 64, 128, 256, 512],
71 nmaxs: [1e9, 1e9],
72 timed: 0,
73 })
74 check('Richardson gives order 4', near(last(r.orders[1]), 4, 0.1), `q = ${last(r.orders[1]).toFixed(3)}`)
75 check('and it beats the plain scheme by >1e4', last(r.E[0]) / last(r.E[1]) > 1e4)
77 const { out: rl } = runLab(ONE, src('chapters/ch02/richardson.m'), {
78 kind: 'lowreg',
79 mode: 'sweep',
80 ns: [64, 128, 256, 512],
81 nmaxs: [1e9, 1e9],
82 timed: 0,
83 })
84 check(
85 'on lowreg the extrapolation buys no extra order',
86 near(last(rl.orders[1]), last(rl.orders[0]), 0.1) && last(rl.orders[1]) < 1.7,
87 `q = ${last(rl.orders[1]).toFixed(2)} against ${last(rl.orders[0]).toFixed(2)}`,
88 )
89}
91console.log('\nChapter 3: the compact fourth-order scheme')
92{
93 const { out } = runLab(ONE, src('chapters/ch03/highorder.m'), {
94 kind: 'smooth',
95 mode: 'sweep',
96 ns: [64, 256, 1024],
97 nmaxs: [1e9, 1e9],
98 timed: 0,
99 })
100 check('Numerov is fourth order', near(last(out.orders[1]), 4, 0.05), `q = ${last(out.orders[1]).toFixed(3)}`)
101 check('at N = 1024 it is below 1e-9', last(out.E[1]) < 1e-9, last(out.E[1]).toExponential(2))
102}
104console.log('\nChapter 4: Chebyshev collocation')
105{
106 const { out } = runLab(ONE, src('chapters/ch04/spectral.m'), {
107 kind: 'smooth',
108 mode: 'sweep',
109 ns: [16, 32, 64, 128],
110 nmaxs: [1e9, 1e9, 1e9],
111 timed: 0,
112 })
113 check('spectral reaches 1e-13 by N = 64', out.E[2][2] < 1e-12, out.E[2][2].toExponential(2))
114 const { out: low } = runLab(ONE, src('chapters/ch04/spectral.m'), {
115 kind: 'lowreg',
116 mode: 'sweep',
117 ns: [64, 128, 256],
118 nmaxs: [1e9, 1e9, 1e9],
119 timed: 0,
120 })
121 check('on lowreg it is no better than the three-point scheme', low.E[2][2] > low.E[0][2])
122}
124console.log('\nChapters 5 and 6: the five-point Laplacian, two solvers')
125{
126 const { out } = runLab(TWO, src('chapters/ch06/fastpoisson.m'), {
127 kind: 'expsin',
128 par: 1,
129 mode: 'sweep',
130 ns: [8, 16, 32, 128, 512],
131 nmaxs: [32, 512],
132 timed: 0,
133 })
134 check('the five-point scheme is second order', near(last(out.orders[1]), 2, 0.02), `q = ${last(out.orders[1]).toFixed(3)}`)
135 for (let j = 0; j < 3; j++) {
136 check(
137 `dense and transform agree at N = ${out.ns[j]}`,
138 Math.abs(out.E[0][j] - out.E[1][j]) < 1e-9 * out.E[0][j],
139 `${out.E[0][j].toExponential(4)} vs ${out.E[1][j].toExponential(4)}`,
140 )
141 }
142 check('at N = 512 the error is 4e-8', near(last(out.E[1]) / 3.99e-8, 1, 0.05), last(out.E[1]).toExponential(3))
143}
145console.log('\nChapter 7: the compact nine-point stencil and 2-D Chebyshev')
146{
147 const { out } = runLab(TWO, src('chapters/ch07/highorder2d.m'), {
148 kind: 'expsin',
149 par: 1,
150 mode: 'sweep',
151 ns: [8, 16, 32, 64],
152 nmaxs: [512, 512, 64],
153 timed: 0,
154 })
155 check('nine-point is sixth order on Laplace', near(out.orders[1][1], 6, 0.1), `q = ${out.orders[1][1].toFixed(2)}`)
156 check('nine-point reaches 1e-13 by N = 32', out.E[1][2] < 1e-13, out.E[1][2].toExponential(2))
157 check('Chebyshev reaches 1e-13 by N = 16', out.E[2][1] < 1e-13, out.E[2][1].toExponential(2))
158}
160console.log('\nChapter 8: iterative solvers')
161{
162 const { out } = runLab(src('chapters/ch08/driver.m'), src('chapters/ch08/iterative.m'), {
163 kind: 'expsin',
164 par: 1,
165 N: 64,
166 nit: 60,
167 tol: 1e-11,
168 })
169 const mg = out.R[4].filter((v) => v != null)
170 check('multigrid converges in fewer than 15 cycles', mg.length < 16, `${mg.length - 1} cycles`)
171 const mgErr = out.E[4].filter((v) => v != null)
172 check(
173 'and it lands on the discretisation error',
174 near(mgErr[mgErr.length - 1] / out.discErr, 1, 0.02),
175 `${mgErr[mgErr.length - 1].toExponential(3)} vs ${out.discErr.toExponential(3)}`,
176 )
177 const jac = out.R[0].filter((v) => v != null)
178 check('Jacobi is nowhere near after 60 sweeps', jac[jac.length - 1] > 1e-2, jac[jac.length - 1].toExponential(2))
179}
181console.log('\nChapter 9: the L-shaped domain')
182{
183 const NSM = [
184 [16, 32, 64],
185 [8, 16, 32],
186 ]
187 const { out: b } = runLab(src('chapters/ch09/driver.m'), src('chapters/ch09/lshape.m'), {
188 kind: 'both',
189 mode: 'sweep',
190 nsm: NSM,
191 })
192 check('the difference scheme loses second order', near(last(b.orders[0]), 1.37, 0.1), `q = ${last(b.orders[0]).toFixed(3)}`)
193 check('the harmonic basis reaches 1e-11', last(b.E[1]) < 1e-11, last(b.E[1]).toExponential(2))
195 const { out: s } = runLab(src('chapters/ch09/driver.m'), src('chapters/ch09/lshape.m'), {
196 kind: 'smooth',
197 mode: 'sweep',
198 nsm: NSM,
199 })
200 check('with smooth data the corner costs nothing', near(last(s.orders[0]), 2, 0.02), `q = ${last(s.orders[0]).toFixed(3)}`)
202 const { out: g } = runLab(src('chapters/ch09/driver.m'), src('chapters/ch09/lshape.m'), {
203 kind: 'sing',
204 mode: 'sweep',
205 nsm: NSM,
206 })
207 check('the singular function is in the basis exactly', g.E[1][0] < 1e-14, g.E[1][0].toExponential(2))
208}
210console.log('\nChapter 10: boundary integral equations')
211{
212 const { out } = runLab(src('chapters/ch10/driver.m'), src('chapters/ch10/bie.m'), {
213 a: 0.3,
214 nlobe: 5,
215 par: 2.2,
216 mode: 'sweep',
217 ns: [64, 128, 256, 512],
218 })
219 check('Nyström reaches machine precision by n = 512', last(out.E[0]) < 1e-14, last(out.E[0]).toExponential(2))
220 check('convergence is geometric, not algebraic', out.E[0][0] / out.E[0][1] > 100)
221 check(
222 'close evaluation is orders of magnitude worse',
223 last(out.Enear[0]) / last(out.E[0]) > 1e6,
224 `${last(out.Enear[0]).toExponential(2)} at 98% radius`,
225 )
226 check('fundamental solutions reach 1e-13', best(out.E[1]) < 1e-13, best(out.E[1]).toExponential(2))
227 check(
228 'and they have no close-evaluation problem',
229 best(out.Enear[1]) < 1e-12,
230 best(out.Enear[1]).toExponential(2),
231 )
232}
234console.log('\nChapter 11: the Pac-Man benchmark')
235{
236 // zeros in a ladder skip that resolution (or the whole method)
237 const { out } = runLab(src('chapters/ch11/driver.m'), src('chapters/ch11/pacman.m'), {
238 kind: 'both',
239 mode: 'sweep',
240 nsm: [
241 [16, 32, 64, 0],
242 [4, 8, 16, 32],
243 [4, 8, 0, 0],
244 [8, 16, 32, 0],
245 ],
246 })
247 const qfd = last(out.orders[0])
248 check('Shortley-Weller is corner-limited near order 4/3', near(qfd, 1.37, 0.12), `q = ${qfd.toFixed(3)}`)
249 check('the harmonic basis reaches 1e-13 with 97 unknowns', out.E[1][3] < 1e-13, out.E[1][3].toExponential(2))
250 check('and converges geometrically', out.E[1][1] / out.E[1][2] > 100)
251 const ratio = out.E[2][0] / out.E[2][1]
252 check(
253 'graded panels gain about 2^(8/3) per four levels of depth',
254 ratio > 4 && ratio < 9,
255 `ratio = ${ratio.toFixed(2)}`,
256 )
257 check('and reach 3e-5 by depth 8', out.E[2][1] < 3e-5, out.E[2][1].toExponential(2))
258 check('polar Chebyshev is algebraic on the singular data', out.E[3][2] > 1e-5 && out.E[3][2] < 3e-4, out.E[3][2].toExponential(2))
260 const { out: sm } = runLab(src('chapters/ch11/driver.m'), src('chapters/ch11/pacman.m'), {
261 kind: 'smooth',
262 mode: 'sweep',
263 nsm: [
264 [16, 32, 64, 0],
265 [0, 0, 0, 0],
266 [0, 0, 0, 0],
267 [6, 16, 40, 0],
268 ],
269 })
270 const qsw = last(sm.orders[0])
271 check('with smooth data Shortley-Weller is second order', near(qsw, 2, 0.15), `q = ${qsw.toFixed(3)}`)
272 check('and polar Chebyshev is geometric again, 1e-9 by N = 40', sm.E[3][2] < 1e-9, sm.E[3][2].toExponential(2))
274 const { out: sg } = runLab(src('chapters/ch11/driver.m'), src('chapters/ch11/pacman.m'), {
275 kind: 'sing',
276 mode: 'sweep',
277 nsm: [
278 [0, 0, 0, 0],
279 [4, 0, 0, 0],
280 [0, 0, 0, 0],
281 [0, 0, 0, 0],
282 ],
283 })
284 check('the singular function is in the basis exactly', sg.E[1][0] < 1e-13, sg.E[1][0].toExponential(2))
285}
287console.log(failures === 0 ? '\nall checks passed' : `\n${failures} check(s) failed`)
288process.exit(failures === 0 ? 0 : 1)