Minimal chaos example and perturbed-twin chaos test with separation plot and growth-rate fit
10 changed files+254−8
README.mdmodified+5−0View file
@@ -45,6 +45,10 @@ For a net that rests flat (all beads at one z), small transverse motion decouple
4545
4646 The "small-strike reference" option runs the same system with every hit scaled down by 1000 and the output scaled back up, which is what the system would do if it were linear. The app plays it for comparison and reports the difference from the actual run.
4747
48+## Chaos test
49+
50+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).
51+
4852 ## Examples
4953
5054 | Example | What it shows |
@@ -52,6 +56,7 @@ The "small-strike reference" option runs the same system with every hit scaled d
5256 | `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. |
5357 | `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. |
5458 | `string` | 120 beads on a string, struck slightly off-axis; fast. |
59+| `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. |
5560
5661 Timings are for a laptop (Chrome, single thread).
5762
matlab/examples/chaos.madded+31−0View file
@@ -0,0 +1,31 @@
1+% A minimal chaotic system: two beads on a short string, pinned at both
2+% ends, free to move in 3D (6 degrees of freedom). Struck hard and
3+% off-axis, the motion is chaotic: turn on "perturbed twin" to run it again
4+% with the strike changed by one part in 1e9 and watch the two runs
5+% separate exponentially until they are unrelated. Scale the strike down
6+% (try 0.05) and the motion is regular: the separation grows only slowly.
7+%
8+% Transverse modes at small amplitude: sqrt(f0/(m*a))/(2*pi) * [1, sqrt(3)]
9+% = 159 and 276 Hz.
10+
11+a = 0.01; % bead spacing (m)
12+m = 1e-4; % bead mass (kg)
13+f0 = 1; % string tension (N)
14+strain = 0.05; % spring stretch at rest
15+
16+[pos, E] = bp_chain(4, 3 * a);
17+sys.pos = pos;
18+sys.mass = m;
19+sys.pinned = [1 4];
20+sys.springs = bp_springs(pos, E, f0 / (strain * a), f0);
21+sys.friction = 2e-5; % light damping: decay time about 10 s
22+
23+s = 1; % strike strength; try 0.05 for regular motion
24+hit.beads = 2;
25+hit.velocity = s * [6 3 20]; % m/s, off-axis
26+sys.hits = hit;
27+
28+sys.listen = struct('beads', 3, 'dir', [0 0 1]);
29+sys.duration = 3;
30+sys.dt = 1 / (48000 * 8); % the automatic step is set at rest; the hard strike
31+ % needs a finer one to keep energy drift small
src/App.tsxmodified+34−6View 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, Spectrogram, Spectrum, Waveform } from "./ui/Plots";
5+import { 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>;
@@ -41,7 +41,8 @@ export default function App() {
4141 const [viewTime, setViewTime] = useState<number | null>(null);
4242 const [dtCheck, setDtCheck] = useState(false);
4343 const [linearRef, setLinearRef] = useState(false);
44- const [playing, setPlaying] = useState<"main" | "linear" | null>(null);
44+ const [twin, setTwin] = useState(false);
45+ const [playing, setPlaying] = useState<"main" | "linear" | "twin" | null>(null);
4546 const [runId, setRunId] = useState(0); // remounts panels with per-run state
4647 const worker = useRef<Worker | null>(null);
4748 const audioCtx = useRef<AudioContext | null>(null);
@@ -93,10 +94,10 @@ export default function App() {
9394 setError(e.message || "Worker error");
9495 setRunning(false);
9596 };
96- const req: RunRequest = { type: "run", script, dtCheck, linear: linearRef };
97+ const req: RunRequest = { type: "run", script, dtCheck, linear: linearRef, twin };
9798 w.postMessage(req);
9899 },
99- [script, stopAudio, dtCheck, linearRef]
100+ [script, stopAudio, dtCheck, linearRef, twin]
100101 );
101102
102103 const stop = () => {
@@ -106,11 +107,16 @@ export default function App() {
106107 setStage("Stopped");
107108 };
108109
109- const play = (which: "main" | "linear") => {
110+ const play = (which: "main" | "linear" | "twin") => {
110111 if (!payload) return;
111112 stopAudio();
112113 const { sampleRate } = payload.result;
113- const audio = which === "linear" && payload.linear ? payload.linear.audio : payload.result.audio;
114+ const audio =
115+ which === "linear" && payload.linear
116+ ? payload.linear.audio
117+ : which === "twin" && payload.twin
118+ ? payload.twin.audio
119+ : payload.result.audio;
114120 const ctx = audioCtx.current ?? new AudioContext();
115121 audioCtx.current = ctx;
116122 // both signals share the main run's gain, so the A/B comparison keeps levels
@@ -196,6 +202,9 @@ export default function App() {
196202 <label title="Run again with every hit scaled by 1e-3 (then scaled back up): what the system would do if it were linear">
197203 <input type="checkbox" checked={linearRef} onChange={e => setLinearRef(e.target.checked)} /> small-strike reference
198204 </label>
205+ <label title="Run again with every hit changed by one part in 1e9 and track how fast the two runs separate (exponential growth means chaos)">
206+ <input type="checkbox" checked={twin} onChange={e => setTwin(e.target.checked)} /> perturbed twin
207+ </label>
199208 </div>
200209 <Editor value={script} onChange={setScript} onRun={() => !running && run()} />
201210 {running && (
@@ -237,6 +246,11 @@ export default function App() {
237246 {playing === "linear" ? "Stop" : "Play small-strike reference"}
238247 </button>
239248 )}
249+ {payload.twin && (
250+ <button onClick={() => (playing === "twin" ? stopAudio() : play("twin"))} title="The run with the strike changed by one part in 1e9">
251+ {playing === "twin" ? "Stop" : "Play twin"}
252+ </button>
253+ )}
240254 <button onClick={download}>Download WAV</button>
241255 <span className="grow" />
242256 <label>
@@ -254,6 +268,12 @@ export default function App() {
254268 </div>
255269 <h3>Listener signal</h3>
256270 <Waveform p={payload} cursor={cursor ?? viewTime} />
271+ {payload.twin && (
272+ <>
273+ <h3>Separation from the perturbed twin (chaos test)</h3>
274+ <Separation p={payload} />
275+ </>
276+ )}
257277 <h3>Spectrogram</h3>
258278 <Spectrogram p={payload} fmax={fmax} />
259279 <h3>Spectrum, two windows (detuning shows as peaks shifting)</h3>
@@ -304,6 +324,14 @@ function Report({ p }: { p: RunPayload }) {
304324 ))}
305325 </tbody>
306326 </table>
327+ {p.twin && (
328+ <div className="caption">
329+ Perturbed twin (hits changed by 1 part in 10^9):{" "}
330+ {p.twin.rate !== null
331+ ? `separation grows exponentially, e-folding time ${(1000 / p.twin.rate).toPrecision(2)} ms (rate ${p.twin.rate.toPrecision(3)} /s). Chaotic.`
332+ : `separation reached ${p.twin.sep.reduce((a, b) => Math.max(a, b), 0).toExponential(1)} of the motion's size; no exponential growth to 1e-2 within the run.`}
333+ </div>
334+ )}
307335 {p.linear && (
308336 <>
309337 <div className="caption">
src/sim/analysis.tsmodified+27−0View file
@@ -139,3 +139,30 @@ export function wavBlob(x: Float32Array, fs: number): Blob {
139139 for (let i = 0; i < x.length; i++) dv.setInt16(44 + 2 * i, Math.round(32767 * Math.max(-1, Math.min(1, g * x[i]))), true);
140140 return new Blob([buf], { type: "audio/wav" });
141141 }
142+
143+// Exponential growth rate of a separation curve: least-squares slope of
144+// ln(sep) over the stretch from 100*eps up to 1e-2 (relative to the motion's
145+// size). Returns null when the separation never reaches 1e-2 or the stretch is
146+// too short to fit.
147+export function growthRate(times: Float64Array, sep: Float64Array, eps: number): number | null {
148+ let i0 = -1;
149+ let i1 = -1;
150+ for (let i = 0; i < sep.length; i++) {
151+ if (i0 < 0 && sep[i] > 100 * eps) i0 = i;
152+ if (i0 >= 0 && sep[i] > 1e-2) {
153+ i1 = i;
154+ break;
155+ }
156+ }
157+ if (i0 < 0 || i1 < 0 || i1 - i0 < 4) return null;
158+ let st = 0, sy = 0, stt = 0, sty = 0;
159+ const n = i1 - i0 + 1;
160+ for (let i = i0; i <= i1; i++) {
161+ const y = Math.log(sep[i]);
162+ st += times[i];
163+ sy += y;
164+ stt += times[i] * times[i];
165+ sty += times[i] * y;
166+ }
167+ return (n * sty - st * sy) / (n * stt - st * st);
168+}
src/sim/solver.tsmodified+3−0View file
@@ -67,6 +67,8 @@ export function frequencyErrorCents(method: Method, omega: number, h: number): n
6767 export interface SimOptions {
6868 modes?: Modes | null;
6969 onProgress?: (t: number) => void;
70+ // called at every frame record with the full-precision state (do not keep the arrays)
71+ onFrame?: (frame: number, t: number, x: Float64Array, v: Float64Array) => void;
7072 maxAnimBytes?: number;
7173 }
7274
@@ -370,6 +372,7 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
370372 energy[frame] = kinetic() + pe - peRest;
371373 dissipatedRec[frame] = dissipated;
372374 recordModes(frame);
375+ opts.onFrame?.(frame, t, x, v);
373376 if (sample % animStride === 0) {
374377 const d = new Float32Array(N3);
375378 for (let q = 0; q < N3; q++) d[q] = x[q] - pos0[q];
src/sim/twin.tsadded+42−0View file
@@ -0,0 +1,42 @@
1+// Perturbed twin: rerun with every hit scaled by (1 + eps) and measure, at
2+// each frame, the largest difference in any bead coordinate, relative to the
3+// largest displacement of the original run.
4+
5+import { simulate } from "./solver";
6+import type { BeadSystem } from "./types";
7+
8+export function runTwin(
9+ sys: BeadSystem,
10+ saved: Float64Array[], // positions of the original run at each frame
11+ dt: number,
12+ eps: number,
13+ onProgress?: (t: number) => void
14+): { 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]));
17+ const sep = new Float64Array(saved.length);
18+ const f = 1 + eps;
19+ const r = simulate(
20+ {
21+ ...sys,
22+ dt,
23+ hits: sys.hits.map(h => ({
24+ ...h,
25+ velocity: [h.velocity[0] * f, h.velocity[1] * f, h.velocity[2] * f],
26+ force: [h.force[0] * f, h.force[1] * f, h.force[2] * f],
27+ })),
28+ },
29+ {
30+ maxAnimBytes: 0,
31+ onProgress,
32+ onFrame: (fr, _t, x) => {
33+ const y = saved[fr];
34+ if (!y) return;
35+ let d = 0;
36+ 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;
38+ },
39+ }
40+ );
41+ return { audio: r.audio, sep };
42+}
src/ui/Plots.tsxmodified+39−0View file
@@ -264,3 +264,42 @@ export function Spectrum({ p, fmax }: { p: RunPayload; fmax: number }) {
264264 </>
265265 );
266266 }
267+
268+export function Separation({ p }: { p: RunPayload }) {
269+ const tw = p.twin!;
270+ const ref = useCanvas(
271+ c => {
272+ const f = setupCanvas(c);
273+ const lo = Math.floor(Math.log10(tw.eps)) - 1;
274+ drawAxes(f, [0, p.sys.duration], [lo, 1], {
275+ xlabel: "time (s)",
276+ ylabel: "log10 separation",
277+ yticks: Array.from({ length: 2 - lo }, (_, i) => lo + i).filter(v => v % 2 === 0 || 1 - lo < 8),
278+ });
279+ const { ctx, x0, y0, pw, ph } = f;
280+ ctx.strokeStyle = f.css("--series");
281+ ctx.lineWidth = 1.5;
282+ ctx.beginPath();
283+ let started = false;
284+ for (let i = 0; i < tw.times.length; i++) {
285+ if (!(tw.sep[i] > 0)) continue;
286+ const px = x0 + (tw.times[i] / p.sys.duration) * pw;
287+ const py = y0 + ph - ((Math.log10(tw.sep[i]) - lo) / (1 - lo)) * ph;
288+ if (started) ctx.lineTo(px, py);
289+ else ctx.moveTo(px, py);
290+ started = true;
291+ }
292+ ctx.stroke();
293+ },
294+ [tw]
295+ );
296+ return (
297+ <>
298+ <canvas ref={ref} className="plot short" />
299+ <div className="legend">
300+ Largest difference in any bead coordinate between the run and its twin, relative to the largest displacement. A straight rise
301+ means exponential growth (chaos); regular motion drifts apart slowly. It levels off near 1 when the runs are unrelated.
302+ </div>
303+ </>
304+ );
305+}
src/worker.tsmodified+28−2View file
@@ -1,15 +1,17 @@
11 // Web Worker: runs the numbl script, computes linear modes, and simulates.
22
3-import { relDiff } from "./sim/analysis";
3+import { growthRate, relDiff } from "./sim/analysis";
44 import { computeModes, isPlanar } from "./sim/modes";
55 import { numblEig, runScript } from "./sim/runScript";
66 import { simulate } from "./sim/solver";
7+import { runTwin } from "./sim/twin";
78 import type { BeadSystem, SimResult } from "./sim/types";
89
910 const helperSources = import.meta.glob("../matlab/helpers/*.m", { query: "?raw", import: "default", eager: true }) as Record<string, string>;
1011 const helpers = Object.entries(helperSources).map(([path, source]) => ({ name: path.split("/").pop()!, source }));
1112
1213 const MAX_MODE_BEADS = 1500;
14+const MAX_TWIN_BYTES = 400e6;
1315
1416 export interface RunPayload {
1517 sys: {
@@ -30,6 +32,9 @@ export interface RunPayload {
3032 check: { dt: number; windows: { t: number; rel: number }[] } | null;
3133 // the same run with every hit scaled by 1e-3, its output scaled back up
3234 linear: { audio: Float32Array; windows: { t: number; rel: number }[] } | null;
35+ // the same run with every hit scaled by (1 + eps): separation of the bead
36+ // positions (max over coordinates, relative to the largest displacement)
37+ twin: { eps: number; audio: Float32Array; times: Float64Array; sep: Float64Array; rate: number | null } | null;
3338 }
3439
3540 export interface RunRequest {
@@ -37,6 +42,7 @@ export interface RunRequest {
3742 script: string;
3843 dtCheck: boolean;
3944 linear: boolean;
45+ twin: boolean;
4046 }
4147
4248 export type FromWorker =
@@ -52,7 +58,7 @@ const windowsOf = (duration: number) =>
5258 [0.01, 0.03, 0.1, 0.3, 1, duration].filter((t, i, a) => t <= duration && a.indexOf(t) === i);
5359
5460 self.onmessage = (ev: MessageEvent<RunRequest>) => {
55- const { script, dtCheck, linear } = ev.data;
61+ const { script, dtCheck, linear, twin } = ev.data;
5662 try {
5763 post({ type: "stage", text: "Running script" });
5864 const { sys, warnings } = runScript(script, helpers, text => post({ type: "log", text }));
@@ -69,10 +75,17 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
6975 }
7076 }
7177
78+ // positions at every frame, kept for the perturbed twin
79+ const saved: Float64Array[] = [];
80+ const savedBytes = Math.round(sys.duration * sys.frameRate + 1) * 3 * sys.n * 8;
81+ const keep = twin && savedBytes <= MAX_TWIN_BYTES;
82+ if (twin && !keep) notes.push(`Perturbed twin skipped: it would store ${(savedBytes / 1e6).toFixed(0)} MB of positions (limit ${MAX_TWIN_BYTES / 1e6} MB). Lower sys.frameRate or sys.duration.`);
83+
7284 post({ type: "stage", text: "Simulating" });
7385 const result = simulate(sys, {
7486 modes,
7587 onProgress: t => post({ type: "progress", t, duration: sys.duration, label: "Simulating" }),
88+ onFrame: keep ? (_f, _t, x) => saved.push(x.slice()) : undefined,
7689 });
7790
7891 let checkOut: RunPayload["check"] = null;
@@ -111,6 +124,17 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
111124 linearOut = { audio, windows };
112125 }
113126
127+ let twinOut: RunPayload["twin"] = null;
128+ if (keep && !result.aborted) {
129+ post({ type: "stage", text: "Simulating the perturbed twin" });
130+ const eps = 1e-9;
131+ const { audio, sep } = runTwin(sys, saved, result.dt, eps, t =>
132+ post({ type: "progress", t, duration: sys.duration, label: "Perturbed twin" })
133+ );
134+ const times = result.frameTimes.slice(0, saved.length);
135+ twinOut = { eps, audio, times, sep, rate: growthRate(times, sep, eps) };
136+ }
137+
114138 const hitBeads = Int32Array.from(new Set(sys.hits.flatMap(h => Array.from(h.beads))));
115139 const payload: RunPayload = {
116140 sys: {
@@ -130,10 +154,12 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
130154 notes,
131155 check: checkOut,
132156 linear: linearOut,
157+ twin: twinOut,
133158 };
134159 const transfer: Transferable[] = [result.audio.buffer, ...result.animFrames.map(f => f.buffer)];
135160 if (result.modeEnergy) transfer.push(result.modeEnergy.buffer);
136161 if (linearOut) transfer.push(linearOut.audio.buffer);
162+ if (twinOut) transfer.push(twinOut.audio.buffer, twinOut.sep.buffer);
137163 post({ type: "result", payload }, transfer);
138164 } catch (e) {
139165 const err = e as { message?: string; line?: number; file?: string };
tests/solver.test.tsmodified+24−0View file
@@ -2,7 +2,9 @@ import { readFileSync, readdirSync } from "fs";
22 import { describe, expect, it } from "vitest";
33 import { computeModes } from "../src/sim/modes";
44 import { numblEig, runScript } from "../src/sim/runScript";
5+import { growthRate } from "../src/sim/analysis";
56 import { frequencyErrorCents, simulate, stabilityLimit } from "../src/sim/solver";
7+import { runTwin } from "../src/sim/twin";
68 import type { BeadSystem } from "../src/sim/types";
79
810 const helpers = readdirSync("matlab/helpers").map(name => ({
@@ -169,3 +171,25 @@ describe("script errors", () => {
169171 expect(() => build(STRING("sys.springs = [1 2 3];"))).toThrow(/4 or 5 columns/);
170172 });
171173 });
174+
175+describe("chaos example", () => {
176+ const src = readFileSync("matlab/examples/chaos.m", "utf8");
177+ const rate = (strength: number) => {
178+ const { sys } = build(src.replace(/^s = 1;/m, `s = ${strength};`));
179+ const s1 = { ...sys, duration: 1, frameRate: 1000 };
180+ const saved: Float64Array[] = [];
181+ const r = simulate(s1, { onFrame: (_f, _t, x) => saved.push(x.slice()) });
182+ const { sep } = runTwin(s1, saved, r.dt, 1e-9);
183+ return growthRate(r.frameTimes, sep, 1e-9);
184+ };
185+
186+ it("is chaotic when struck hard", () => {
187+ const g = rate(1);
188+ expect(g).not.toBeNull();
189+ expect(g!).toBeGreaterThan(50); // e-folding time under 20 ms
190+ });
191+
192+ it("is regular when struck softly", () => {
193+ expect(rate(0.05)).toBeNull();
194+ });
195+});
tests/twin.tsadded+21−0View file
@@ -0,0 +1,21 @@
1+// Separation of two runs whose hits differ by a factor (1 + eps).
2+// Usage: tsx tests/twin.ts script.m duration [eps]
3+import { readFileSync, readdirSync } from "fs";
4+import { runScript } from "../src/sim/runScript";
5+import { simulate } from "../src/sim/solver";
6+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+sys.duration = Number(process.argv[3]);
9+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); } });
18+const out: string[] = [];
19+const step = Math.max(1, Math.floor(sep.length / 14));
20+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(" ")}`);