Dome example: chaos that keeps its pitch. Constant load with static relaxation (FIRE), full 3D linear modes for curved nets, spectrum-based default frequency range, screening script
14 changed files+436−124
README.mdmodified+7−1View file
@@ -20,6 +20,7 @@ Each bead has a mass and a position in 3D. The script defines:
2020 - *Dashpots* between pairs of beads, with force c (v_j - v_i) on bead i. Short-wavelength motion has larger relative velocities between neighbours, so this damping grows with frequency.
2121 - *Friction* on each bead, force -c0 v, which damps all frequencies alike.
2222 - *Pinned* beads, held fixed.
23+- A constant *load* on the beads (like air pressure under a drum head, or gravity). With a load, the solver first moves the beads to static equilibrium (`sys.relax`, by the FIRE minimizer), so a net can rest in a curved shape held by tension alone.
2324 - *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.
2425 - A *listener*: a bead (or the average of several) and a direction; the displacement along that direction is the audio signal.
2526
@@ -43,7 +44,7 @@ Each run reports:
4344
4445 ## Linear modes
4546
46-For a net that rests flat (all beads at one z), small transverse motion decouples from in-plane motion, and its stiffness matrix is a graph Laplacian weighted by spring tension over length. The app computes its eigenpairs with numbl's `eig` (up to 1500 free beads; about 1 s for the 721 free beads of the drum example) and records the energy of every mode over time. The "energy by linear mode" plot shows the share of the quadratic energy in log-spaced frequency bands: in a linear system each band only decays, so any band that gains energy is fed by the nonlinearity.
47+For a net that rests flat (all beads at one z), small transverse motion decouples from in-plane motion, and its stiffness matrix is a graph Laplacian weighted by spring tension over length. For any other rest shape (a dome, a 3D structure) the app uses the full 3n-by-3n stiffness of the springs. It computes the eigenpairs with numbl's `eig` (up to 1500 degrees of freedom; about 1 s for the 721 free beads of the drum example) and records the energy of every mode over time. The "energy by linear mode" plot shows the share of the quadratic energy in log-spaced frequency bands: in a linear system each band only decays, so any band that gains energy is fed by the nonlinearity.
4748
4849 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.
4950
@@ -58,6 +59,7 @@ The "perturbed twin" option runs the system again with every hit changed by one
5859 | `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. |
5960 | `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. |
6061 | `string` | 120 beads on a string, struck slightly off-axis; fast. |
62+| `dome` | Chaos that keeps its pitch. A small drum head (127 beads) bulged 4.6 mm into a dome by a steady load and struck hard: the twin separates with an e-folding time of 8 ms, while the main partial (582 Hz) moves by only about 20 cents and the spectrum stays tonal. Flat nets of the same kind become chaotic only when struck hard enough to move the pitch by about a semitone. About real time. |
6163 | `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. |
6264
6365 Timings are for a laptop (Chrome, single thread).
@@ -68,3 +70,7 @@ Timings are for a laptop (Chrome, single thread).
6870 - The solver is plain TypeScript in one Web Worker. Cost scales with the number of springs times the number of steps; nets of several thousand beads will need WASM or WebGPU.
6971 - A bead net has an upper frequency set by its spacing (about 4 kHz for the 1 cm drum example), so a dense high-frequency spectrum needs many beads.
7072 - Linear modes are computed only for flat nets.
73+
74+## What makes chaos sound good here
75+
76+Measured over several nets (2026-10-06): a net held only by tension (flat, or a string) becomes chaotic at moderate size only when the strike is strong enough to raise its pitch by 90 cents or more, and the two-bead `chaos` example, which needs a very hard strike, shifts by 3 semitones or more and sounds harsh. The reason appears to be that a flat net's nonlinearity is cubic (stretching stiffens it), so the coupling that moves energy between modes comes with a pitch shift of the same order. A curved net has a quadratic coupling as well, through which two modes drive a third directly, as in gongs and cymbals; in the `dome` example the motion is chaotic with a pitch change of about 20 cents. `tests/screen.ts` prints the measurements (twin growth rate, nonlinear pitch shift against a small-strike run, spectral flatness) for any script.
matlab/examples/dome.madded+41−0View file
@@ -0,0 +1,41 @@
1+% Chaos that keeps its pitch: a small drum head bulged into a shallow dome
2+% by a steady load (like air pressure underneath), struck hard.
3+%
4+% A flat net only becomes chaotic when struck hard enough to bend its pitch
5+% by about a semitone, and then it sounds harsh. Curvature changes that:
6+% in a dome, two modes can drive a third directly (as in gongs and
7+% cymbals), so energy moves between modes at much smaller amplitude. Here
8+% the motion is chaotic (turn on "perturbed twin") while the main partial
9+% moves by only about 20 cents; what changes from strike to strike is the
10+% shimmer of the partials, not the pitch. Try force 5 (regular) and 15
11+% (chaotic), and compare with load = 0 (flat).
12+
13+R = 0.06; % radius (m)
14+a = 0.01; % bead spacing (m)
15+sigma = 1; % areal density (kg/m^2)
16+T = 2000; % membrane tension (N/m)
17+strain = 0.005; % spring stretch at rest
18+
19+[pos, tri, rim] = bp_hexdisk(R, a);
20+E = bp_edges(tri);
21+f0 = T * a / sqrt(3);
22+k = f0 / (strain * a);
23+
24+sys.pos = pos;
25+sys.mass = sigma * (sqrt(3)/2) * a^2;
26+sys.pinned = rim;
27+sys.springs = bp_springs(pos, E, k, f0);
28+[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 2.5; 2000 1.5]);
29+
30+% Steady upward force on every bead (N). The solver first finds the domed
31+% shape where springs and load balance (sys.relax, on by default with a
32+% load); the sag is reported with the run.
33+sys.load = [0 0 2];
34+
35+hit.beads = bp_near(pos, [0.3*R 0.1*R 0], 1.5*a);
36+hit.force = [0 0 -15]; % peak force per bead (N)
37+hit.duration = 2e-3;
38+sys.hits = hit;
39+
40+sys.listen = struct('beads', bp_near(pos, [-0.3*R 0.2*R 0], 0), 'dir', [0 0 1]);
41+sys.duration = 2.5;
src/App.tsxmodified+21−4View file
@@ -1,5 +1,5 @@
11 import { useCallback, useEffect, useRef, useState } from "react";
2-import { wavBlob } from "./sim/analysis";
2+import { spectrum, wavBlob } from "./sim/analysis";
33 import { BeadView } from "./ui/BeadView";
44 import { Editor } from "./ui/Editor";
55 import { DecayTimes, Energy, ModeEnergy, Separation, Spectrogram, Spectrum, Waveform } from "./ui/Plots";
@@ -12,6 +12,23 @@ const examples = Object.entries(exampleSources)
1212 .sort((a, b) => (a.name === "drum" ? -1 : b.name === "drum" ? 1 : a.name.localeCompare(b.name)));
1313 const helpers = Object.entries(helperSources).map(([path, source]) => ({ name: path.split("/").pop()!, source }));
1414
15+// Smallest display range holding 99.9% of the power in the first half second.
16+function defaultFmax(audio: Float32Array, fs: number): number {
17+ const s = spectrum(audio, 0, Math.min(audio.length, Math.round(0.5 * fs)));
18+ const p = Array.from(s.db, d => 10 ** (d / 10));
19+ const total = p.reduce((a, b) => a + b, 0);
20+ let acc = 0;
21+ let f99 = fs / 2;
22+ for (let b = 0; b < p.length; b++) {
23+ acc += p[b];
24+ if (acc >= 0.999 * total) {
25+ f99 = b * s.df * fs;
26+ break;
27+ }
28+ }
29+ return [1000, 2000, 5000, 10000, 20000].find(f => f >= 1.2 * f99 && f <= fs / 2) ?? Math.min(20000, fs / 2);
30+}
31+
1532 const STORE_KEY = "bead-percussion:script";
1633 const load = () => {
1734 try {
@@ -83,9 +100,7 @@ export default function App() {
83100 } else if (m.type === "result") {
84101 setPayload(m.payload);
85102 setRunId(i => i + 1);
86- const fsHalf = m.payload.result.sampleRate / 2;
87- const top = m.payload.modeFreqs ? m.payload.modeFreqs[m.payload.modeFreqs.length - 1] * 1.15 : 10000;
88- setFmax([1000, 2000, 5000, 10000, 20000].find(f => f >= top && f <= fsHalf) ?? Math.min(20000, fsHalf));
103+ setFmax(defaultFmax(m.payload.result.audio, m.payload.result.sampleRate));
89104 setRunning(false);
90105 w.terminate();
91106 }
@@ -388,6 +403,8 @@ Elements
388403 [sys.friction, sys.dashpots] = bp_decay(pos, mass, sys.springs, [f1 T1; f2 T2])
389404 sets both from decay times (T60, s) at two frequencies (Hz)
390405 sys.pinned bead indices (or logical mask) held fixed
406+ sys.load constant force on every bead: 1-by-3 (same for all) or n-by-3 (N)
407+ sys.relax move beads to static equilibrium first (default: true when there is a load)
391408 Excitation
392409 sys.hits struct or struct array, each with
393410 beads, t (default 0), and either
src/sim/forces.tsadded+74−0View file
@@ -0,0 +1,74 @@
1+// Spring forces and potential energy, shared by the time stepper and the
2+// static relaxation.
3+
4+import type { BeadSystem } from "./types";
5+
6+export function springForces(sys: BeadSystem, x: Float64Array, F: Float64Array): number {
7+ const { si, sj, sk, sL0, sk3 } = sys;
8+ const ns = si.length;
9+ F.fill(0);
10+ let pe = 0;
11+ for (let s = 0; s < ns; s++) {
12+ const i3 = 3 * si[s];
13+ const j3 = 3 * sj[s];
14+ const dx = x[j3] - x[i3];
15+ const dy = x[j3 + 1] - x[i3 + 1];
16+ const dz = x[j3 + 2] - x[i3 + 2];
17+ const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
18+ const e = L - sL0[s];
19+ const k = sk[s];
20+ const k3 = sk3[s];
21+ const e2 = e * e;
22+ pe += 0.5 * k * e2 + 0.25 * k3 * e2 * e2;
23+ const g = (k * e + k3 * e2 * e) / L;
24+ F[i3] += g * dx;
25+ F[i3 + 1] += g * dy;
26+ F[i3 + 2] += g * dz;
27+ F[j3] -= g * dx;
28+ F[j3 + 1] -= g * dy;
29+ F[j3 + 2] -= g * dz;
30+ }
31+ return pe;
32+}
33+
34+// Constant load: adds it to F and returns its potential energy, -load . (x - pos0).
35+export function loadForces(sys: BeadSystem, x: Float64Array, F: Float64Array): number {
36+ const load = sys.load;
37+ if (!load) return 0;
38+ let pe = 0;
39+ for (let q = 0; q < load.length; q++) {
40+ F[q] += load[q];
41+ pe -= load[q] * (x[q] - sys.pos0[q]);
42+ }
43+ return pe;
44+}
45+
46+// Per-spring 3x3 stiffness blocks at the rest positions,
47+// B = keff d d^T + (f/L)(I - d d^T), stored as [xx, yy, zz, xy, xz, yz].
48+// The stiffness matrix is the sum over springs of B on the (i,i) and (j,j)
49+// blocks and -B on the (i,j) and (j,i) blocks.
50+export function springBlocks(sys: BeadSystem): Float64Array {
51+ const { si, sj, sk, sL0, sk3, pos0 } = sys;
52+ const ns = si.length;
53+ const B = new Float64Array(ns * 6);
54+ for (let s = 0; s < ns; s++) {
55+ const i3 = 3 * si[s];
56+ const j3 = 3 * sj[s];
57+ const dx = pos0[j3] - pos0[i3];
58+ const dy = pos0[j3 + 1] - pos0[i3 + 1];
59+ const dz = pos0[j3 + 2] - pos0[i3 + 2];
60+ const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
61+ const e = L - sL0[s];
62+ const keff = sk[s] + 3 * sk3[s] * e * e;
63+ const g = (sk[s] * e + sk3[s] * e * e * e) / L;
64+ const u = [dx / L, dy / L, dz / L];
65+ const a = keff - g;
66+ B[6 * s] = a * u[0] * u[0] + g;
67+ B[6 * s + 1] = a * u[1] * u[1] + g;
68+ B[6 * s + 2] = a * u[2] * u[2] + g;
69+ B[6 * s + 3] = a * u[0] * u[1];
70+ B[6 * s + 4] = a * u[0] * u[2];
71+ B[6 * s + 5] = a * u[1] * u[2];
72+ }
73+ return B;
74+}
src/sim/modes.tsmodified+64−43View file
@@ -1,11 +1,13 @@
1-// Linear transverse modes of a flat bead system.
1+// Linear modes about the rest positions.
22 //
3-// When every bead rests in a plane z = const and every spring lies in that
4-// plane, small transverse (z) motion decouples from in-plane motion, and its
5-// stiffness matrix is a graph Laplacian weighted by spring tension over length:
6-// K = sum_s (f_s / L_s) (e_i - e_j)(e_i - e_j)^T. We return the eigenpairs of
7-// M^{-1/2} K M^{-1/2} over the free beads.
3+// For a flat net (all beads at one z, all springs in that plane), small
4+// transverse (z) motion decouples from in-plane motion, and its stiffness is
5+// a graph Laplacian weighted by spring tension over length,
6+// K = sum_s (f_s / L_s) (e_i - e_j)(e_i - e_j)^T; we use that n-by-n problem.
7+// Otherwise we use the full 3n-by-3n stiffness of the springs. Either way we
8+// return the eigenpairs of M^{-1/2} K M^{-1/2} over the free degrees of freedom.
89
10+import { springBlocks } from "./forces";
911 import type { BeadSystem, Modes } from "./types";
1012
1113 export type SymEig = (A: Float64Array, n: number) => { values: Float64Array; vectors: Float64Array };
@@ -23,65 +25,84 @@ export function isPlanar(sys: BeadSystem): boolean {
2325 return zmax - zmin <= 1e-12 * Math.max(span, 1e-30);
2426 }
2527
26-export function transverseStiffness(sys: BeadSystem): { K: Float64Array; free: Int32Array } {
27- const { n, pinned, si, sj, sk, sL0, sk3, pos0, mass } = sys;
28- const freeIdx = new Int32Array(n).fill(-1);
29- const free: number[] = [];
30- for (let i = 0; i < n; i++) if (!pinned[i]) freeIdx[i] = free.push(i) - 1;
31- const nf = free.length;
32- const K = new Float64Array(nf * nf);
28+// Number of degrees of freedom computeModes would use.
29+export function modeSize(sys: BeadSystem): number {
30+ const nFree = sys.pinned.reduce((a, p) => a + (p ? 0 : 1), 0);
31+ return isPlanar(sys) ? nFree : 3 * nFree;
32+}
33+
34+function stiffness(sys: BeadSystem): { K: Float64Array; dof: Int32Array } {
35+ const { n, pinned, si, sj, mass } = sys;
36+ const planar = isPlanar(sys);
37+ const comps = planar ? [2] : [0, 1, 2];
38+ const row = new Int32Array(3 * n).fill(-1);
39+ const dof: number[] = [];
40+ for (let i = 0; i < n; i++) if (!pinned[i]) for (const c of comps) row[3 * i + c] = dof.push(3 * i + c) - 1;
41+ const nd = dof.length;
42+ const K = new Float64Array(nd * nd);
43+ const B = springBlocks(sys);
44+ const idx = [
45+ [0, 3, 4],
46+ [3, 1, 5],
47+ [4, 5, 2],
48+ ];
49+ const add = (r: number, c: number, v: number) => {
50+ if (r >= 0 && c >= 0) K[r * nd + c] += v;
51+ };
3352 for (let s = 0; s < si.length; s++) {
3453 const i3 = 3 * si[s];
3554 const j3 = 3 * sj[s];
36- const L = Math.hypot(pos0[j3] - pos0[i3], pos0[j3 + 1] - pos0[i3 + 1], pos0[j3 + 2] - pos0[i3 + 2]);
37- const e = L - sL0[s];
38- const g = (sk[s] * e + sk3[s] * e * e * e) / L;
39- const a = freeIdx[si[s]];
40- const b = freeIdx[sj[s]];
41- if (a >= 0) K[a * nf + a] += g;
42- if (b >= 0) K[b * nf + b] += g;
43- if (a >= 0 && b >= 0) {
44- K[a * nf + b] -= g;
45- K[b * nf + a] -= g;
55+ for (const a of comps) {
56+ for (const b of comps) {
57+ const v = B[6 * s + idx[a][b]];
58+ add(row[i3 + a], row[i3 + b], v);
59+ add(row[j3 + a], row[j3 + b], v);
60+ add(row[i3 + a], row[j3 + b], -v);
61+ add(row[j3 + a], row[i3 + b], -v);
62+ }
4663 }
4764 }
48- for (let a = 0; a < nf; a++) {
49- for (let b = 0; b < nf; b++) K[a * nf + b] /= Math.sqrt(mass[free[a]] * mass[free[b]]);
65+ for (let r = 0; r < nd; r++) {
66+ for (let c = 0; c < nd; c++) K[r * nd + c] /= Math.sqrt(mass[Math.floor(dof[r] / 3)] * mass[Math.floor(dof[c] / 3)]);
5067 }
51- return { K, free: Int32Array.from(free) };
68+ return { K, dof: Int32Array.from(dof) };
5269 }
5370
5471 export function computeModes(sys: BeadSystem, eig: SymEig): Modes {
55- const { K, free } = transverseStiffness(sys);
56- const nf = free.length;
57- const { values, vectors } = eig(K, nf);
72+ const { K, dof } = stiffness(sys);
73+ const nd = dof.length;
74+ const { values, vectors } = eig(K, nd);
5875 const order = Array.from(values.keys()).sort((p, q) => values[p] - values[q]);
59- const freqs = new Float64Array(nf);
60- const psi = new Float64Array(nf * nf);
76+ const freqs = new Float64Array(nd);
77+ const psi = new Float64Array(nd * nd);
6178 order.forEach((src, k) => {
6279 freqs[k] = Math.sqrt(Math.max(values[src], 0)) / (2 * Math.PI);
63- psi.set(vectors.subarray(src * nf, (src + 1) * nf), k * nf);
80+ psi.set(vectors.subarray(src * nd, (src + 1) * nd), k * nd);
6481 });
65- return { freqs, psi, free };
82+ return { freqs, psi, dof, kind: isPlanar(sys) ? "transverse" : "full" };
6683 }
6784
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
85+// Small-amplitude decay rate of each mode (amplitude, 1/s), from the friction
86+// and dashpots: sigma_k = (1/2) u^T C u with u the mode shape normalized so
87+// u^T M u = 1. Exact when the damping matrix is a combination of M and K (as
88+// bp_decay builds it for flat nets); otherwise it neglects the coupling of
7289 // modes by damping, which is small for light damping.
7390 export function modalDamping(sys: BeadSystem, modes: Modes): Float64Array {
7491 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);
92+ const { dof, psi, freqs } = modes;
93+ const nd = dof.length;
94+ const u = new Float64Array(3 * n);
7895 const out = new Float64Array(freqs.length);
7996 for (let k = 0; k < freqs.length; k++) {
8097 u.fill(0);
81- for (let r = 0; r < nf; r++) u[free[r]] = psi[k * nf + r] / Math.sqrt(mass[free[r]]);
98+ for (let r = 0; r < nd; r++) u[dof[r]] = psi[k * nd + r] / Math.sqrt(mass[Math.floor(dof[r] / 3)]);
8299 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;
100+ for (let i = 0; i < n; i++) c += friction[i] * (u[3 * i] ** 2 + u[3 * i + 1] ** 2 + u[3 * i + 2] ** 2);
101+ for (let d = 0; d < di.length; d++) {
102+ const i3 = 3 * di[d];
103+ const j3 = 3 * dj[d];
104+ c += dc[d] * ((u[i3] - u[j3]) ** 2 + (u[i3 + 1] - u[j3 + 1]) ** 2 + (u[i3 + 2] - u[j3 + 2]) ** 2);
105+ }
85106 out[k] = c / 2;
86107 }
87108 return out;
src/sim/parse.tsmodified+15−2View file
@@ -87,7 +87,7 @@ function fields(v: RuntimeValue | undefined, name: string): Map<string, RuntimeV
8787
8888 const KNOWN = new Set([
8989 "pos", "mass", "pinned", "springs", "dashpots", "friction", "hits", "listen",
90- "duration", "sampleRate", "dt", "safety", "method", "frameRate", "modes",
90+ "duration", "sampleRate", "dt", "safety", "method", "frameRate", "modes", "load", "relax",
9191 ]);
9292
9393 export function parseSystem(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warnings: string[] } {
@@ -170,6 +170,18 @@ function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warn
170170 warnings.push(`No damping given: using the default, every mode falls 60 dB in ${DEFAULT_T60} s. Set sys.friction = 0 for none.`);
171171 }
172172
173+ let load: Float64Array | null = null;
174+ const Lm = asMat(f.get("load"), "sys.load");
175+ if (Lm && Lm.data.length) {
176+ load = new Float64Array(3 * n);
177+ if (Lm.data.length === 3) {
178+ for (let i = 0; i < n; i++) for (let c = 0; c < 3; c++) load[3 * i + c] = Lm.data[c];
179+ } else if (Lm.rows === n && Lm.cols === 3) {
180+ for (let i = 0; i < n; i++) for (let c = 0; c < 3; c++) load[3 * i + c] = Lm.data[i + c * n];
181+ } else throw new SysError("sys.load must be 1-by-3 (the same force on every bead) or n-by-3");
182+ for (let i = 0; i < n; i++) if (pinned[i]) load.fill(0, 3 * i, 3 * i + 3);
183+ }
184+
173185 const hits: Hit[] = fields(f.get("hits"), "sys.hits").map((h, q) => {
174186 const nm = `sys.hits(${q + 1})`;
175187 const beads = indices(h.get("beads"), `${nm}.beads`, n);
@@ -202,7 +214,8 @@ function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warn
202214 if (method !== "verlet" && method !== "yoshida4") throw new SysError("sys.method must be 'verlet' or 'yoshida4'");
203215
204216 const sys: BeadSystem = {
205- n, pos0, mass, pinned, si, sj, sk, sL0, sk3, di, dj, dc, friction, hits, listenBeads,
217+ n, pos0, mass, pinned, si, sj, sk, sL0, sk3, di, dj, dc, friction, load, hits, listenBeads,
218+ relax: scalar(f.get("relax"), "sys.relax", load ? 1 : 0) !== 0,
206219 listenDir: [dir[0] / dn, dir[1] / dn, dir[2] / dn],
207220 duration: scalar(f.get("duration"), "sys.duration", 2),
208221 sampleRate: scalar(f.get("sampleRate"), "sys.sampleRate", 48000),
src/sim/relax.tsadded+83−0View file
@@ -0,0 +1,83 @@
1+// Static equilibrium of the springs under the constant load, by FIRE
2+// (Bitzek et al. 2006): damped dynamics that keeps the velocity aligned with
3+// the force and lengthens its step while the energy keeps falling.
4+
5+import { loadForces, springForces } from "./forces";
6+import { estimateOmegaMax } from "./solver";
7+import type { BeadSystem } from "./types";
8+
9+export function relax(sys: BeadSystem, maxIter = 200000): { pos: Float64Array; iterations: number; residual: number; maxTension: number } {
10+ const { n, pinned, mass } = sys;
11+ const N3 = 3 * n;
12+ const x = new Float64Array(sys.pos0);
13+ const v = new Float64Array(N3);
14+ const F = new Float64Array(N3);
15+ const invm = new Float64Array(n);
16+ for (let i = 0; i < n; i++) invm[i] = pinned[i] ? 0 : 1 / mass[i];
17+ const force = () => {
18+ springForces(sys, x, F);
19+ loadForces(sys, x, F);
20+ for (let i = 0; i < n; i++) if (pinned[i]) F.fill(0, 3 * i, 3 * i + 3);
21+ };
22+ const residual = () => {
23+ let r = 0;
24+ for (let i = 0; i < n; i++) r = Math.max(r, Math.hypot(F[3 * i], F[3 * i + 1], F[3 * i + 2]));
25+ return r;
26+ };
27+ let scale = 0; // force scale for the tolerance: largest spring or load force
28+ for (let s = 0; s < sys.si.length; s++) scale = Math.max(scale, sys.sk[s] * Math.abs(sys.sL0[s]) * 1e-3);
29+ if (sys.load) for (const l of sys.load) scale = Math.max(scale, Math.abs(l));
30+ const tol = 1e-10 * Math.max(scale, 1e-30);
31+
32+ const dt0 = 0.2 / estimateOmegaMax(sys);
33+ let dt = dt0;
34+ const dtMax = 10 * dt0;
35+ let alpha = 0.1;
36+ let positive = 0;
37+ force();
38+ let it = 0;
39+ for (; it < maxIter; it++) {
40+ if (residual() < tol) break;
41+ let P = 0, vn = 0, fn = 0;
42+ for (let q = 0; q < N3; q++) {
43+ P += F[q] * v[q];
44+ vn += v[q] * v[q];
45+ fn += F[q] * F[q];
46+ }
47+ vn = Math.sqrt(vn);
48+ fn = Math.sqrt(fn);
49+ if (P > 0) {
50+ for (let q = 0; q < N3; q++) v[q] = (1 - alpha) * v[q] + (alpha * vn * F[q]) / (fn || 1);
51+ if (++positive > 5) {
52+ dt = Math.min(dt * 1.1, dtMax);
53+ alpha *= 0.99;
54+ }
55+ } else {
56+ v.fill(0);
57+ dt *= 0.5;
58+ alpha = 0.1;
59+ positive = 0;
60+ }
61+ for (let i = 0; i < n; i++) {
62+ const h = 0.5 * dt * invm[i];
63+ v[3 * i] += h * F[3 * i];
64+ v[3 * i + 1] += h * F[3 * i + 1];
65+ v[3 * i + 2] += h * F[3 * i + 2];
66+ }
67+ for (let q = 0; q < N3; q++) x[q] += dt * v[q];
68+ force();
69+ for (let i = 0; i < n; i++) {
70+ const h = 0.5 * dt * invm[i];
71+ v[3 * i] += h * F[3 * i];
72+ v[3 * i + 1] += h * F[3 * i + 1];
73+ v[3 * i + 2] += h * F[3 * i + 2];
74+ }
75+ }
76+ let maxTension = 0;
77+ for (let s = 0; s < sys.si.length; s++) {
78+ const i3 = 3 * sys.si[s], j3 = 3 * sys.sj[s];
79+ const L = Math.hypot(x[j3] - x[i3], x[j3 + 1] - x[i3 + 1], x[j3 + 2] - x[i3 + 2]);
80+ maxTension = Math.max(maxTension, Math.abs(sys.sk[s] * (L - sys.sL0[s])));
81+ }
82+ return { pos: x, iterations: it, residual: residual(), maxTension };
83+}
src/sim/runScript.tsmodified+17−0View file
@@ -3,6 +3,7 @@
33
44 import { executeCode, RTV } from "numbl";
55 import { parseSystem } from "./parse";
6+import { relax } from "./relax";
67 import type { SymEig } from "./modes";
78 import type { BeadSystem } from "./types";
89
@@ -20,6 +21,22 @@ export function runScript(
2021 return parseSystem(result.variableValues.sys);
2122 }
2223
24+// Steps that change the system before it is simulated (now: static relaxation
25+// under the load). Returns notes for the run report.
26+export function prepareSystem(sys: BeadSystem): string[] {
27+ const notes: string[] = [];
28+ if (sys.relax) {
29+ const r = relax(sys);
30+ let sag = 0;
31+ for (let q = 0; q < r.pos.length; q++) sag = Math.max(sag, Math.abs(r.pos[q] - sys.pos0[q]));
32+ sys.pos0 = r.pos;
33+ notes.push(
34+ `Relaxed to equilibrium in ${r.iterations} iterations (largest net force ${r.residual.toExponential(1)} N; beads moved up to ${(sag * 1000).toPrecision(3)} mm).`
35+ );
36+ }
37+ return notes;
38+}
39+
2340 // Symmetric eigendecomposition through numbl's eig.
2441 export const numblEig: SymEig = (A, n) => {
2542 const r = executeCode("[V, D] = eig(A); d = diag(D);", {
src/sim/solver.tsmodified+17−59View file
@@ -9,6 +9,7 @@
99 // which is small compared with the oscillation frequencies; the conservative
1010 // part is second order (Verlet) or fourth order (Yoshida) in the step h.
1111
12+import { loadForces, springBlocks, springForces } from "./forces";
1213 import type { BeadSystem, Method, Modes, SimResult } from "./types";
1314
1415 const YOSHIDA_W1 = 1 / (2 - Math.cbrt(2));
@@ -86,31 +87,8 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
8687 // ---- forces ------------------------------------------------------------
8788 const { si, sj, sk, sL0, sk3 } = sys;
8889 const ns = si.length;
89- const internalForces = (): number => {
90- F.fill(0);
91- let pe = 0;
92- for (let s = 0; s < ns; s++) {
93- const i3 = 3 * si[s];
94- const j3 = 3 * sj[s];
95- const dx = x[j3] - x[i3];
96- const dy = x[j3 + 1] - x[i3 + 1];
97- const dz = x[j3 + 2] - x[i3 + 2];
98- const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
99- const e = L - sL0[s];
100- const k = sk[s];
101- const k3 = sk3[s];
102- const e2 = e * e;
103- pe += 0.5 * k * e2 + 0.25 * k3 * e2 * e2;
104- const g = (k * e + k3 * e2 * e) / L;
105- F[i3] += g * dx;
106- F[i3 + 1] += g * dy;
107- F[i3 + 2] += g * dz;
108- F[j3] -= g * dx;
109- F[j3 + 1] -= g * dy;
110- F[j3 + 2] -= g * dz;
111- }
112- return pe;
113- };
90+ // spring forces plus the constant load; returns their potential energy
91+ const internalForces = (): number => springForces(sys, x, F) + loadForces(sys, x, F);
11492 const pulses = sys.hits.filter(h => h.duration > 0);
11593 const externalForces = (t: number) => {
11694 for (const p of pulses) {
@@ -145,7 +123,7 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
145123 if (residualForce > 1e-6 * Math.max(maxTension, 1e-30)) {
146124 warnings.push(
147125 `The rest positions are not an equilibrium: the largest net force on a free bead is ${residualForce.toExponential(2)} N ` +
148- `(largest spring force ${maxTension.toExponential(2)} N). The beads will move before any strike.`
126+ `(largest spring force ${maxTension.toExponential(2)} N). The beads will move before any strike. Set sys.relax = true to find the equilibrium first.`
149127 );
150128 }
151129
@@ -296,27 +274,27 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
296274 let bv: Float64Array | null = null;
297275 let om2: Float64Array | null = null;
298276 if (modes) {
299- bz = new Float64Array(modes.free.length);
300- bv = new Float64Array(modes.free.length);
277+ bz = new Float64Array(modes.dof.length);
278+ bv = new Float64Array(modes.dof.length);
301279 om2 = new Float64Array(nModes);
302280 for (let k = 0; k < nModes; k++) om2[k] = (2 * Math.PI * modes.freqs[k]) ** 2;
303281 }
304282 const recordModes = (f: number) => {
305283 if (!modes || !modeEnergy || !bz || !bv || !om2) return;
306- const free = modes.free;
307- const nf = free.length;
308- for (let r = 0; r < nf; r++) {
309- const b = free[r];
310- const sm = Math.sqrt(mass[b]);
311- bz[r] = sm * (x[3 * b + 2] - pos0[3 * b + 2]);
312- bv[r] = sm * v[3 * b + 2];
284+ const dof = modes.dof;
285+ const nr = dof.length;
286+ for (let r = 0; r < nr; r++) {
287+ const q = dof[r];
288+ const sm = Math.sqrt(mass[Math.floor(q / 3)]);
289+ bz[r] = sm * (x[q] - pos0[q]);
290+ bv[r] = sm * v[q];
313291 }
314292 const psi = modes.psi;
315293 for (let k = 0; k < nModes; k++) {
316294 let q = 0;
317295 let qd = 0;
318- const off = k * nf;
319- for (let r = 0; r < nf; r++) {
296+ const off = k * nr;
297+ for (let r = 0; r < nr; r++) {
320298 q += psi[off + r] * bz[r];
321299 qd += psi[off + r] * bv[r];
322300 }
@@ -461,29 +439,9 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
461439 // Highest frequency (rad/s) of the system linearized about the rest positions,
462440 // by power iteration on M^{-1/2} H M^{-1/2}, with a small margin.
463441 export function estimateOmegaMax(sys: BeadSystem): number {
464- const { n, si, sj, sk, sL0, sk3, pos0, mass, pinned } = sys;
442+ const { n, si, sj, mass, pinned } = sys;
465443 const ns = si.length;
466- // per-spring 3x3 blocks B = keff d d^T + (f/L)(I - d d^T)
467- const B = new Float64Array(ns * 6);
468- for (let s = 0; s < ns; s++) {
469- const i3 = 3 * si[s];
470- const j3 = 3 * sj[s];
471- const dx = pos0[j3] - pos0[i3];
472- const dy = pos0[j3 + 1] - pos0[i3 + 1];
473- const dz = pos0[j3 + 2] - pos0[i3 + 2];
474- const L = Math.sqrt(dx * dx + dy * dy + dz * dz);
475- const e = L - sL0[s];
476- const keff = sk[s] + 3 * sk3[s] * e * e;
477- const g = (sk[s] * e + sk3[s] * e * e * e) / L;
478- const u = [dx / L, dy / L, dz / L];
479- const a = keff - g;
480- B[6 * s] = a * u[0] * u[0] + g;
481- B[6 * s + 1] = a * u[1] * u[1] + g;
482- B[6 * s + 2] = a * u[2] * u[2] + g;
483- B[6 * s + 3] = a * u[0] * u[1];
484- B[6 * s + 4] = a * u[0] * u[2];
485- B[6 * s + 5] = a * u[1] * u[2];
486- }
444+ const B = springBlocks(sys);
487445 const ism = new Float64Array(n);
488446 for (let i = 0; i < n; i++) ism[i] = pinned[i] ? 0 : 1 / Math.sqrt(mass[i]);
489447 let w = new Float64Array(3 * n);
src/sim/types.tsmodified+6−3View file
@@ -29,6 +29,8 @@ export interface BeadSystem {
2929 dj: Int32Array;
3030 dc: Float64Array;
3131 friction: Float64Array; // n, F_i = -c0_i v_i (N s/m)
32+ load: Float64Array | null; // 3n, constant force on each bead (N), e.g. air pressure or gravity
33+ relax: boolean; // move the beads to static equilibrium (springs + load) before simulating
3234 hits: Hit[];
3335 listenBeads: Int32Array;
3436 listenDir: [number, number, number]; // unit vector
@@ -48,10 +50,11 @@ export interface SimProgress {
4850
4951 export interface Modes {
5052 freqs: Float64Array; // Hz, ascending
51- // psi: orthonormal eigenvectors of M^{-1/2} K M^{-1/2} over free beads,
52- // column-major (nFree x nModes)
53+ // psi: orthonormal eigenvectors of M^{-1/2} K M^{-1/2} over the degrees of
54+ // freedom in dof, column-major (dof.length x nModes)
5355 psi: Float64Array;
54- free: Int32Array; // bead index of each row of psi
56+ dof: Int32Array; // index into the 3n coordinates (3*bead + component) of each row of psi
57+ kind: "transverse" | "full"; // z motion of a flat net only, or all motion
5558 }
5659
5760 export interface SimResult {
src/ui/Plots.tsxmodified+3−3View file
@@ -184,12 +184,12 @@ export function ModeEnergy({ p }: { p: RunPayload }) {
184184 },
185185 [bands, tmax]
186186 );
187- if (!bands) return <div className="note">Mode energies are shown when linear modes are computed (flat systems only).</div>;
187+ if (!bands) return <div className="note">Mode energies are shown when linear modes are computed.</div>;
188188 return (
189189 <>
190190 <canvas ref={ref} className="plot" />
191191 <div className="legend">
192- Share of the linear (quadratic) energy in each band of transverse modes, 8 decades of color. Show the first{" "}
192+ Share of the linear (quadratic) energy in each band of {p.modeKind === "full" ? "modes (3D)" : "transverse modes"}, 8 decades of color. Show the first{" "}
193193 <select value={tmax} onChange={e => setTmax(+e.target.value)}>
194194 {[0.03, 0.1, 0.3, 1]
195195 .filter(t => t < p.sys.duration)
@@ -347,7 +347,7 @@ export function DecayTimes({ p }: { p: RunPayload }) {
347347 <>
348348 <canvas ref={ref} className="plot short" />
349349 <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
350+ Time for each {p.modeKind === "full" ? "mode (3D)" : "transverse mode"} to fall 60 dB at small amplitude, from sys.friction and sys.dashpots. Use bp_decay to set it from two
351351 (frequency, T60) pairs.
352352 </div>
353353 </>
src/worker.tsmodified+12−8View file
@@ -1,8 +1,8 @@
11 // Web Worker: runs the numbl script, computes linear modes, and simulates.
22
33 import { growthRate, relDiff } from "./sim/analysis";
4-import { computeModes, isPlanar, modalDamping } from "./sim/modes";
5-import { numblEig, runScript } from "./sim/runScript";
4+import { computeModes, isPlanar, modalDamping, modeSize } from "./sim/modes";
5+import { numblEig, prepareSystem, runScript } from "./sim/runScript";
66 import { simulate } from "./sim/solver";
77 import { runTwin } from "./sim/twin";
88 import type { BeadSystem, SimResult } from "./sim/types";
@@ -10,7 +10,7 @@ import type { BeadSystem, SimResult } from "./sim/types";
1010 const helperSources = import.meta.glob("../matlab/helpers/*.m", { query: "?raw", import: "default", eager: true }) as Record<string, string>;
1111 const helpers = Object.entries(helperSources).map(([path, source]) => ({ name: path.split("/").pop()!, source }));
1212
13-const MAX_MODE_BEADS = 1500;
13+const MAX_MODE_DOF = 1500;
1414 const MAX_TWIN_BYTES = 400e6;
1515
1616 export interface RunPayload {
@@ -28,7 +28,8 @@ 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)
31+ modeDecay: Float64Array | null;
32+ modeKind: "transverse" | "full" | null; // small-amplitude amplitude decay rate of each mode (1/s)
3233 notes: string[];
3334 check: { dt: number; windows: { t: number; rel: number }[] } | null;
3435 // the same run with every hit scaled by 1e-3, its output scaled back up
@@ -64,14 +65,16 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
6465 post({ type: "stage", text: "Running script" });
6566 const { sys, warnings } = runScript(script, helpers, text => post({ type: "log", text }));
6667 const notes = [...warnings];
68+ if (sys.relax) post({ type: "stage", text: "Relaxing to equilibrium" });
69+ notes.push(...prepareSystem(sys));
6770
6871 let modes = null;
69- const nFree = sys.pinned.reduce((a, p) => a + (p ? 0 : 1), 0);
7072 if (sys.modes) {
71- if (!isPlanar(sys)) notes.push("Linear modes skipped: the beads do not rest in a plane z = const.");
72- else if (nFree > MAX_MODE_BEADS) notes.push(`Linear modes skipped: ${nFree} free beads (limit ${MAX_MODE_BEADS}).`);
73+ const size = modeSize(sys);
74+ const what = isPlanar(sys) ? "transverse modes of a flat net" : "3D modes";
75+ if (size > MAX_MODE_DOF) notes.push(`Linear modes skipped: ${size} degrees of freedom (${what}; limit ${MAX_MODE_DOF}).`);
7376 else {
74- post({ type: "stage", text: `Computing linear modes (${nFree} free beads)` });
77+ post({ type: "stage", text: `Computing linear modes (${size} degrees of freedom)` });
7578 modes = computeModes(sys, numblEig);
7679 }
7780 }
@@ -153,6 +156,7 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
153156 result,
154157 modeFreqs: modes ? modes.freqs : null,
155158 modeDecay: modes ? modalDamping(sys, modes) : null,
159+ modeKind: modes ? modes.kind : null,
156160 notes,
157161 check: checkOut,
158162 linear: linearOut,
tests/screen.tsadded+34−0View file
@@ -0,0 +1,34 @@
1+// Screen a script for musical chaos: twin growth rate, pitch shift of the
2+// strongest low peak vs the lowest linear mode, and spectral flatness.
3+// Usage: tsx tests/screen.ts script.m duration
4+import { readFileSync, readdirSync } from "fs";
5+import { prepareSystem, runScript } from "../src/sim/runScript";
6+import { simulate } from "../src/sim/solver";
7+import { runTwin } from "../src/sim/twin";
8+import { growthRate, spectrum } from "../src/sim/analysis";
9+const helpers = readdirSync("matlab/helpers").map(name => ({ name, source: readFileSync("matlab/helpers/" + name, "utf8") }));
10+const { sys } = runScript(readFileSync(process.argv[2], "utf8"), helpers, () => {});
11+sys.duration = Number(process.argv[3]);
12+const notes = prepareSystem(sys);
13+if (notes.length) console.log(notes.join(" "));
14+const saved: Float64Array[] = [];
15+const a = simulate(sys, { onFrame: (_f, _t, x) => saved.push(x.slice()) });
16+const { sep } = runTwin(sys, saved, a.dt, 1e-9);
17+const g = growthRate(a.frameTimes, sep, 1e-9);
18+const fs = a.sampleRate;
19+const sp = spectrum(a.audio, Math.round(0.2 * fs), Math.round(0.5 * fs));
20+const hz = (b: number) => b * sp.df * fs;
21+// flatness over 50 Hz..5 kHz
22+let lg = 0, ar = 0, cnt = 0;
23+for (let b = 0; b < sp.db.length; b++) { if (hz(b) < 50 || hz(b) > 5000) continue; const p = 10 ** (sp.db[b] / 10); lg += Math.log(p); ar += p; cnt++; }
24+const flatness = Math.exp(lg / cnt) / (ar / cnt);
25+// nonlinear pitch shift: strongest peak of the same strike scaled down by 1e3, vs the nearest peak of the actual run
26+const small = simulate({ ...sys, hits: sys.hits.map(h => ({ ...h, velocity: h.velocity.map(v => v * 1e-3) as any, force: h.force.map(v => v * 1e-3) as any })) }, { maxAnimBytes: 0 });
27+const spS = spectrum(small.audio, Math.round(0.2 * fs), Math.round(0.5 * fs));
28+let bestS = -Infinity, fS = 0;
29+for (let b = 1; b < spS.db.length; b++) if (hz(b) > 50 && spS.db[b] > bestS) (bestS = spS.db[b]), (fS = hz(b));
30+let best = -Infinity, fb = 0;
31+for (let b = 1; b < sp.db.length; b++) if (hz(b) > fS * 0.89 && hz(b) < fS * 1.12 && sp.db[b] > best) (best = sp.db[b]), (fb = hz(b));
32+const cents = 1200 * Math.log2(fb / fS);
33+const tmax = sep.findIndex(s => s > 0.1);
34+console.log(`main peak ${fS.toFixed(0)} Hz, nonlinear shift ${cents.toFixed(0)} cents | chaos rate ${g?.toFixed(0) ?? "-"} /s, sep>0.1 at ${tmax >= 0 ? a.frameTimes[tmax].toFixed(2) + " s" : "never"} | flatness ${flatness.toExponential(1)} | drift ${a.energyDrift.toExponential(1)} | ${(a.elapsedMs / 1000).toFixed(1)} s`);
tests/solver.test.tsmodified+42−1View file
@@ -1,7 +1,7 @@
11 import { readFileSync, readdirSync } from "fs";
22 import { describe, expect, it } from "vitest";
33 import { computeModes, modalDamping } from "../src/sim/modes";
4-import { numblEig, runScript } from "../src/sim/runScript";
4+import { numblEig, prepareSystem, runScript } from "../src/sim/runScript";
55 import { growthRate } from "../src/sim/analysis";
66 import { frequencyErrorCents, simulate, stabilityLimit } from "../src/sim/solver";
77 import { runTwin } from "../src/sim/twin";
@@ -258,3 +258,44 @@ describe("default damping", () => {
258258 expect(Array.from(sys.friction).every(c => c === 0)).toBe(true);
259259 });
260260 });
261+
262+describe("load and relaxation", () => {
263+ it("finds the analytic equilibrium of one bead under a sideways load", () => {
264+ // one bead between two pinned beads at distance a, springs k with rest
265+ // length L0, load P along z: 2 k (L - L0) w / L = P with L = sqrt(a^2 + w^2)
266+ const { sys } = build(STRING("sys.load = [0 0 3]; sys.friction = 0; sys.hits.beads = 2; sys.listen.beads = 2;").replace("N = 15", "N = 1"));
267+ prepareSystem(sys);
268+ const w = sys.pos0[3 * 1 + 2];
269+ const a = 0.01, k = 1e4, L0 = a - 10 / k;
270+ const L = Math.hypot(a, w);
271+ expect(Math.abs((2 * k * (L - L0) * w) / L - 3)).toBeLessThan(1e-8);
272+ const r = simulate({ ...sys, duration: 0.001 });
273+ expect(r.residualForce).toBeLessThan(1e-8);
274+ });
275+
276+ it("gives 3D modes that contain the transverse ones of a nearly flat net", () => {
277+ const { sys } = build(STRING(""));
278+ const flat = computeModes(sys, numblEig);
279+ const lifted = { ...sys, pos0: Float64Array.from(sys.pos0) };
280+ lifted.pos0[2] = 1e-12; // a pinned end lifted by 1e-12 m: no longer planar
281+ const full = computeModes(lifted, numblEig);
282+ expect(full.kind).toBe("full");
283+ expect(full.freqs.length).toBe(3 * flat.freqs.length);
284+ for (const f of flat.freqs) {
285+ const nearest = Math.min(...Array.from(full.freqs, g => Math.abs(g - f)));
286+ expect(nearest / f).toBeLessThan(1e-6);
287+ }
288+ });
289+});
290+
291+describe("dome example", () => {
292+ it("is chaotic with a small pitch change", () => {
293+ const { sys } = build(readFileSync("matlab/examples/dome.m", "utf8"));
294+ prepareSystem(sys);
295+ const s1 = { ...sys, duration: 0.6 };
296+ const saved: Float64Array[] = [];
297+ const r = simulate(s1, { onFrame: (_f, _t, x) => saved.push(x.slice()) });
298+ const { sep } = runTwin(s1, saved, r.dt, 1e-9);
299+ expect(Math.max(...sep)).toBeGreaterThan(0.1);
300+ });
301+});