1% A plucked dulcimer string, and the air it sounds in. Two wave equations,
2% advanced together, both compiled to WebGPU kernels.
3%
4% The string carries transverse displacement u (metres) on ns nodes along x:
5%
6% u_tt = cs^2 u_xx - kap^2 u_xxxx - 2 sig0 u_t + 2 sig1 (u_xx)_t
7%
8% cs = 2 Ls f0 is the wave speed a string of length Ls needs to sound the
9% fundamental f0. The u_xxxx term is bending stiffness, which sharpens the
10% upper partials (inharmonicity B: partial n lands near n*f0*sqrt(1 + B n^2)).
11% sig0 is plain decay — 6.91/t60 makes the amplitude fall 60 dB in t60
12% seconds — and sig1 damps high frequencies faster than low, which is why a
13% plucked note starts bright and mellows as it rings. The pluck is an initial
14% condition: init below shapes the string into a triangle and releases it
15% from rest. Both ends are pinned by the `pin` mask.
16%
17% The air carries pressure p on the box grid:
18%
19% p_tt + 2 sig p_t = c^2 lapw(p, wall) + s
20%
21% lapw is the wall-masked Laplacian: the dulcimer body's shell enters as the
22% mask `wall` (0 in the shell, 1 in air), which makes the shell rigid — the
23% Neumann condition, not a fast material, so it costs no timestep. The
24% source s is the string, coupled two ways with an independent gain on each:
25%
26% gline the string radiates directly: its acceleration, spread along its
27% line of cells (a thin string is a poor radiator in reality, so
28% this is the idealized version of that);
29% gbridge the string's pull on the bridge drives a patch of air above the
30% top plate — the instrument's actual mechanism, minus the plate's
31% own resonances, which this model does not carry.
32%
33% Either gain at zero switches that route off entirely. The absolute source
34% strength is arbitrary (the equation is linear, and playback is normalized);
35% the bridge branch carries a cs^2/Ls factor so that equal gains are of
36% comparable loudness rather than one branch drowning the other.
38function [u, um, p, pm] = init(xs, npts, Ls, pluckpos, amp)
39 % A triangle peaked at the pluck point, `amp` metres high, released from
40 % rest (um = u). Written via abs because min(a, b) = (a + b - |a - b|)/2.
41 a = xs / (pluckpos * Ls);
42 b = (Ls - xs) / (Ls - pluckpos * Ls);
43 u = amp * 0.5 * (a + b - abs(a - b));
44 um = u;
45 p = zeros(npts, 1);
46 pm = zeros(npts, 1);
47end
49function [un, uold, pn, pold] = step(u, um, p, pm, c, sig, wall, lineprof, boardprof, pin, dt, Ls, f0, B, t60, sig1, gline, gbridge)
50 % -- the string ---------------------------------------------------------
51 cs = 2 * Ls * f0; % tension, expressed as a wave speed
52 kap2 = B * (cs * Ls / pi)^2; % stiffness, expressed as inharmonicity
53 sig0 = 6.91 / t60;
55 lu = dxx(u);
56 lum = dxx(um);
57 l4 = dxxxx(u);
58 d = sig0 * dt;
59 un = pin .* ((2*u - (1 - d)*um + (dt*dt)*((cs*cs)*lu - kap2*l4) + (2*sig1*dt)*(lu - lum)) ./ (1 + d));
61 % -- the coupling -------------------------------------------------------
62 % The string's acceleration, sampled at each air cell's own x; and the
63 % slope at the bridge (the pull of the tension on it), broadcast. The
64 % scene says where each acts, through lineprof and boardprof. cs^2/Ls
65 % puts the bridge branch on the same footing as the acceleration branch.
66 acc = (un - 2*u + um) / (dt*dt);
67 sl = spread(acc);
68 fb = bridge(un);
69 s = gline * (lineprof .* sl) - (gbridge * (cs*cs) / Ls) * (boardprof .* fb);
71 % -- the air ------------------------------------------------------------
72 lp = lapw(p, wall);
73 sd = sig * dt;
74 pn = wall .* ((2*p - (1 - sd) .* pm + (c*dt).^2 .* lp + (dt*dt) * s) ./ (1 + sd));
76 % This step's fields become the next step's history. Lines of their own,
77 % so each plans as the copy it is.
78 uold = u;
79 pold = p;
80end