Revert "Dome example: chaos that keeps its pitch..."
This reverts commit a8b5c65.
14 changed files+124−436
README.mdmodified+1−7View file
@@ -20,7 +20,6 @@ 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.
2423 - *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.
2524 - A *listener*: a bead (or the average of several) and a direction; the displacement along that direction is the audio signal.
2625
@@ -44,7 +43,7 @@ Each run reports:
4443
4544 ## Linear modes
4645
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.
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.
4847
4948 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.
5049
@@ -59,7 +58,6 @@ The "perturbed twin" option runs the system again with every hit changed by one
5958 | `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. |
6059 | `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. |
6160 | `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. |
6361 | `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. |
6462
6563 Timings are for a laptop (Chrome, single thread).
@@ -70,7 +68,3 @@ Timings are for a laptop (Chrome, single thread).
7068 - 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.
7169 - 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.
7270 - 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.mdeleted+0−41View file
@@ -1,41 +0,0 @@
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+4−21View file
@@ -1,5 +1,5 @@
11 import { useCallback, useEffect, useRef, useState } from "react";
2-import { spectrum, wavBlob } from "./sim/analysis";
2+import { 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,23 +12,6 @@ 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-
3215 const STORE_KEY = "bead-percussion:script";
3316 const load = () => {
3417 try {
@@ -100,7 +83,9 @@ export default function App() {
10083 } else if (m.type === "result") {
10184 setPayload(m.payload);
10285 setRunId(i => i + 1);
103- setFmax(defaultFmax(m.payload.result.audio, m.payload.result.sampleRate));
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));
10489 setRunning(false);
10590 w.terminate();
10691 }
@@ -403,8 +388,6 @@ Elements
403388 [sys.friction, sys.dashpots] = bp_decay(pos, mass, sys.springs, [f1 T1; f2 T2])
404389 sets both from decay times (T60, s) at two frequencies (Hz)
405390 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)
408391 Excitation
409392 sys.hits struct or struct array, each with
410393 beads, t (default 0), and either
src/sim/forces.tsdeleted+0−74View file
@@ -1,74 +0,0 @@
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+43−64View file
@@ -1,13 +1,11 @@
1-// Linear modes about the rest positions.
1+// Linear transverse modes of a flat bead system.
22 //
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.
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.
98
10-import { springBlocks } from "./forces";
119 import type { BeadSystem, Modes } from "./types";
1210
1311 export type SymEig = (A: Float64Array, n: number) => { values: Float64Array; vectors: Float64Array };
@@ -25,84 +23,65 @@ export function isPlanar(sys: BeadSystem): boolean {
2523 return zmax - zmin <= 1e-12 * Math.max(span, 1e-30);
2624 }
2725
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- };
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);
5233 for (let s = 0; s < si.length; s++) {
5334 const i3 = 3 * si[s];
5435 const j3 = 3 * sj[s];
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- }
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;
6346 }
6447 }
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)]);
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]]);
6750 }
68- return { K, dof: Int32Array.from(dof) };
51+ return { K, free: Int32Array.from(free) };
6952 }
7053
7154 export function computeModes(sys: BeadSystem, eig: SymEig): Modes {
72- const { K, dof } = stiffness(sys);
73- const nd = dof.length;
74- const { values, vectors } = eig(K, nd);
55+ const { K, free } = transverseStiffness(sys);
56+ const nf = free.length;
57+ const { values, vectors } = eig(K, nf);
7558 const order = Array.from(values.keys()).sort((p, q) => values[p] - values[q]);
76- const freqs = new Float64Array(nd);
77- const psi = new Float64Array(nd * nd);
59+ const freqs = new Float64Array(nf);
60+ const psi = new Float64Array(nf * nf);
7861 order.forEach((src, k) => {
7962 freqs[k] = Math.sqrt(Math.max(values[src], 0)) / (2 * Math.PI);
80- psi.set(vectors.subarray(src * nd, (src + 1) * nd), k * nd);
63+ psi.set(vectors.subarray(src * nf, (src + 1) * nf), k * nf);
8164 });
82- return { freqs, psi, dof, kind: isPlanar(sys) ? "transverse" : "full" };
65+ return { freqs, psi, free };
8366 }
8467
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
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
8972 // modes by damping, which is small for light damping.
9073 export function modalDamping(sys: BeadSystem, modes: Modes): Float64Array {
9174 const { n, mass, friction, di, dj, dc } = sys;
92- const { dof, psi, freqs } = modes;
93- const nd = dof.length;
94- const u = new Float64Array(3 * n);
75+ const { free, psi, freqs } = modes;
76+ const nf = free.length;
77+ const u = new Float64Array(n);
9578 const out = new Float64Array(freqs.length);
9679 for (let k = 0; k < freqs.length; k++) {
9780 u.fill(0);
98- for (let r = 0; r < nd; r++) u[dof[r]] = psi[k * nd + r] / Math.sqrt(mass[Math.floor(dof[r] / 3)]);
81+ for (let r = 0; r < nf; r++) u[free[r]] = psi[k * nf + r] / Math.sqrt(mass[free[r]]);
9982 let c = 0;
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- }
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;
10685 out[k] = c / 2;
10786 }
10887 return out;
src/sim/parse.tsmodified+2−15View 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", "load", "relax",
90+ "duration", "sampleRate", "dt", "safety", "method", "frameRate", "modes",
9191 ]);
9292
9393 export function parseSystem(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warnings: string[] } {
@@ -170,18 +170,6 @@ 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-
185173 const hits: Hit[] = fields(f.get("hits"), "sys.hits").map((h, q) => {
186174 const nm = `sys.hits(${q + 1})`;
187175 const beads = indices(h.get("beads"), `${nm}.beads`, n);
@@ -214,8 +202,7 @@ function parseInner(sysValue: RuntimeValue | undefined): { sys: BeadSystem; warn
214202 if (method !== "verlet" && method !== "yoshida4") throw new SysError("sys.method must be 'verlet' or 'yoshida4'");
215203
216204 const sys: BeadSystem = {
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,
205+ n, pos0, mass, pinned, si, sj, sk, sL0, sk3, di, dj, dc, friction, hits, listenBeads,
219206 listenDir: [dir[0] / dn, dir[1] / dn, dir[2] / dn],
220207 duration: scalar(f.get("duration"), "sys.duration", 2),
221208 sampleRate: scalar(f.get("sampleRate"), "sys.sampleRate", 48000),
src/sim/relax.tsdeleted+0−83View file
@@ -1,83 +0,0 @@
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+0−17View file
@@ -3,7 +3,6 @@
33
44 import { executeCode, RTV } from "numbl";
55 import { parseSystem } from "./parse";
6-import { relax } from "./relax";
76 import type { SymEig } from "./modes";
87 import type { BeadSystem } from "./types";
98
@@ -21,22 +20,6 @@ export function runScript(
2120 return parseSystem(result.variableValues.sys);
2221 }
2322
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-
4023 // Symmetric eigendecomposition through numbl's eig.
4124 export const numblEig: SymEig = (A, n) => {
4225 const r = executeCode("[V, D] = eig(A); d = diag(D);", {
src/sim/solver.tsmodified+59−17View file
@@ -9,7 +9,6 @@
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";
1312 import type { BeadSystem, Method, Modes, SimResult } from "./types";
1413
1514 const YOSHIDA_W1 = 1 / (2 - Math.cbrt(2));
@@ -87,8 +86,31 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
8786 // ---- forces ------------------------------------------------------------
8887 const { si, sj, sk, sL0, sk3 } = sys;
8988 const ns = si.length;
90- // spring forces plus the constant load; returns their potential energy
91- const internalForces = (): number => springForces(sys, x, F) + loadForces(sys, x, F);
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+ };
92114 const pulses = sys.hits.filter(h => h.duration > 0);
93115 const externalForces = (t: number) => {
94116 for (const p of pulses) {
@@ -123,7 +145,7 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
123145 if (residualForce > 1e-6 * Math.max(maxTension, 1e-30)) {
124146 warnings.push(
125147 `The rest positions are not an equilibrium: the largest net force on a free bead is ${residualForce.toExponential(2)} N ` +
126- `(largest spring force ${maxTension.toExponential(2)} N). The beads will move before any strike. Set sys.relax = true to find the equilibrium first.`
148+ `(largest spring force ${maxTension.toExponential(2)} N). The beads will move before any strike.`
127149 );
128150 }
129151
@@ -274,27 +296,27 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
274296 let bv: Float64Array | null = null;
275297 let om2: Float64Array | null = null;
276298 if (modes) {
277- bz = new Float64Array(modes.dof.length);
278- bv = new Float64Array(modes.dof.length);
299+ bz = new Float64Array(modes.free.length);
300+ bv = new Float64Array(modes.free.length);
279301 om2 = new Float64Array(nModes);
280302 for (let k = 0; k < nModes; k++) om2[k] = (2 * Math.PI * modes.freqs[k]) ** 2;
281303 }
282304 const recordModes = (f: number) => {
283305 if (!modes || !modeEnergy || !bz || !bv || !om2) return;
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];
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];
291313 }
292314 const psi = modes.psi;
293315 for (let k = 0; k < nModes; k++) {
294316 let q = 0;
295317 let qd = 0;
296- const off = k * nr;
297- for (let r = 0; r < nr; r++) {
318+ const off = k * nf;
319+ for (let r = 0; r < nf; r++) {
298320 q += psi[off + r] * bz[r];
299321 qd += psi[off + r] * bv[r];
300322 }
@@ -439,9 +461,29 @@ export function simulate(sys: BeadSystem, opts: SimOptions = {}): SimResult {
439461 // Highest frequency (rad/s) of the system linearized about the rest positions,
440462 // by power iteration on M^{-1/2} H M^{-1/2}, with a small margin.
441463 export function estimateOmegaMax(sys: BeadSystem): number {
442- const { n, si, sj, mass, pinned } = sys;
464+ const { n, si, sj, sk, sL0, sk3, pos0, mass, pinned } = sys;
443465 const ns = si.length;
444- const B = springBlocks(sys);
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+ }
445487 const ism = new Float64Array(n);
446488 for (let i = 0; i < n; i++) ism[i] = pinned[i] ? 0 : 1 / Math.sqrt(mass[i]);
447489 let w = new Float64Array(3 * n);
src/sim/types.tsmodified+3−6View file
@@ -29,8 +29,6 @@ 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
3432 hits: Hit[];
3533 listenBeads: Int32Array;
3634 listenDir: [number, number, number]; // unit vector
@@ -50,11 +48,10 @@ export interface SimProgress {
5048
5149 export interface Modes {
5250 freqs: Float64Array; // Hz, ascending
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)
51+ // psi: orthonormal eigenvectors of M^{-1/2} K M^{-1/2} over free beads,
52+ // column-major (nFree x nModes)
5553 psi: Float64Array;
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
54+ free: Int32Array; // bead index of each row of psi
5855 }
5956
6057 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.</div>;
187+ if (!bands) return <div className="note">Mode energies are shown when linear modes are computed (flat systems only).</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 {p.modeKind === "full" ? "modes (3D)" : "transverse modes"}, 8 decades of color. Show the first{" "}
192+ Share of the linear (quadratic) energy in each band of 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 {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
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
351351 (frequency, T60) pairs.
352352 </div>
353353 </>
src/worker.tsmodified+8−12View 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, modeSize } from "./sim/modes";
5-import { numblEig, prepareSystem, runScript } from "./sim/runScript";
4+import { computeModes, isPlanar, modalDamping } from "./sim/modes";
5+import { numblEig, 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_DOF = 1500;
13+const MAX_MODE_BEADS = 1500;
1414 const MAX_TWIN_BYTES = 400e6;
1515
1616 export interface RunPayload {
@@ -28,8 +28,7 @@ export interface RunPayload {
2828 };
2929 result: SimResult;
3030 modeFreqs: Float64Array | null;
31- modeDecay: Float64Array | null;
32- modeKind: "transverse" | "full" | null; // small-amplitude amplitude decay rate of each mode (1/s)
31+ modeDecay: Float64Array | null; // small-amplitude amplitude decay rate of each mode (1/s)
3332 notes: string[];
3433 check: { dt: number; windows: { t: number; rel: number }[] } | null;
3534 // the same run with every hit scaled by 1e-3, its output scaled back up
@@ -65,16 +64,14 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
6564 post({ type: "stage", text: "Running script" });
6665 const { sys, warnings } = runScript(script, helpers, text => post({ type: "log", text }));
6766 const notes = [...warnings];
68- if (sys.relax) post({ type: "stage", text: "Relaxing to equilibrium" });
69- notes.push(...prepareSystem(sys));
7067
7168 let modes = null;
69+ const nFree = sys.pinned.reduce((a, p) => a + (p ? 0 : 1), 0);
7270 if (sys.modes) {
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}).`);
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}).`);
7673 else {
77- post({ type: "stage", text: `Computing linear modes (${size} degrees of freedom)` });
74+ post({ type: "stage", text: `Computing linear modes (${nFree} free beads)` });
7875 modes = computeModes(sys, numblEig);
7976 }
8077 }
@@ -156,7 +153,6 @@ self.onmessage = (ev: MessageEvent<RunRequest>) => {
156153 result,
157154 modeFreqs: modes ? modes.freqs : null,
158155 modeDecay: modes ? modalDamping(sys, modes) : null,
159- modeKind: modes ? modes.kind : null,
160156 notes,
161157 check: checkOut,
162158 linear: linearOut,
tests/screen.tsdeleted+0−34View file
@@ -1,34 +0,0 @@
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+1−42View 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, prepareSystem, runScript } from "../src/sim/runScript";
4+import { numblEig, 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,44 +258,3 @@ 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-});