Blame analysis/chaos-check.ts at main - proof-of-concept/toy-instruments

ae42884 Anti-aliased pickup: solve at R fs and low-pass filter the pickup down to fs (both solvers); CACHE_VERSION 3; chaos-check analysis scriptsJeremy Magland 12 days ago1// Is the large-pluck motion of the 5-bead string chaotic, or a solver artifact?

2//

3// Runs the post's large pluck (n=5, tension 5.001 N, bead 2, 10 mm, 2 s) with the app's own

4// model and CPU solver (double precision), and measures

5// 1. step size: runs at cfl 0.1, 0.05, ..., halving; each against the next finer one.

6// For chaos with a convergent solver, the runs agree up to a time that grows by about

7// ln(4)/lambda per halving (velocity Verlet is second order, so the error drops 4x).

8// 2. perturbations: pluck amplitude times (1 + eps), at two step sizes. The separation

9// should grow exponentially at the same rate lambda (the largest Lyapunov exponent)

10// at both step sizes, and its divergence time should grow like ln(1/eps)/lambda.

11// 3. energy balance error of every run (velocity Verlet nearly conserves energy).

12// 4. statistics: band levels of the sound (octave bands) at every step size.

13// 5. the gentle pluck (0.1 mm) against the linear modes.

14// 6. aliasing: the solver point-samples the pickup at fs (48 kHz) without a low-pass

15// filter, so sound above fs/2 folds back into the spectrogram. Runs at 16 fs and

16// reports the fraction of the pickup's power above 24 kHz.

17//

18// Positions are saved at every audio sample (frameRate = fs), since the runs separate

19// within milliseconds; the replay rate does not affect the dynamics.

20//

21// npx tsx analysis/chaos-check.ts [out.json] (about a minute per fine run)

23import { writeFileSync } from 'node:fs';

24import { CONFIGS } from '../src/instruments/configs';

25import { parseSettings, simOptions } from '../src/ui/settings';

26import { simulate } from '../src/sim/solver';

27import { spectrum } from '../src/ui/fft';

28import type { SimResult } from '../src/sim/types';

30const BASE = 'n=5&tension=5.001&excBead=2';

31const AMP = 0.01;

33interface Run {

34 label: string;

35 cfl: number;

36 amp: number;

37 r: SimResult;

38}

40function run(label: string, amp: number, cfl: number, over: { fs?: number; frameRate?: number } = { frameRate: 48000 }): Run {

41 const cfg = CONFIGS['bead-string'];

42 const s = parseSettings(cfg, `${BASE}&pluckAmp=${amp}&cfl=${cfl}`);

43 // set directly: a link only accepts the app's sample rates

44 Object.assign(s.params, over);

45 const built = cfg.setup(s.params).built!;

46 const r = simulate(built.sys, built.exc, simOptions(s.params, built));

47 const emax = Math.max(...r.energy);

48 let bal = 0;

49 for (const b of r.balance) bal = Math.max(bal, Math.abs(b));

50 console.error(`${label}: dt ${r.dt.toExponential(3)} s (${r.stepsPerSample}/sample), ${r.wallSeconds.toFixed(1)} s wall, energy error ${((bal / emax) * 100).toExponential(2)}% of peak`);

51 return { label, cfl, amp, r };

52}

54// RMS over beads of the displacement difference at each frame (m)

55function separation(a: SimResult, b: SimResult): Float64Array {

56 const nf = a.frameTimes.length;

57 const w = a.frames.length / nf;

58 const d = new Float64Array(nf);

59 for (let f = 0; f < nf; f++) {

60 let s = 0;

61 for (let i = 0; i < w; i++) {

62 const x = a.frames[f * w + i] - b.frames[f * w + i];

63 s += x * x;

64 }

65 d[f] = Math.sqrt(s / (w / 3));

66 }

67 return d;

68}

70// first time the separation exceeds thr (s), or null

71function crossing(t: Float64Array, d: Float64Array, thr: number): number | null {

72 for (let i = 0; i < d.length; i++) if (d[i] > thr) return t[i];

73 return null;

74}

76// growth rate of the running maximum of ln d between lo and hi (1/s): least squares

77function growthRate(t: Float64Array, d: Float64Array, lo: number, hi: number): { lambda: number; t0: number; t1: number } | null {

78 const xs: number[] = [];

79 const ys: number[] = [];

80 let m = 0;

81 for (let i = 0; i < d.length; i++) {

82 m = Math.max(m, d[i]);

83 if (m >= lo && m <= hi) {

84 xs.push(t[i]);

85 ys.push(Math.log(m));

86 }

87 }

88 if (xs.length < 10) return null;

89 const n = xs.length;

90 const mx = xs.reduce((s, x) => s + x, 0) / n;

91 const my = ys.reduce((s, y) => s + y, 0) / n;

92 let sxy = 0;

93 let sxx = 0;

94 for (let i = 0; i < n; i++) {

95 sxy += (xs[i] - mx) * (ys[i] - my);

96 sxx += (xs[i] - mx) ** 2;

97 }

98 return { lambda: sxy / sxx, t0: xs[0], t1: xs[n - 1] };

99}

100

101// sound level in octave bands (dB re the loudest band), from 1 s on (the settled motion)

102function bands(r: SimResult): { f: number; db: number }[] {

103 const start = Math.round(1 * r.fs);

104 const s = spectrum(r.audio.subarray(start), r.fs);

105 const out: { f: number; db: number }[] = [];

106 for (let fc = 125; fc <= 16000; fc *= 2) {

107 let p = 0;

108 for (let k = 0; k < s.freq.length; k++) if (s.freq[k] >= fc / Math.SQRT2 && s.freq[k] < fc * Math.SQRT2) p += 10 ** (s.db[k] / 10);

109 out.push({ f: fc, db: 10 * Math.log10(p) });

110 }

111 const top = Math.max(...out.map((b) => b.db));

112 return out.map((b) => ({ f: b.f, db: b.db - top }));

113}

114

115const downsample = (x: Float64Array, k: number) => Array.from(x.filter((_, i) => i % k === 0));

116

117const out: Record<string, unknown> = { settings: `${BASE}&pluckAmp=${AMP}` };

118

119// 1. step size

120const CFLS = [0.1, 0.05, 0.025, 0.0125, 0.00625, 0.003125];

121const byCfl = CFLS.map((c) => run(`cfl ${c}`, AMP, c));

122const t = byCfl[0].r.frameTimes;

123const DS = 24; // keep every 24th frame (0.5 ms) in the output

124out.time = downsample(t, DS);

125out.stepSize = byCfl.slice(0, -1).map((a, i) => {

126 const b = byCfl[i + 1];

127 const d = separation(a.r, b.r);

128 return {

129 cfl: a.cfl,

130 against: b.cfl,

131 dt: a.r.dt,

132 divergeAt1um: crossing(t, d, 1e-6),

133 divergeAt1mm: crossing(t, d, 1e-3),

134 d: downsample(d, DS),

135 };

136});

137

138// 2. perturbations at two step sizes

139const EPS = [1e-5, 1e-7, 1e-9, 1e-11];

140const pert: Record<string, unknown>[] = [];

141for (const [ci, cfl] of [[0, 0.1], [2, 0.025]] as const) {

142 const base = byCfl[ci];

143 for (const eps of EPS) {

144 const p = run(`cfl ${cfl}, eps ${eps}`, AMP * (1 + eps), cfl);

145 const d = separation(base.r, p.r);

146 pert.push({

147 cfl,

148 eps,

149 divergeAt1mm: crossing(t, d, 1e-3),

150 growth: growthRate(t, d, 1e-8, 1e-4),

151 bands: bands(p.r),

152 d: downsample(d, DS),

153 });

154 }

155}

156out.perturbation = pert;

157

158// 3. energy balance and 4. sound statistics per step size

159out.runs = byCfl.map(({ cfl, r }) => {

160 const emax = Math.max(...r.energy);

161 let bal = 0;

162 for (const b of r.balance) bal = Math.max(bal, Math.abs(b));

163 return { cfl, dt: r.dt, stepsPerSample: r.stepsPerSample, wallSeconds: r.wallSeconds, energyError: bal / emax, maxStretch: r.maxStretch, bands: bands(r) };

164});

165

166// 5. gentle pluck vs the linear modes: f_k = (1/pi) sqrt(T / (m a)) sin(k pi / 12)

167const gentle = run('gentle 0.1 mm', 0.0001, 0.1);

168{

169 const s = spectrum(gentle.r.audio, gentle.r.fs);

170 const T = 5.001;

171 const m = 5.08e-5;

172 const a = 0.0127;