d8098edAn interactive tutorial on the wave equation in acoustics: strings, drums and rooms, measured by earJeremy Magland 1// Figures for the introduction and Part I: the three settings, d'Alembert's
2// travelling halves, folding, the leapfrog stencil, the dispersion curves and
3// the delay loop. Every curve that looks computed is computed, here, from
4// the same closed forms the chapters use.
5import type { ReactElement } from 'react'
6import { Arrowheads, Figure, Label, Node, arrow } from './svg.tsx'
7import { faded, familyColor, ink, method } from '../plot/palette.ts'
9const pts = (xy: [number, number][]) => xy.map(([x, y], i) => `${i ? 'L' : 'M'}${x.toFixed(1)} ${y.toFixed(1)}`).join('')
11/** The three settings of the tutorial, with the boundary each one has. */
12export function Settings() {
13 // a plucked string
14 const s0 = 40
15 const s1 = 210
16 const sy = 110
17 const apex = s0 + 0.3 * (s1 - s0)
18 // a square membrane, drawn with the sign pattern of mode (2, 1)
19 const q0 = 282
20 const qs = 132
21 const qy = 44
22 const cells: ReactElement[] = []
23 const n = 12
24 for (let i = 0; i < n; i++) {
25 for (let j = 0; j < n; j++) {
26 const x = (j + 0.5) / n
27 const y = (i + 0.5) / n
28 const v = Math.sin(2 * Math.PI * x) * Math.sin(Math.PI * y)
29 const c = v > 0 ? method.boundary : method.fd2
30 cells.push(
31 <rect
32 key={`${i}-${j}`}
33 x={q0 + (j * qs) / n}
34 y={qy + (i * qs) / n}
35 width={qs / n + 0.4}
36 height={qs / n + 0.4}
37 fill={faded(c, 1 - Math.min(1, Math.abs(v)) * 0.75)}
38 />,
39 )
40 }
41 }
42 // a room in oblique projection
43 const r0 = 500
44 const ry = 72
45 const rw = 150
46 const rh = 96
47 const dx = 44
48 const dy = -30
49 const S: [number, number] = [r0 + 36, ry + 70]
50 const R: [number, number] = [r0 + 128, ry + 52]
51 const hit: [number, number] = [r0 + 84, ry + rh]
52 return (
53 <Figure
54 width={716}
55 height={236}
56 caption="The three settings. A string, fixed at both ends; a drumhead, clamped at its rim, shown here in one of its modes with the two signs of the displacement in two colours; and a room, whose walls reflect sound, with the direct path from a source to a listener and one reflection from the floor."
57 >
58 <Arrowheads id="set" colors={[ink.muted]} />
59 {/* string */}
60 <line x1={s0} y1={sy - 18} x2={s0} y2={sy + 18} stroke={ink.secondary} strokeWidth={3} />
61 <line x1={s1} y1={sy - 18} x2={s1} y2={sy + 18} stroke={ink.secondary} strokeWidth={3} />
62 <line x1={s0} y1={sy} x2={s1} y2={sy} stroke={ink.axis} strokeWidth={1} strokeDasharray="3 3" />
63 <path d={pts([[s0, sy], [apex, sy - 46], [s1, sy]])} fill="none" stroke={method.modal} strokeWidth={2.2} />
64 <Label x={(s0 + s1) / 2} y={196} size={12.5} fill={ink.primary}>
65 string
66 </Label>
67 <Label x={(s0 + s1) / 2} y={215} size={11.5} fill={ink.muted}>
68 u = 0 at both ends
69 </Label>
70 {/* drum */}
71 {cells}
72 <rect x={q0} y={qy} width={qs} height={qs} fill="none" stroke={ink.secondary} strokeWidth={2.5} />
73 <Label x={q0 + qs / 2} y={196} size={12.5} fill={ink.primary}>
74 drum
75 </Label>
76 <Label x={q0 + qs / 2} y={215} size={11.5} fill={ink.muted}>
77 u = 0 on the rim
78 </Label>
79 {/* room */}
80 <g stroke={ink.secondary} strokeWidth={1.4} fill="none">
81 <rect x={r0} y={ry} width={rw} height={rh} />
82 <path d={`M${r0} ${ry}l${dx} ${dy}h${rw}v${rh}l${-dx} ${-dy}`} />
83 <path d={`M${r0 + rw} ${ry}l${dx} ${dy}`} />
84 </g>
85 <path d={`M${r0} ${ry + rh}l${dx} ${dy}h${rw}`} fill="none" stroke={ink.axis} strokeWidth={1} strokeDasharray="3 3" />
86 <line x1={S[0]} y1={S[1]} x2={R[0]} y2={R[1]} stroke={method.boundary} strokeWidth={1.6} />
87 <path d={pts([S, hit, R])} fill="none" stroke={faded(method.boundary, 0.35)} strokeWidth={1.4} strokeDasharray="4 3" />
88 <Node x={S[0]} y={S[1]} r={4.5} fill={method.waveguide} />
89 <Node x={R[0]} y={R[1]} r={4.5} fill={ink.node} />
90 <Label x={S[0] - 8} y={S[1] - 10} size={11} fill={method.waveguide} anchor="end">
91 source
92 </Label>
93 <Label x={R[0] + 8} y={R[1] - 10} size={11} fill={ink.secondary} anchor="start">
94 listener
95 </Label>
96 <Label x={r0 + (rw + dx) / 2} y={196} size={12.5} fill={ink.primary}>
97 room
98 </Label>
99 <Label x={r0 + (rw + dx) / 2} y={215} size={11.5} fill={ink.muted}>
100 ∂p/∂n = 0 on rigid walls
101 </Label>
102 </Figure>
103 )
104}
106/** d'Alembert: the plucked shape as two halves travelling apart. */
107export function Dalembert() {
108 const beta = 0.3
109 const x0 = 70
110 const x1 = 646
111 const W = x1 - x0
112 const phi = (x: number) => {
113 let r = ((x % 2) + 2) % 2
114 let sg = 1
115 if (r > 1) {
116 r = 2 - r
117 sg = -1
118 }
119 return sg * Math.min(r / beta, (1 - r) / (1 - beta))
120 }
121 const rows = [0, 0.12, 0.3, 0.55, 1.0]
122 const rh = 64
123 const top = 18
124 const amp = 26
125 const curve = (f: (x: number) => number, y0: number) => {
126 const out: [number, number][] = []
127 for (let i = 0; i <= 200; i++) {
128 const x = i / 200
129 out.push([x0 + x * W, y0 - amp * f(x)])
130 }
131 return pts(out)
132 }
133 return (
134 <Figure
135 width={716}
136 height={top + rows.length * rh + 24}
137 caption="d'Alembert's solution for a string plucked at β = 0.3. The shape (solid) is the sum of two copies of the initial triangle at half height (dashed), one moving right and one moving left, each reflected inverted at the ends. After half a period (bottom row) the string is the mirror image of its starting shape, upside down."
138 >
139 {rows.map((ct, k) => {
140 const y0 = top + k * rh + rh / 2
141 return (
142 <g key={k}>
143 <line x1={x0} y1={y0} x2={x1} y2={y0} stroke={ink.grid} strokeWidth={1} />
144 <line x1={x0} y1={y0 - 22} x2={x0} y2={y0 + 22} stroke={ink.axis} strokeWidth={2} />
145 <line x1={x1} y1={y0 - 22} x2={x1} y2={y0 + 22} stroke={ink.axis} strokeWidth={2} />
146 <path d={curve((x) => 0.5 * phi(x - ct), y0)} fill="none" stroke={method.fd2} strokeWidth={1.2} strokeDasharray="4 3" />
147 <path d={curve((x) => 0.5 * phi(x + ct), y0)} fill="none" stroke={method.boundary} strokeWidth={1.2} strokeDasharray="4 3" />
148 <path d={curve((x) => 0.5 * (phi(x - ct) + phi(x + ct)), y0)} fill="none" stroke={method.modal} strokeWidth={2.2} />
149 <Label x={x0 - 12} y={y0} size={11.5} fill={ink.muted} anchor="end">
150 {`ct = ${ct}`}
151 </Label>
152 </g>
153 )
154 })}
155 <Label x={x0 + W / 2} y={top + rows.length * rh + 10} size={11.5} fill={ink.muted}>
156 right-moving half dashed blue, left-moving half dashed pink, their sum solid
157 </Label>
158 </Figure>
159 )
160}
162/** Partials above the Nyquist frequency fold back below it. */
163export function Folding() {
164 const x0 = 50
165 const x1 = 670
166 const fs = 1
167 const f0 = fs / 7.35
168 const fmax = 1.02 * fs
169 const sx = (f: number) => x0 + ((x1 - x0) * f) / fmax
170 const base = 140
171 const h = 118
172 const parts = Array.from({ length: 7 }, (_, i) => i + 1)
173 const amp = (n: number) => h / n
174 const fold = (f: number) => {
175 const r = f % fs
176 return r > fs / 2 ? fs - r : r
177 }
178 const ay = base + 26
179 return (
180 <Figure
181 width={716}
182 height={262}
183 caption="Aliasing. A note's partials (green) lie at whole multiples of its fundamental. Sampling at fs cannot tell a partial above the Nyquist frequency fs/2 from one at the same distance below it, so each partial above fs/2 is reflected about fs/2 (arcs) and heard at a frequency that is generally not a multiple of the fundamental (red)."
184 >
185 <line x1={x0} y1={base} x2={x1} y2={base} stroke={ink.axis} strokeWidth={1.2} />
186 <line x1={sx(fs / 2)} y1={ay + 70} x2={sx(fs / 2)} y2={14} stroke={ink.muted} strokeWidth={1.2} strokeDasharray="4 4" />
187 <Label x={sx(fs / 2)} y={8} size={11.5} fill={ink.muted}>
188 fs / 2
189 </Label>
190 <line x1={sx(fs)} y1={base + 6} x2={sx(fs)} y2={base - 6} stroke={ink.muted} strokeWidth={1.2} />
191 <Label x={sx(fs)} y={base + 16} size={11.5} fill={ink.muted}>
192 fs
193 </Label>
194 <Label x={sx(0)} y={base + 16} size={11.5} fill={ink.muted}>
195 0
196 </Label>
197 {parts.map((n) => {
198 const f = n * f0
199 const above = f > fs / 2
200 const g = fold(f)
201 const rx = (sx(f) - sx(g)) / 2
202 return (
203 <g key={n}>
204 <line
205 x1={sx(f)}
206 y1={base}
207 x2={sx(f)}
208 y2={base - amp(n)}
209 stroke={above ? faded(method.modal, 0.55) : method.modal}
210 strokeWidth={3}
211 />
212 <Label x={sx(f)} y={base + 16} size={10.5} fill={ink.muted}>
213 {n === 1 ? 'f₀' : `${n}f₀`}
214 </Label>
215 {above && (
216 <>
217 <line x1={sx(g)} y1={base} x2={sx(g)} y2={base - amp(n)} stroke={method.alt} strokeWidth={3} />
218 <path
219 d={`M${sx(f)} ${ay}A${rx} ${rx * 0.28} 0 0 1 ${sx(g)} ${ay}`}
220 fill="none"
221 stroke={faded(method.alt, 0.25)}
222 strokeWidth={1.3}
223 />
224 <line x1={sx(f)} y1={base + 22} x2={sx(f)} y2={ay} stroke={faded(method.alt, 0.25)} strokeWidth={1.3} />
225 <line x1={sx(g)} y1={base + 4} x2={sx(g)} y2={ay} stroke={faded(method.alt, 0.25)} strokeWidth={1.3} />
226 </>
227 )}
228 </g>
229 )
230 })}
231 </Figure>
232 )
233}
235/** The leapfrog stencil in space and time, and the two domains of dependence. */
236export function Stencil() {
237 const gx = 60
238 const gy = 30
239 const ox = 90
240 const oy = 190
241 const node = (j: number, n: number): [number, number] => [ox + j * gx, oy - n * gy * 1.5]
242 const lam = 0.6
243 const P = node(4, 3)
244 const nodes: ReactElement[] = []
245 for (let n = 0; n <= 3; n++) {
246 for (let j = 0; j <= 8; j++) {
247 const [x, y] = node(j, n)
248 nodes.push(<Node key={`${j}-${n}`} x={x} y={y} r={2.4} fill={ink.axis} />)
249 }
250 }
251 const used: [number, number, string][] = [
252 [3, 2, method.fd2],
253 [4, 2, method.fd2],
254 [5, 2, method.fd2],
255 [4, 1, method.fd2],
256 ]
257 const cone = (w: number) => pts([node(4 - 3 * w, 0), P, node(4 + 3 * w, 0)])
258 return (
259 <Figure
260 width={716}
261 height={236}
262 caption="The leapfrog stencil. The new value at the top (orange) is computed from three values at the previous step and one at the step before (blue). Information therefore moves at most one node per step, which bounds the numerical domain of dependence (grey); the physical one (dashed) spreads at c Δt per step, λ nodes. The scheme is stable only if the numerical domain contains the physical one, λ ≤ 1."
263 >
264 <path d={`${cone(1)}Z`} fill={faded(ink.muted, 0.85)} stroke={ink.muted} strokeWidth={1} />
265 <path d={cone(lam)} fill="none" stroke={method.modal} strokeWidth={1.6} strokeDasharray="5 4" />
266 {nodes}
267 {used.map(([j, n, c]) => {
268 const [x, y] = node(j, n)
269 return (
270 <g key={`${j}-${n}`}>
271 <line x1={x} y1={y} x2={P[0]} y2={P[1]} stroke={faded(method.fd2, 0.4)} strokeWidth={1.3} />
272 <Node x={x} y={y} r={5} fill={c} />
273 </g>
274 )
275 })}
276 <Node x={P[0]} y={P[1]} r={6} fill={method.fd4} />
277 <Label x={P[0] + 10} y={P[1] - 14} size={11.5} fill={method.fd4} anchor="start">
278 u(j, n+1)
279 </Label>
280 {[0, 1, 2, 3].map((n) => (
281 <Label key={n} x={ox - 16} y={node(0, n)[1]} size={11} fill={ink.muted} anchor="end">
282 {n === 3 ? 'n + 1' : n === 2 ? 'n' : n === 1 ? 'n − 1' : 'n − 2'}
283 </Label>
284 ))}
285 <Label x={ox + 4 * gx} y={oy + 20} size={11} fill={ink.muted}>
286 x, one node every h
287 </Label>
288 <Label x={596} y={70} size={11.5} fill={ink.secondary} anchor="start">
289 numerical:
290 </Label>
291 <Label x={596} y={86} size={11.5} fill={ink.secondary} anchor="start">
292 1 node per step
293 </Label>
294 <Label x={596} y={116} size={11.5} fill={method.modal} anchor="start">
295 physical:
296 </Label>
297 <Label x={596} y={132} size={11.5} fill={method.modal} anchor="start">
298 λ nodes per step
299 </Label>
300 </Figure>
301 )
302}
304/** Phase speed against wavenumber for several Courant numbers. */
305export function Dispersion() {
306 const x0 = 70
307 const x1 = 560
308 const y0 = 196
309 const y1 = 24
310 const vmin = 0.55
311 const sx = (t: number) => x0 + (x1 - x0) * t
312 const sy = (v: number) => y0 - ((y0 - y1) * (v - vmin)) / (1 - vmin)
313 const lams = [1, 0.75, 0.5, 0.25]
314 const curve = (lam: number) => {
315 const out: [number, number][] = []
316 for (let i = 1; i <= 200; i++) {
317 const kh = (Math.PI * i) / 200
318 const r = (2 * Math.asin(lam * Math.sin(kh / 2))) / (lam * kh)
319 out.push([sx(i / 200), sy(r)])
320 }
321 return pts(out)
322 }
323 const ticks = [0.55, 0.7, 0.85, 1]
324 const ppw = [20, 10, 5, 4, 3, 2]
325 return (
326 <Figure
327 width={716}
328 height={250}
329 caption="The leapfrog dispersion relation, as the speed of each wave on the grid relative to the true speed, against its wavenumber. At λ = 1 every wave travels at exactly the right speed; below it, short waves travel slowly, and at two points per wavelength (right edge) they travel at 2 arcsin(λ) / (λπ) of the right speed."
330 >
331 <line x1={x0} y1={y0} x2={x1} y2={y0} stroke={ink.axis} />
332 <line x1={x0} y1={y0} x2={x0} y2={y1 - 6} stroke={ink.axis} />
333 {ticks.map((v) => (
334 <g key={v}>
335 <line x1={x0} x2={x1} y1={sy(v)} y2={sy(v)} stroke={ink.grid} />
336 <Label x={x0 - 8} y={sy(v)} size={11} fill={ink.muted} anchor="end">
337 {v.toFixed(2)}
338 </Label>
339 </g>
340 ))}
341 {ppw.map((P) => (
342 <g key={P}>
343 <line x1={sx(2 / P)} x2={sx(2 / P)} y1={y0} y2={y0 + 5} stroke={ink.axis} />
344 <Label x={sx(2 / P)} y={y0 + 16} size={11} fill={ink.muted}>
345 {P}
346 </Label>
347 </g>
348 ))}
349 <Label x={x1} y={y0 + 34} size={11} fill={ink.muted} anchor="end">
350 points per wavelength
351 </Label>
352 <Label x={x0 - 44} y={y1 - 12} size={11} fill={ink.muted} anchor="start">
353 speed on the grid / true speed
354 </Label>
355 {lams.map((l, k) => (
356 <g key={l}>
357 <path d={curve(l)} fill="none" stroke={familyColor(k, lams.length)} strokeWidth={2} />
358 <Label
359 x={x1 + 10}
360 y={sy((2 * Math.asin(l)) / (l * Math.PI)) + [0, -4, 0, 9][k]}
361 size={11.5}
362 fill={familyColor(k, lams.length)}
363 anchor="start"
364 >
365 {`λ = ${l}`}
366 </Label>
367 </g>
368 ))}
369 </Figure>
370 )
371}
373/** The Karplus-Strong loop as a block diagram. */
374export function Loop() {
375 const y = 90
376 const box = (x: number, w: number, label: string, sub: string, c: string) => (
377 <g>
378 <rect x={x} y={y - 24} width={w} height={48} rx={6} fill={faded(c, 0.82)} stroke={c} strokeWidth={1.4} />
379 <Label x={x + w / 2} y={y - 6} size={12.5} fill={ink.primary}>
380 {label}
381 </Label>
382 <Label x={x + w / 2} y={y + 11} size={10.5} fill={ink.muted}>
383 {sub}
384 </Label>
385 </g>
386 )
387 return (
388 <Figure
389 width={716}
390 height={200}
391 caption="The Karplus–Strong string. A delay of P samples, a two-point average that damps each partial a little on every trip and delays the loop by half a sample, and a fractional delay that makes up the length the whole samples cannot. The output can be taken anywhere in the loop; the pitch is fs divided by the total delay."
392 >
393 <Arrowheads id="ks" colors={[ink.secondary]} />
394 <circle cx={70} cy={y} r={15} fill="none" stroke={ink.secondary} strokeWidth={1.4} />
395 <Label x={70} y={y + 1} size={15} fill={ink.secondary}>
396 +
397 </Label>
398 <Label x={70} y={y - 40} size={11} fill={ink.muted}>
399 pluck
400 </Label>
401 <line x1={70} y1={y - 32} x2={70} y2={y - 17} stroke={ink.secondary} strokeWidth={1.3} markerEnd={arrow('ks', 0)} />
402 <line x1={85} y1={y} x2={122} y2={y} stroke={ink.secondary} strokeWidth={1.3} markerEnd={arrow('ks', 0)} />
403 {box(124, 150, 'delay z⁻ᴾ', 'P whole samples', method.waveguide)}
404 <line x1={274} y1={y} x2={312} y2={y} stroke={ink.secondary} strokeWidth={1.3} markerEnd={arrow('ks', 0)} />
405 {box(314, 160, '½ (1 + z⁻¹)', 'loss, and ½ sample', method.fd2)}
406 <line x1={474} y1={y} x2={512} y2={y} stroke={ink.secondary} strokeWidth={1.3} markerEnd={arrow('ks', 0)} />
407 {box(514, 140, 'fractional delay', 'η of a sample', method.boundary)}
408 <path d={`M654 ${y}H690V${y + 62}H70V${y + 17}`} fill="none" stroke={ink.secondary} strokeWidth={1.3} markerEnd={arrow('ks', 0)} />
409 <line x1={690} y1={y} x2={690} y2={y - 42} stroke={ink.secondary} strokeWidth={1.3} markerEnd={arrow('ks', 0)} />
410 <Label x={690} y={y - 52} size={11} fill={ink.muted}>
411 out
412 </Label>
413 <Label x={380} y={y + 76} size={11} fill={ink.muted}>
414 total delay P + ½ + η = fs / f₀
415 </Label>
416 </Figure>
417 )
418}