Damping by decay times: bp_decay helper, per-mode decay-time plot, examples converted
12 changed files+203−16
README.mdmodified+3−1View file
@@ -23,9 +23,11 @@ Each bead has a mass and a position in 3D. The script defines:
2323 - *Hits*: an instant velocity given to a set of beads, or a raised-cosine force pulse of given peak and duration. Any number, at any times.
2424 - A *listener*: a bead (or the average of several) and a direction; the displacement along that direction is the audio signal.
2525
26+Damping is easiest to set by its effect. `[sys.friction, sys.dashpots] = bp_decay(pos, mass, springs, [f1 T1; f2 T2])` chooses them so that every transverse mode of a flat net falls 60 dB in time T(f) following the law sigma(f) = sigma0 + sigma1 (2 pi f)^2 for the amplitude decay rate, through the two given points (one point gives the same decay time at every frequency). It does this by giving each bead friction 2 sigma0 m and each spring a parallel dashpot 2 sigma1 times its tension per length, which makes the damping matrix of the transverse motion 2 sigma0 M + 2 sigma1 K; the modes then stay uncoupled and the law holds exactly at small amplitude. Large motion and the in-plane modes are damped by the same elements, but not by this law. The page plots the small-amplitude decay time of every mode for whatever damping the script sets.
27+
2628 Bending stiffness (for plates and cymbals) is not in yet.
2729
28-The `matlab/helpers/bp_*.m` functions build common geometries (`bp_hexdisk`: a triangular lattice clipped to a disk, with its rim beads; `bp_chain`: a line of beads), list edges (`bp_edges`), make springs with a given tension at rest (`bp_springs`) and find beads near a point (`bp_near`). They are ordinary scripts; the app lists them, and a script can ignore them and build `sys` directly. The full list of `sys` fields is in the app under "sys fields".
30+The `matlab/helpers/bp_*.m` functions build common geometries (`bp_hexdisk`: a triangular lattice clipped to a disk, with its rim beads; `bp_chain`: a line of beads), list edges (`bp_edges`), make springs with a given tension at rest (`bp_springs`), set damping from decay times (`bp_decay`) and find beads near a point (`bp_near`). They are ordinary scripts; the app lists them, and a script can ignore them and build `sys` directly. The full list of `sys` fields is in the app under "sys fields".
2931
3032 ## Time integration
3133
matlab/examples/cascade.mmodified+5−5View file
@@ -3,9 +3,10 @@
33 %
44 % Compared with drum.m, the springs are stretched only 0.5% at rest, so the
55 % stretching caused by the motion itself is large relative to the rest
6-% tension, and the dashpots are 3x weaker, so high modes keep their energy
7-% long enough to pass it on. Turn on "small-strike reference" to hear and
8-% measure the difference from the same strike at small amplitude.
6+% tension, and the high modes ring about twice as long (T60 1.6 s rather
7+% than 0.75 s near 2 kHz), so they keep their energy long enough to pass
8+% it on. Turn on "small-strike reference" to hear and measure the
9+% difference from the same strike at small amplitude.
910
1011 R = 0.15;
1112 a = 0.01;
@@ -22,8 +23,7 @@ sys.pos = pos;
2223 sys.mass = sigma * (sqrt(3)/2) * a^2;
2324 sys.pinned = rim;
2425 sys.springs = bp_springs(pos, E, k, f0);
25-sys.dashpots = [E, 3e-5 * ones(size(E, 1), 1)];
26-sys.friction = 1e-4;
26+[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 3; 2000 1.6]);
2727
2828 hit.t = 0;
2929 hit.beads = bp_near(pos, [0.05 0 0], 0.03);
matlab/examples/chaos.mmodified+2−1View file
@@ -18,7 +18,8 @@ sys.pos = pos;
1818 sys.mass = m;
1919 sys.pinned = [1 4];
2020 sys.springs = bp_springs(pos, E, f0 / (strain * a), f0);
21-sys.friction = 2e-5; % light damping: decay time about 10 s
21+% light damping: every mode falls 60 dB in 60 s
22+[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 60]);
2223
2324 s = 1; % strike strength; try 0.05 for regular motion
2425 hit.beads = 2;
matlab/examples/drum.mmodified+5−2View file
@@ -20,8 +20,11 @@ sys.pos = pos;
2020 sys.mass = sigma * (sqrt(3)/2) * a^2; % lattice area per bead
2121 sys.pinned = rim;
2222 sys.springs = bp_springs(pos, E, k, f0);
23-sys.dashpots = [E, 1e-4 * ones(size(E, 1), 1)]; % damping that grows with frequency
24-sys.friction = 1e-4; % damping that does not
23+% Damping, set by decay times (T60: seconds to fall 60 dB) at two
24+% frequencies: low modes ring about 3 s, modes near 2 kHz about 0.75 s.
25+% bp_decay turns this into sys.friction and sys.dashpots; the page plots
26+% the decay time each mode gets.
27+[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 3; 2000 0.75]);
2528
2629 % A soft mallet: a 1 ms raised-cosine push on the beads within 3 cm.
2730 hit.t = 0;
matlab/examples/string.mmodified+1−2View file
@@ -19,8 +19,7 @@ sys.pos = pos;
1919 sys.mass = mu * a;
2020 sys.pinned = [1, N + 2];
2121 sys.springs = bp_springs(pos, E, k, T);
22-sys.dashpots = [E, 2e-4 * ones(size(E, 1), 1)];
23-sys.friction = 2e-5;
22+[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 7; 2000 3]); % T60 (s)
2423
2524 hit.t = 0;
2625 hit.beads = bp_near(pos, [0.12 0 0], 0.01);
matlab/helpers/bp_decay.madded+49−0View file
@@ -0,0 +1,49 @@
1+function [friction, dashpots, law] = bp_decay(pos, mass, springs, T60)
2+% BP_DECAY Friction and dashpots that give chosen decay times.
3+% [friction, dashpots] = bp_decay(pos, mass, springs, T60) returns values
4+% for sys.friction and sys.dashpots such that, at small amplitude, each
5+% transverse mode of a flat net decays at the amplitude rate
6+% sigma(f) = sigma0 + sigma1 * (2*pi*f)^2,
7+% passing through the decay times in T60 = [f1 T1; f2 T2] (Hz, seconds to
8+% fall by 60 dB). With a single row [f1 T1], every mode gets the same
9+% decay time. T1 must be at least T2 when f1 < f2 (higher modes cannot
10+% ring longer under this law).
11+%
12+% Each bead gets friction 2*sigma0*m, and each spring gets a dashpot
13+% 2*sigma1*g in parallel, where g is the spring's tension divided by its
14+% length at rest. The damping is then 2*sigma0*M + 2*sigma1*K for the
15+% transverse motion, so the modes stay uncoupled and the law is exact for
16+% them. Larger motion, and the in-plane modes, are damped by the same
17+% friction and dashpots but not by this law.
18+%
19+% law = [sigma0 sigma1] is the third output.
20+n = size(pos, 1);
21+s = log(1000) ./ T60(:, 2); % amplitude decay rates (1/s)
22+w2 = (2 * pi * T60(:, 1)).^2;
23+if size(T60, 1) == 1
24+ sigma0 = s(1);
25+ sigma1 = 0;
26+elseif size(T60, 1) == 2
27+ sigma1 = (s(2) - s(1)) / (w2(2) - w2(1));
28+ sigma0 = s(1) - sigma1 * w2(1);
29+else
30+ error('bp_decay: T60 must be [f1 T1] or [f1 T1; f2 T2]');
31+end
32+if sigma1 < 0
33+ error('bp_decay: the higher frequency must not ring longer than the lower one');
34+end
35+if sigma0 < 0
36+ error('bp_decay: decay at the higher frequency is too fast for sigma0 + sigma1*w^2; raise its T60 or lower the other');
37+end
38+friction = 2 * sigma0 * mass(:) .* ones(n, 1);
39+i = springs(:, 1);
40+j = springs(:, 2);
41+L = vecnorm(pos(i, :) - pos(j, :), 2, 2);
42+e = L - springs(:, 4);
43+f = springs(:, 3) .* e;
44+if size(springs, 2) > 4
45+ f = f + springs(:, 5) .* e.^3;
46+end
47+dashpots = [i, j, 2 * sigma1 * f ./ L];
48+law = [sigma0, sigma1];
49+end
src/App.tsxmodified+9−1View file
@@ -2,7 +2,7 @@ import { useCallback, useEffect, useRef, useState } from "react";
22 import { wavBlob } from "./sim/analysis";
33 import { BeadView } from "./ui/BeadView";
44 import { Editor } from "./ui/Editor";
5-import { Energy, ModeEnergy, Separation, Spectrogram, Spectrum, Waveform } from "./ui/Plots";
5+import { DecayTimes, Energy, ModeEnergy, Separation, Spectrogram, Spectrum, Waveform } from "./ui/Plots";
66 import type { FromWorker, RunPayload, RunRequest } from "./worker";
77
88 const exampleSources = import.meta.glob("../matlab/examples/*.m", { query: "?raw", import: "default", eager: true }) as Record<string, string>;
@@ -280,6 +280,12 @@ export default function App() {
280280 <Spectrum key={`spectrum-${runId}`} p={payload} fmax={fmax} />
281281 <h3>Energy by linear mode (the cascade)</h3>
282282 <ModeEnergy key={`modes-${runId}`} p={payload} />
283+ {payload.modeFreqs && payload.modeDecay && (
284+ <>
285+ <h3>Decay time by mode</h3>
286+ <DecayTimes p={payload} />
287+ </>
288+ )}
283289 <h3>Energy</h3>
284290 <Energy p={payload} />
285291 <h3>Beads</h3>
@@ -378,6 +384,8 @@ Elements
378384 sys.springs rows [i j k L0] or [i j k L0 k3]; force k*e + k3*e^3, e = length - L0
379385 sys.dashpots rows [i j c]; force c*(v_j - v_i) on bead i
380386 sys.friction scalar or n-by-1; force -c0*v
387+ [sys.friction, sys.dashpots] = bp_decay(pos, mass, sys.springs, [f1 T1; f2 T2])
388+ sets both from decay times (T60, s) at two frequencies (Hz)
381389 sys.pinned bead indices (or logical mask) held fixed
382390 Excitation
383391 sys.hits struct or struct array, each with
src/sim/modes.tsmodified+22−0View file
@@ -64,3 +64,25 @@ export function computeModes(sys: BeadSystem, eig: SymEig): Modes {
6464 });
6565 return { freqs, psi, free };
6666 }
67+
68+// Small-amplitude decay rate of each transverse mode (amplitude, 1/s), from
69+// the friction and dashpots: sigma_k = (1/2) u^T C u with u the mode shape
70+// normalized so u^T M u = 1. Exact when the damping matrix is a combination
71+// of M and K (as bp_decay builds it); otherwise it neglects the coupling of
72+// modes by damping, which is small for light damping.
73+export function modalDamping(sys: BeadSystem, modes: Modes): Float64Array {
74+ const { n, mass, friction, di, dj, dc } = sys;
75+ const { free, psi, freqs } = modes;
76+ const nf = free.length;
77+ const u = new Float64Array(n);
78+ const out = new Float64Array(freqs.length);
79+ for (let k = 0; k < freqs.length; k++) {
80+ u.fill(0);
81+ for (let r = 0; r < nf; r++) u[free[r]] = psi[k * nf + r] / Math.sqrt(mass[free[r]]);
82+ let c = 0;
83+ for (let r = 0; r < nf; r++) c += friction[free[r]] * u[free[r]] ** 2;
84+ for (let d = 0; d < di.length; d++) c += dc[d] * (u[di[d]] - u[dj[d]]) ** 2;
85+ out[k] = c / 2;
86+ }
87+ return out;
88+}
src/sim/solver.tsmodified+1−1View file
@@ -183,7 +183,7 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
183183 for (let i = 0; i < n; i++) normA = Math.max(normA, row[i] * invm[i]);
184184 }
185185 const hasDamping = normA > 0;
186- const frictionOnly = nd === 0;
186+ const frictionOnly = dc.every(c => c === 0); // then exp(-tau A) is diagonal
187187 const term = new Float64Array(N3);
188188 const tmp = new Float64Array(N3);
189189 const fa = new Float64Array(n);
src/ui/Plots.tsxmodified+51−1View file
@@ -1,7 +1,7 @@
11 import { useEffect, useMemo, useRef, useState } from "react";
22 import { modeBands, spectrogram, spectrum } from "../sim/analysis";
33 import type { RunPayload } from "../worker";
4-import { drawAxes, heat, niceTicks, setupCanvas } from "./plot";
4+import { drawAxes, fmt, heat, niceTicks, setupCanvas } from "./plot";
55
66 function useCanvas(draw: (c: HTMLCanvasElement) => void, deps: unknown[]) {
77 const ref = useRef<HTMLCanvasElement>(null);
@@ -303,3 +303,53 @@ export function Separation({ p }: { p: RunPayload }) {
303303 </>
304304 );
305305 }
306+
307+export function DecayTimes({ p }: { p: RunPayload }) {
308+ const freqs = p.modeFreqs!;
309+ const sig = p.modeDecay!;
310+ const ref = useCanvas(
311+ c => {
312+ const f = setupCanvas(c);
313+ const T60 = Array.from(sig, s => (s > 0 ? Math.log(1000) / s : Infinity));
314+ const finite = T60.filter(Number.isFinite);
315+ const fx0 = Math.log10(freqs[0] * 0.9);
316+ const fx1 = Math.log10(freqs[freqs.length - 1] * 1.1);
317+ let ty0 = Math.floor(Math.log10(Math.min(...finite, 1)) * 2) / 2;
318+ let ty1 = Math.ceil(Math.log10(Math.max(...finite, 1)) * 2) / 2;
319+ if (ty1 - ty0 < 1) (ty0 -= 0.5), (ty1 += 0.5);
320+ const logTicks = (lo: number, hi: number) => {
321+ const out: number[] = [];
322+ for (let d = Math.floor(lo); d <= hi; d++) for (const m of [1, 2, 5]) if (d + Math.log10(m) >= lo && d + Math.log10(m) <= hi) out.push(d + Math.log10(m));
323+ return out;
324+ };
325+ drawAxes(f, [fx0, fx1], [ty0, ty1], {
326+ xlabel: "mode frequency (Hz)",
327+ ylabel: "T60 (s)",
328+ xticks: logTicks(fx0, fx1),
329+ yticks: logTicks(ty0, ty1),
330+ xfmt: v => fmt(Math.pow(10, v)),
331+ yfmt: v => fmt(Math.pow(10, v)),
332+ });
333+ const { ctx, x0, y0, pw, ph } = f;
334+ ctx.fillStyle = f.css("--series");
335+ for (let k = 0; k < freqs.length; k++) {
336+ if (!Number.isFinite(T60[k])) continue;
337+ const px = x0 + ((Math.log10(freqs[k]) - fx0) / (fx1 - fx0)) * pw;
338+ const py = y0 + ph - ((Math.log10(T60[k]) - ty0) / (ty1 - ty0)) * ph;
339+ ctx.fillRect(px - 1.5, py - 1.5, 3, 3);
340+ }
341+ },
342+ [p]
343+ );
344+ const undamped = Array.from(sig).every(s => s === 0);
345+ if (undamped) return <div className="note">No damping: every mode rings indefinitely.</div>;
346+ return (
347+ <>
348+ <canvas ref={ref} className="plot short" />
349+ <div className="legend">
350+ Time for each transverse mode to fall 60 dB at small amplitude, from sys.friction and sys.dashpots. Use bp_decay to set it from two
351+ (frequency, T60) pairs.
352+ </div>
353+ </>
354+ );
355+}
src/worker.tsmodified+3−1View file
@@ -1,7 +1,7 @@
11 // Web Worker: runs the numbl script, computes linear modes, and simulates.
22
33 import { growthRate, relDiff } from "./sim/analysis";
4-import { computeModes, isPlanar } from "./sim/modes";
4+import { computeModes, isPlanar, modalDamping } from "./sim/modes";
55 import { numblEig, runScript } from "./sim/runScript";
66 import { simulate } from "./sim/solver";
77 import { runTwin } from "./sim/twin";
@@ -28,6 +28,7 @@ export interface RunPayload {
2828 };
2929 result: SimResult;
3030 modeFreqs: Float64Array | null;
31+ modeDecay: Float64Array | null; // small-amplitude amplitude decay rate of each mode (1/s)
3132 notes: string[];
3233 check: { dt: number; windows: { t: number; rel: number }[] } | null;
3334 // the same run with every hit scaled by 1e-3, its output scaled back up
@@ -151,6 +152,7 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
151152 },
152153 result,
153154 modeFreqs: modes ? modes.freqs : null,
155+ modeDecay: modes ? modalDamping(sys, modes) : null,
154156 notes,
155157 check: checkOut,
156158 linear: linearOut,
tests/solver.test.tsmodified+52−1View file
@@ -1,6 +1,6 @@
11 import { readFileSync, readdirSync } from "fs";
22 import { describe, expect, it } from "vitest";
3-import { computeModes } from "../src/sim/modes";
3+import { computeModes, modalDamping } from "../src/sim/modes";
44 import { numblEig, runScript } from "../src/sim/runScript";
55 import { growthRate } from "../src/sim/analysis";
66 import { frequencyErrorCents, simulate, stabilityLimit } from "../src/sim/solver";
@@ -193,3 +193,54 @@ describe("chaos example", () => {
193193 expect(rate(0.05)).toBeNull();
194194 });
195195 });
196+
197+describe("bp_decay", () => {
198+ // String struck in the shape of mode p; its amplitude should decay at
199+ // sigma0 + sigma1 w_p^2, the law through T60 = 0.5 s at 100 Hz and 0.1 s at 1 kHz.
200+ const shaped = (p: number) =>
201+ build(
202+ STRING(`
203+sys.hits = [];
204+for i = 1:N
205+ sys.hits(i).t = 0;
206+ sys.hits(i).beads = i + 1;
207+ sys.hits(i).velocity = [0 0 1e-5 * sin(${p} * pi * i / (N + 1))];
208+end
209+[sys.friction, sys.dashpots, law] = bp_decay(pos, m, sys.springs, [100 0.5; 1000 0.1]);
210+sys.listen.beads = 2:N+1;
211+sys.listen.dir = [0 0 1];
212+`).replace("sys.listen = struct('beads', 4);", "")
213+ );
214+
215+ it("gives each mode the decay rate of the chosen law", () => {
216+ const s = [Math.log(1000) / 0.5, Math.log(1000) / 0.1];
217+ const w2 = [(2 * Math.PI * 100) ** 2, (2 * Math.PI * 1000) ** 2];
218+ const sigma1 = (s[1] - s[0]) / (w2[1] - w2[0]);
219+ const sigma0 = s[0] - sigma1 * w2[0];
220+ const { sys } = shaped(1);
221+ const modes = computeModes(sys, numblEig);
222+ const md = modalDamping(sys, modes);
223+ for (let k = 0; k < modes.freqs.length; k++) {
224+ const expected = sigma0 + sigma1 * (2 * Math.PI * modes.freqs[k]) ** 2;
225+ expect(Math.abs(md[k] - expected) / expected).toBeLessThan(1e-9);
226+ }
227+ for (const p of [1, 3]) {
228+ const { sys: sp } = shaped(p);
229+ const r = simulate({ ...sp, duration: 0.3, frameRate: 1000 });
230+ // energy decays at 2 sigma; compare energy at 0.05 s and 0.25 s (averaged over a period)
231+ const avg = (t: number) => {
232+ let e = 0;
233+ let c = 0;
234+ for (let i = 0; i < r.frameTimes.length; i++) if (Math.abs(r.frameTimes[i] - t) < 0.01) (e += r.energy[i]), c++;
235+ return e / c;
236+ };
237+ const rate = Math.log(avg(0.05) / avg(0.25)) / 0.2 / 2;
238+ const expected = sigma0 + sigma1 * (2 * Math.PI * modes.freqs[p - 1]) ** 2;
239+ expect(Math.abs(rate - expected) / expected).toBeLessThan(0.02);
240+ }
241+ });
242+
243+ it("rejects a law the dashpots cannot make", () => {
244+ expect(() => build(STRING("[f, d] = bp_decay(pos, m, sys.springs, [100 0.1; 1000 0.5]);"))).toThrow(/must not ring longer/);
245+ });
246+});