Default damping: every mode falls 60 dB in 2.5 s; examples at 2.5 s; twin separation relative to recent displacement; stronger chaos strike
11 changed files+59−32
README.mdmodified+3−3View file
@@ -23,7 +23,7 @@ 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.
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. A script that sets neither `sys.friction` nor `sys.dashpots` gets a default under which every mode falls 60 dB in 2.5 s; `sys.friction = 0` turns damping off.
2727
2828 Bending stiffness (for plates and cymbals) is not in yet.
2929
@@ -49,7 +49,7 @@ The "small-strike reference" option runs the same system with every hit scaled d
4949
5050 ## Chaos test
5151
52-The "perturbed twin" option runs the system again with every hit changed by one part in 10^9 and plots, over time, the largest difference in any bead coordinate between the two runs, relative to the largest displacement. In chaotic motion this separation grows exponentially (a straight line on the log plot) until the runs are unrelated; in regular motion it grows slowly, roughly in proportion to time. The app fits the exponential rate over the stretch where the separation rises from 100 times the perturbation to 10^-2 and reports the e-folding time. Positions are compared at full precision, which needs the original run's positions at every frame in memory (limit 400 MB).
52+The "perturbed twin" option runs the system again with every hit changed by one part in 10^9 and plots, over time, the largest difference in any bead coordinate between the two runs, relative to the largest displacement of the original run over the preceding 50 ms (so the measure stays meaningful as the motion decays). In chaotic motion this separation grows exponentially (a straight line on the log plot) until the runs are unrelated; in regular motion it grows slowly, roughly in proportion to time. The app fits the exponential rate over the stretch where the separation rises from 100 times the perturbation to 10^-2 and reports the e-folding time. Positions are compared at full precision, which needs the original run's positions at every frame in memory (limit 400 MB).
5353
5454 ## Examples
5555
@@ -58,7 +58,7 @@ The "perturbed twin" option runs the system again with every hit changed by one
5858 | `drum` | A 15 cm drum head (823 beads, fundamental 237 Hz) struck by a soft mallet. Nearly linear at first (2% from the small-strike reference over 10 ms); the slight pitch change at large amplitude accumulates to 25% by 0.3 s. About 7x real time. |
5959 | `cascade` | The same net with a 0.5% pre-stretch, weaker dashpots and a hard strike. The share of energy in modes above 1 kHz grows from 4% (small strike) to about 25% within 30 ms. About 10x real time. |
6060 | `string` | 120 beads on a string, struck slightly off-axis; fast. |
61-| `chaos` | The minimal chaotic example: two beads on a short string, free to move in 3D (6 degrees of freedom), struck hard and off-axis. With the perturbed twin, the separation e-folds every 7 ms (10 ms at 0.3 of the strike, 105 ms at 0.1); at 0.05 of the strike the motion is regular and the separation stays near 4e-8 for 3 s. Runs in a fraction of a second. |
61+| `chaos` | The minimal chaotic example: two beads on a short string, free to move in 3D (6 degrees of freedom), struck hard and off-axis, with every mode decaying 60 dB in 2.5 s. With the perturbed twin, the separation e-folds every 2 ms and the runs are unrelated within 0.15 s; at 0.025 of the strike the motion is regular and the separation stays below 1e-8 for 3 s. With half the strike, the damping pulls the motion out of the chaotic range before the runs fully separate (transient chaos). Runs in a fraction of a second. |
6262
6363 Timings are for a laptop (Chrome, single thread).
6464
matlab/examples/cascade.mmodified+1−1View file
@@ -23,7 +23,7 @@ sys.pos = pos;
2323 sys.mass = sigma * (sqrt(3)/2) * a^2;
2424 sys.pinned = rim;
2525 sys.springs = bp_springs(pos, E, k, f0);
26-[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 3; 2000 1.6]);
26+[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 2.5; 2000 1.6]);
2727
2828 hit.t = 0;
2929 hit.beads = bp_near(pos, [0.05 0 0], 0.03);
matlab/examples/chaos.mmodified+5−5View file
@@ -3,7 +3,7 @@
33 % off-axis, the motion is chaotic: turn on "perturbed twin" to run it again
44 % with the strike changed by one part in 1e9 and watch the two runs
55 % separate exponentially until they are unrelated. Scale the strike down
6-% (try 0.05) and the motion is regular: the separation grows only slowly.
6+% (try 0.025) and the motion is regular: the separation grows only slowly.
77 %
88 % Transverse modes at small amplitude: sqrt(f0/(m*a))/(2*pi) * [1, sqrt(3)]
99 % = 159 and 276 Hz.
@@ -18,12 +18,12 @@ sys.pos = pos;
1818 sys.mass = m;
1919 sys.pinned = [1 4];
2020 sys.springs = bp_springs(pos, E, f0 / (strain * a), f0);
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]);
21+% every mode falls 60 dB in 2.5 s
22+[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 2.5]);
2323
24-s = 1; % strike strength; try 0.05 for regular motion
24+s = 1; % strike strength; try 0.025 for regular motion
2525 hit.beads = 2;
26-hit.velocity = s * [6 3 20]; % m/s, off-axis
26+hit.velocity = s * [12 6 40]; % m/s, off-axis
2727 sys.hits = hit;
2828
2929 sys.listen = struct('beads', 3, 'dir', [0 0 1]);
matlab/examples/drum.mmodified+2−2View file
@@ -21,10 +21,10 @@ 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);
2323 % 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.
24+% frequencies: low modes ring about 2.5 s, modes near 2 kHz about 0.75 s.
2525 % bp_decay turns this into sys.friction and sys.dashpots; the page plots
2626 % the decay time each mode gets.
27-[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 3; 2000 0.75]);
27+[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 2.5; 2000 0.75]);
2828
2929 % A soft mallet: a 1 ms raised-cosine push on the beads within 3 cm.
3030 hit.t = 0;
matlab/examples/string.mmodified+1−1View file
@@ -19,7 +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.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 7; 2000 3]); % T60 (s)
22+[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 2.5; 2000 1.5]); % T60 (s)
2323
2424 hit.t = 0;
2525 hit.beads = bp_near(pos, [0.12 0 0], 0.01);
src/App.tsxmodified+1−0View file
@@ -384,6 +384,7 @@ Elements
384384 sys.springs rows [i j k L0] or [i j k L0 k3]; force k*e + k3*e^3, e = length - L0
385385 sys.dashpots rows [i j c]; force c*(v_j - v_i) on bead i
386386 sys.friction scalar or n-by-1; force -c0*v
387+ (if neither friction nor dashpots is given: every mode falls 60 dB in 2.5 s)
387388 [sys.friction, sys.dashpots] = bp_decay(pos, mass, sys.springs, [f1 T1; f2 T2])
388389 sets both from decay times (T60, s) at two frequencies (Hz)
389390 sys.pinned bead indices (or logical mask) held fixed
src/sim/parse.tsmodified+10−1View file
@@ -8,6 +8,8 @@ type Mat = { rows: number; cols: number; data: Float64Array };
88
99 class SysError extends Error {}
1010
11+export const DEFAULT_T60 = 2.5; // s
12+
1113 function asMat(v: RuntimeValue | undefined, name: string): Mat | null {
1214 if (v === undefined) return null;
1315 if (typeof v === "number") return { rows: 1, cols: 1, data: Float64Array.of(v) };
@@ -159,7 +161,14 @@ function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warn
159161 dc[d] = col(2);
160162 if (dc[d] < 0) throw new SysError(`sys.dashpots row ${d + 1}: c must be nonnegative`);
161163 }
162- const friction = perBead(f.get("friction"), "sys.friction", n, 0);
164+ let friction = perBead(f.get("friction"), "sys.friction", n, 0);
165+ if (!f.has("friction") && !f.has("dashpots")) {
166+ // Default damping: every mode falls 60 dB in DEFAULT_T60 seconds
167+ // (amplitude rate sigma = ln(1000)/T60, from friction 2 sigma m per bead).
168+ const sigma = Math.log(1000) / DEFAULT_T60;
169+ friction = mass.map(m => 2 * sigma * m);
170+ warnings.push(`No damping given: using the default, every mode falls 60 dB in ${DEFAULT_T60} s. Set sys.friction = 0 for none.`);
171+ }
163172
164173 const hits: Hit[] = fields(f.get("hits"), "sys.hits").map((h, q) => {
165174 const nm = `sys.hits(${q + 1})`;
src/sim/twin.tsmodified+10−4View file
@@ -1,6 +1,7 @@
11 // Perturbed twin: rerun with every hit scaled by (1 + eps) and measure, at
22 // each frame, the largest difference in any bead coordinate, relative to the
3-// largest displacement of the original run.
3+// largest displacement of the original run over the preceding 50 ms (so that
4+// the measure stays meaningful while the motion decays).
45
56 import { simulate } from "./solver";
67 import type { BeadSystem } from "./types";
@@ -12,8 +13,13 @@ export function runTwin(
1213 eps: number,
1314 onProgress?: (t: number) => void
1415 ): { audio: Float32Array; sep: Float64Array } {
15- let scale = 0;
16- for (const x of saved) for (let q = 0; q < x.length; q++) scale = Math.max(scale, Math.abs(x[q] - sys.pos0[q]));
16+ const disp = saved.map(x => {
17+ let d = 0;
18+ for (let q = 0; q < x.length; q++) d = Math.max(d, Math.abs(x[q] - sys.pos0[q]));
19+ return d;
20+ });
21+ const win = Math.max(1, Math.round(0.05 * sys.frameRate));
22+ const scale = disp.map((_, f) => Math.max(...disp.slice(Math.max(0, f - win), f + 1)));
1723 const sep = new Float64Array(saved.length);
1824 const f = 1 + eps;
1925 const r = simulate(
@@ -34,7 +40,7 @@ export function runTwin(
3440 if (!y) return;
3541 let d = 0;
3642 for (let q = 0; q < x.length; q++) d = Math.max(d, Math.abs(x[q] - y[q]));
37- sep[fr] = scale > 0 ? d / scale : 0;
43+ sep[fr] = scale[fr] > 0 ? d / scale[fr] : 0;
3844 },
3945 }
4046 );
src/ui/Plots.tsxmodified+1−1View file
@@ -297,7 +297,7 @@ export function Separation({ p }: { p: RunPayload }) {
297297 <>
298298 <canvas ref={ref} className="plot short" />
299299 <div className="legend">
300- Largest difference in any bead coordinate between the run and its twin, relative to the largest displacement. A straight rise
300+ Largest difference in any bead coordinate between the run and its twin, relative to the largest displacement over the preceding 50 ms. A straight rise
301301 means exponential growth (chaos); regular motion drifts apart slowly. It levels off near 1 when the runs are unrelated.
302302 </div>
303303 </>
tests/solver.test.tsmodified+16−2View file
@@ -82,7 +82,7 @@ describe("string", () => {
8282 });
8383
8484 describe("energy", () => {
85- const { sys } = build(STRING("sys.hits.velocity = [0 0.3 2];"));
85+ const { sys } = build(STRING("sys.hits.velocity = [0 0.3 2]; sys.friction = 0;"));
8686 const drift = (method: "verlet" | "yoshida4", dt: number) =>
8787 simulate({ ...sys, method, dt, duration: 0.05 }).energyDrift;
8888
@@ -190,7 +190,7 @@ describe("chaos example", () => {
190190 });
191191
192192 it("is regular when struck softly", () => {
193- expect(rate(0.05)).toBeNull();
193+ expect(rate(0.025)).toBeNull();
194194 });
195195 });
196196
@@ -244,3 +244,17 @@ sys.listen.dir = [0 0 1];
244244 expect(() => build(STRING("[f, d] = bp_decay(pos, m, sys.springs, [100 0.1; 1000 0.5]);"))).toThrow(/must not ring longer/);
245245 });
246246 });
247+
248+describe("default damping", () => {
249+ it("makes every mode fall 60 dB in 2.5 s when the script sets none", () => {
250+ const { sys, warnings } = build(STRING(""));
251+ expect(warnings.join(" ")).toMatch(/default/);
252+ const modes = computeModes(sys, numblEig);
253+ for (const s of modalDamping(sys, modes)) expect(Math.log(1000) / s).toBeCloseTo(2.5, 9);
254+ });
255+
256+ it("is off when the script sets friction = 0", () => {
257+ const { sys } = build(STRING("sys.friction = 0;"));
258+ expect(Array.from(sys.friction).every(c => c === 0)).toBe(true);
259+ });
260+});
tests/twin.tsmodified+9−12View file
@@ -1,21 +1,18 @@
1-// Separation of two runs whose hits differ by a factor (1 + eps).
2-// Usage: tsx tests/twin.ts script.m duration [eps]
1+// Twin separation summary. Usage: tsx tests/twin.ts script.m duration [eps]
32 import { readFileSync, readdirSync } from "fs";
43 import { runScript } from "../src/sim/runScript";
54 import { simulate } from "../src/sim/solver";
5+import { runTwin } from "../src/sim/twin";
6+import { growthRate } from "../src/sim/analysis";
67 const helpers = readdirSync("matlab/helpers").map(name => ({ name, source: readFileSync("matlab/helpers/" + name, "utf8") }));
7-const { sys } = runScript(readFileSync(process.argv[2], "utf8"), helpers, s => process.stdout.write(s));
8+const { sys } = runScript(readFileSync(process.argv[2], "utf8"), helpers, () => {});
89 sys.duration = Number(process.argv[3]);
910 const eps = Number(process.argv[4] ?? 1e-9);
10-sys.frameRate = 1000;
11-const xs: Float64Array[] = [];
12-const a = simulate(sys, { onFrame: (_f, _t, x) => xs.push(x.slice()) });
13-const sep: number[] = [];
14-let scale = 0;
15-for (const x of xs) for (let q = 0; q < x.length; q++) scale = Math.max(scale, Math.abs(x[q] - sys.pos0[q]));
16-simulate({ ...sys, hits: sys.hits.map(h => ({ ...h, velocity: h.velocity.map(v => v * (1 + eps)) as any, force: h.force.map(v => v * (1 + eps)) as any })) },
17- { onFrame: (f, _t, x) => { let d = 0; for (let q = 0; q < x.length; q++) d = Math.max(d, Math.abs(x[q] - xs[f][q])); sep.push(d / scale); } });
11+const saved: Float64Array[] = [];
12+const a = simulate(sys, { onFrame: (_f, _t, x) => saved.push(x.slice()) });
13+const { sep } = runTwin(sys, saved, a.dt, eps);
1814 const out: string[] = [];
1915 const step = Math.max(1, Math.floor(sep.length / 14));
2016 for (let i = 0; i < sep.length; i += step) out.push(`${a.frameTimes[i].toFixed(2)}:${Math.log10(sep[i] + 1e-30).toFixed(1)}`);
21-console.log(`dt ${a.dt.toExponential(2)} drift ${a.energyDrift.toExponential(1)} ${(a.elapsedMs/1000).toFixed(1)}s | log10 sep: ${out.join(" ")}`);
17+const g = growthRate(a.frameTimes, sep, eps);
18+console.log(`rate ${g?.toFixed(0) ?? "-"} drift ${a.energyDrift.toExponential(1)} | log10 sep: ${out.join(" ")}`);