Lyapunov method in src/common/lyapunov.ts; each app supplies its Shadow
12 changed files+206−210
README.mdmodified+1−1View file
@@ -4,7 +4,7 @@ Interactive apps for "Listening for chaos": struck dynamical systems heard as so
44
55 - `src/double-pendulum/`: double pendulum ([README](src/double-pendulum/README.md)), served at `/double-pendulum/`.
66 - `src/fput-chain/`: FPUT-beta chain ([README](src/fput-chain/README.md)), served at `/fput-chain/`.
7-- `src/common/`: code shared by the apps (help dialogs, chart cursor, WAV export, λ vs E scatter, strip charts, whole-note spectrogram) and `style.css`, the shared styles. Each app's `main.ts` imports `../common/style.css` and then its own `style.css`, so app rules come after the shared ones in dev and in the build.
7+- `src/common/`: code shared by the apps (help dialogs, chart cursor, WAV export, λ vs E scatter, strip charts, whole-note spectrogram, Lyapunov exponent by the two-trajectory method with a per-app `Shadow`) and `style.css`, the shared styles. Each app's `main.ts` imports `../common/style.css` and then its own `style.css`, so app rules come after the shared ones in dev and in the build.
88 - `src/index.html`: landing page linking the apps.
99
1010 Live site: https://chaos-apps--listening-for-chaos.vault1-sites.magland.org/ (built and deployed by `.github/workflows/pages.yml` on every push to main).
src/common/lyapunov.tsadded+108−0View file
@@ -0,0 +1,108 @@
1+// Largest Lyapunov exponent by the two-trajectory (Benettin) method. A shadow
2+// copy of the system is kept at distance D0 from it: after every step the
3+// separation is measured, its logarithmic growth recorded, and the shadow
4+// moved back to distance D0 along the current separation direction. The
5+// exponent is the average growth rate: over a sliding window (finite-time
6+// exponent, which follows a damped ring-down) and since the start.
7+//
8+// Each app supplies the shadow (a Shadow): how it is stepped, how far it is
9+// from the main state (in the app's scaled coordinates, so that positions
10+// and momenta count comparably), and how it is moved back. The long-time
11+// exponent does not depend on the scaling; finite-time values depend on it a
12+// little.
13+
14+export const D0 = 1e-8; // shadow distance
15+
16+export interface Shadow {
17+ // Advance the shadow by h (called right after the main state's step).
18+ step(h: number): void;
19+ // Distance from the main state, in scaled coordinates.
20+ separation(): number;
21+ // Move the shadow to main + r (shadow - main), plus whatever the system
22+ // needs after a jump (forces at the new positions, energy bookkeeping).
23+ pullBack(r: number): void;
24+}
25+
26+export class Lyapunov {
27+ private window: number; // s
28+ private shadow: Shadow | null = null;
29+ // Ring of (time, log growth) per step, for the sliding window (grown as
30+ // needed; not kept at all for an infinite window).
31+ private ts = new Float64Array(1024);
32+ private gs = new Float64Array(1024);
33+ private head = 0;
34+ private tail = 0;
35+ private windowSum = 0;
36+ private totalSum = 0;
37+ private t0 = 0; // start time
38+ private t = 0; // time of the last step
39+
40+ constructor(window: number) {
41+ this.window = window;
42+ }
43+
44+ // Start (or restart) measuring with a fresh shadow (at distance D0) at time t.
45+ start(shadow: Shadow, t: number): void {
46+ this.shadow = shadow;
47+ this.head = this.tail = 0;
48+ this.windowSum = this.totalSum = 0;
49+ this.t0 = this.t = t;
50+ }
51+
52+ stop(): void {
53+ this.shadow = null;
54+ }
55+
56+ // Advance the shadow by h (call right after the main state's step) and
57+ // measure; t is the main state's time after the step.
58+ step(h: number, t: number): void {
59+ const s = this.shadow;
60+ if (!s) return;
61+ s.step(h);
62+ const d = s.separation();
63+ s.pullBack(D0 / d);
64+
65+ const g = Math.log(d / D0);
66+ this.t = t;
67+ this.totalSum += g;
68+ if (this.window === Infinity) return;
69+ this.windowSum += g;
70+ this.push(t, g);
71+ while (this.tail !== this.head && this.ts[this.tail] <= t - this.window) {
72+ this.windowSum -= this.gs[this.tail];
73+ this.tail = (this.tail + 1) % this.ts.length;
74+ }
75+ }
76+
77+ // Exponent over the last `window` seconds (or since the start, if shorter), 1/s.
78+ windowed(): number {
79+ if (this.window === Infinity) return this.mean();
80+ const span = Math.min(this.window, this.t - this.t0);
81+ return this.shadow && span > 0 ? this.windowSum / span : NaN;
82+ }
83+
84+ // Exponent since the start, 1/s.
85+ mean(): number {
86+ const span = this.t - this.t0;
87+ return this.shadow && span > 0 ? this.totalSum / span : NaN;
88+ }
89+
90+ private push(t: number, g: number): void {
91+ const n = this.ts.length;
92+ if ((this.head + 1) % n === this.tail) {
93+ // Full: copy into rings twice the size, oldest first.
94+ const ts = new Float64Array(2 * n), gs = new Float64Array(2 * n);
95+ for (let i = 0, j = this.tail; j !== this.head; i++, j = (j + 1) % n) {
96+ ts[i] = this.ts[j];
97+ gs[i] = this.gs[j];
98+ }
99+ this.tail = 0;
100+ this.head = n - 1;
101+ this.ts = ts;
102+ this.gs = gs;
103+ }
104+ this.ts[this.head] = t;
105+ this.gs[this.head] = g;
106+ this.head = (this.head + 1) % this.ts.length;
107+ }
108+}
src/double-pendulum/README.mdmodified+2−2View file
@@ -14,10 +14,10 @@ Live site: https://chaos-apps--listening-for-chaos.vault1-sites.magland.org/doub
1414 - Energy is zero at hanging rest.
1515 - Damping: linear drag on each bob, force -gamma m_i v_i (proportional to its mass; equal for equal bobs), so the generalized forces are -gamma p and every mode's amplitude decays with time constant tau = 2/gamma. Default tau = 540 s: about 10% amplitude per minute in the animation, about 2 s in the sound at x262. The state carries the dissipated energy W, so E + W is checked for conservation. With damping the energy slider (labelled "Strike energy") sets the strike energy and the current energy is shown at the top left of the animation; with tau = 0 it is labelled "Energy" and sets the energy live. The section's controls are in two rows (energy and Follow; Color by λ and Click to strike), and "computing…" shows in the plot's top-right legend.
1616 - Live energy change: both momenta are scaled by one factor at fixed angles (a push along the current motion). A target below the current potential energy can't be reached at once; until it can, kinetic energy is removed at a rate of 6 per second.
17-- Lyapunov exponent: two-trajectory (Benettin) method in `lyapunov.ts`. A hidden shadow pendulum is kept 1e-8 away (angles in rad, momenta in units of sqrt((m1 + m2) l1^2 E_top)), moved back after every step, and the log growth averaged. The shadow keeps its own energy (forcing it onto the pendulum's energy by rescaling momenta inflates the separation near turning points). Checked (no damping, 300 s): 2 J and 10 J fall toward 0 (0.01/s), 30 J ~1.0/s, 50 J ~1.4/s; unchanged for shadow distances 1e-6 and 1e-10. With damping (tau 540 s) a 50 J strike's 10 s value falls from ~1.4 to ~0.1/s by 300 s.
17+- Lyapunov exponent: two-trajectory (Benettin) method in `../common/lyapunov.ts`, with the pendulum's shadow in `lyapunov.ts`. A hidden shadow pendulum is kept 1e-8 away (angles in rad, momenta in units of sqrt((m1 + m2) l1^2 E_top)), moved back after every step, and the log growth averaged. The shadow keeps its own energy (forcing it onto the pendulum's energy by rescaling momenta inflates the separation near turning points). Checked (no damping, 300 s): 2 J and 10 J fall toward 0 (0.01/s), 30 J ~1.0/s, 50 J ~1.4/s; unchanged for shadow distances 1e-6 and 1e-10. With damping (tau 540 s) a 50 J strike's 10 s value falls from ~1.4 to ~0.1/s by 300 s.
1818 - Poincare section (`poincare-core.ts`): on the section th1 = 0, a point (th2, p2) is reachable at energy E exactly when p2^2 <= 2 m2 l2^2 (E - V(th2)), V(th2) = m2 g l2 (1 - cos th2); p1 is the root with the inner arm moving right. Crossings are located by linear interpolation within a step. The step is 4x the animation's: over 300 crossings the energy error stays below ~2e-5, the regular/chaotic classification is unchanged, and regular orbits move by ~1e-4. Each orbit's exponent is the shadow method's mean over its run (regular orbits ~0.01/s or less).
1919 - Integrator: 2-stage Gauss-Legendre (order 4, symplectic), fixed step h chosen per parameter set by a short calibration run so the relative energy error stays below 1e-6 up to the top of the energy slider (twice the both-inverted energy); 3.13 ms for the defaults, fixed-point stage solve (started from the previous step's stage derivatives extrapolated in time; tolerance 1e-12, below which the energy error does not change). Angles are wrapped to (-pi, pi] after each step. The physics costs well under 1 ms of CPU per second; drawing dominates, so traces draw at most a min/max pair per pixel column.
2020
2121 Audio: a Web Worker (`audio-worker.ts`) integrates the audio copy with substeps of at most 3.2 ms of simulated time per sample, high-passes at 20 Hz, and streams 1024-sample chunks through a MessageChannel to an AudioWorklet (`audio-worklet.ts`) that plays them from a ring buffer, keeping about 0.25 s queued. The page shows how many times faster than real time the worker computes, and any dropouts. The worker also sends each chunk to the page, which feeds it to the audio spectrogram as the audio clock reaches it (2048-sample window, new column every 512 samples). The worker runs the same Lyapunov shadow on the note and sends one exponent and one energy value per 256 samples with each chunk; the page plots them in step with the spectrogram. The whole-note spectrogram keeps every column (pairs are averaged past 4096 columns) and draws only new pixel columns unless the axis or color scale changes.
2222
23-Code: `pendulum.ts` (physics), `sim.ts` (running simulation with energy target), `audio.ts` (audio setup), `signals.ts` (audible quantities), `lyapunov.ts` (Lyapunov exponent), `../common/scatter.ts` (lambda vs energy), `../common/cursor.ts` (shared hover/pin cursor), `poincare-core.ts`, `poincare-worker.ts`, `poincare.ts` (Poincare section), `atlas.ts` (cached sections for following the sound), `audio-shared.ts` (constants shared with the worker), `../common/traces.ts` (strip charts), `../common/spectrogram.ts` (whole-note spectrogram), `main.ts` (page).
23+Code: `pendulum.ts` (physics), `sim.ts` (running simulation with energy target), `audio.ts` (audio setup), `signals.ts` (audible quantities), `lyapunov.ts` (the Lyapunov shadow; the method is `../common/lyapunov.ts`), `../common/scatter.ts` (lambda vs energy), `../common/cursor.ts` (shared hover/pin cursor), `poincare-core.ts`, `poincare-worker.ts`, `poincare.ts` (Poincare section), `atlas.ts` (cached sections for following the sound), `audio-shared.ts` (constants shared with the worker), `../common/traces.ts` (strip charts), `../common/spectrogram.ts` (whole-note spectrogram), `main.ts` (page).
src/double-pendulum/audio-worker.tsmodified+4−3View file
@@ -13,7 +13,8 @@
1313
1414 import { energy, normalModes, stepSize, topEnergy, type Params } from "./pendulum";
1515 import { END_RATIO, LAMBDA_EVERY } from "./audio-shared";
16-import { Lyapunov } from "./lyapunov";
16+import { Lyapunov } from "../common/lyapunov";
17+import { shadowOf } from "./lyapunov";
1718 import { signalByKey, type Signal } from "./signals";
1819 import { Sim } from "./sim";
1920
@@ -96,7 +97,7 @@ onmessage = (e: MessageEvent<ToWorker>) => {
9697 session++;
9798 sim.setState(msg.y, msg.E0, msg.target);
9899 eStart = energy(sim.P, sim.y);
99- lyap.start(sim);
100+ lyap.start(shadowOf(sim), sim.t);
100101 produced = consumed = underruns = 0;
101102 computeMs = computedSamples = 0;
102103 lastReport = performance.now();
@@ -140,7 +141,7 @@ function pump(): void {
140141 for (let k = 0; k < n; k++) {
141142 const a = sim.y[0], b = sim.y[1], q = sim.y[3];
142143 sim.step(h);
143- lyap.step(sim, h);
144+ lyap.step(h, sim.t);
144145 if (a < 0 && sim.y[0] >= 0 && sim.y[0] - a < 1) {
145146 const f = -a / (sim.y[0] - a);
146147 const t2 = b + f * Math.atan2(Math.sin(sim.y[1] - b), Math.cos(sim.y[1] - b));
src/double-pendulum/lyapunov.tsmodified+37−94View file
@@ -1,112 +1,55 @@
1-// Largest Lyapunov exponent by the two-trajectory (Benettin) method. A shadow
2-// copy of the pendulum is kept at distance D0 from it: after every step the
3-// separation is measured, its logarithmic growth recorded, and the shadow
4-// moved back to distance D0 along the current separation direction. The
5-// exponent is the average growth rate: over a sliding window (finite-time
6-// exponent, which follows a damped ring-down) and since the start.
1+// The pendulum's Lyapunov shadow (see ../common/lyapunov.ts): a second Sim
2+// kept D0 from the pendulum.
73 //
84 // Distance: angles in rad (difference wrapped to (-pi, pi]) and momenta in
95 // units of pRef = sqrt((m1 + m2) l1^2 E_top), so both count comparably at high
10-// energy. The long-time exponent does not depend on this choice; finite-time
11-// values depend on it a little.
6+// energy.
127
8+import { D0, type Shadow } from "../common/lyapunov";
139 import { energy, topEnergy } from "./pendulum";
1410 import { Sim } from "./sim";
1511
16-const D0 = 1e-8; // shadow distance
17-const CAP = 1 << 16; // per-step growth values kept (> window / step)
18-
1912 const wrap = (a: number) => Math.atan2(Math.sin(a), Math.cos(a));
2013
21-export class Lyapunov {
22- private window: number; // s
23- private shadow: Sim | null = null;
24- private pRef = 1;
25- // Ring of (time, log growth) per step, for the sliding window.
26- private ts = new Float64Array(CAP);
27- private gs = new Float64Array(CAP);
28- private head = 0;
29- private tail = 0;
30- private windowSum = 0;
31- private totalSum = 0;
32- private t0 = 0; // start time
33- private t = 0; // time of the last step
34-
35- constructor(window: number) {
36- this.window = window;
37- }
38-
39- // Start (or restart) measuring from the pendulum's current state.
40- start(sim: Sim): void {
41- const P = sim.P;
42- this.pRef = Math.sqrt((P.m1 + P.m2) * P.l1 * P.l1 * topEnergy(P));
43- const y = Float64Array.from(sim.y);
44- // Initial separation: equal parts in all four coordinates.
45- y[0] += D0 / 2;
46- y[1] += D0 / 2;
47- y[2] += (D0 / 2) * this.pRef;
48- y[3] += (D0 / 2) * this.pRef;
49- this.shadow = new Sim(P, 0, 0.5);
50- this.shadow.setState(y, 0, 0);
51- this.matchEnergy(sim);
52- this.head = this.tail = 0;
53- this.windowSum = this.totalSum = 0;
54- this.t0 = this.t = sim.t;
55- }
56-
57- stop(): void {
58- this.shadow = null;
59- }
60-
61- // Advance the shadow by h (call right after sim.step(h)) and measure.
62- step(sim: Sim, h: number): void {
63- const s = this.shadow;
64- if (!s) return;
65- s.step(h);
66- const a = sim.y, b = s.y, pr = this.pRef;
67- const d0 = wrap(b[0] - a[0]), d1 = wrap(b[1] - a[1]);
68- const d2 = (b[2] - a[2]) / pr, d3 = (b[3] - a[3]) / pr;
69- const d = Math.hypot(d0, d1, d2, d3);
70- const r = D0 / d;
71- b[0] = a[0] + d0 * r;
72- b[1] = a[1] + d1 * r;
73- b[2] = a[2] + d2 * r * pr;
74- b[3] = a[3] + d3 * r * pr;
75- this.matchEnergy(sim);
76-
77- const g = Math.log(d / D0);
78- this.t = sim.t;
79- this.totalSum += g;
80- this.windowSum += g;
81- this.ts[this.head] = sim.t;
82- this.gs[this.head] = g;
83- this.head = (this.head + 1) % CAP;
84- while (this.tail !== this.head && this.ts[this.tail] <= sim.t - this.window) {
85- this.windowSum -= this.gs[this.tail];
86- this.tail = (this.tail + 1) % CAP;
87- }
88- }
89-
90- // Exponent over the last `window` seconds (or since the start, if shorter), 1/s.
91- windowed(): number {
92- const span = Math.min(this.window, this.t - this.t0);
93- return this.shadow && span > 0 ? this.windowSum / span : NaN;
94- }
95-
96- // Exponent since the start, 1/s.
97- mean(): number {
98- const span = this.t - this.t0;
99- return this.shadow && span > 0 ? this.totalSum / span : NaN;
100- }
14+// A shadow of the pendulum's current state, for Lyapunov.start.
15+export function shadowOf(sim: Sim): Shadow {
16+ const P = sim.P;
17+ const pr = Math.sqrt((P.m1 + P.m2) * P.l1 * P.l1 * topEnergy(P));
18+ const y = Float64Array.from(sim.y);
19+ // Initial separation: equal parts in all four coordinates.
20+ y[0] += D0 / 2;
21+ y[1] += D0 / 2;
22+ y[2] += (D0 / 2) * pr;
23+ y[3] += (D0 / 2) * pr;
24+ const s = new Sim(P, 0, 0.5);
25+ s.setState(y, 0, 0);
10126
10227 // The shadow keeps its own energy (forcing it onto the pendulum's exact
10328 // energy would rescale momenta, which near a turning point, T ~ 0, injects
10429 // a large artificial separation). Without damping it follows changes of the
10530 // pendulum's energy target by the same amount.
106- private matchEnergy(sim: Sim): void {
107- const s = this.shadow!;
31+ const matchEnergy = () => {
10832 const E = energy(s.P, s.y);
10933 s.E0 = E;
11034 s.target = s.P.tau > 0 ? E : E + (sim.target - sim.E0);
111- }
35+ };
36+ matchEnergy();
37+
38+ return {
39+ step: (h) => s.step(h),
40+ separation: () => {
41+ const a = sim.y, b = s.y;
42+ return Math.hypot(wrap(b[0] - a[0]), wrap(b[1] - a[1]), (b[2] - a[2]) / pr, (b[3] - a[3]) / pr);
43+ },
44+ pullBack: (r) => {
45+ const a = sim.y, b = s.y;
46+ const d0 = wrap(b[0] - a[0]), d1 = wrap(b[1] - a[1]);
47+ const d2 = (b[2] - a[2]) / pr, d3 = (b[3] - a[3]) / pr;
48+ b[0] = a[0] + d0 * r;
49+ b[1] = a[1] + d1 * r;
50+ b[2] = a[2] + d2 * r * pr;
51+ b[3] = a[3] + d3 * r * pr;
52+ matchEnergy();
53+ },
54+ };
11255 }
src/double-pendulum/main.tsmodified+4−3View file
@@ -5,7 +5,8 @@ import { END_RATIO, LAMBDA_EVERY } from "./audio-shared";
55 import { newCursor } from "../common/cursor";
66 import { setupHelp } from "../common/help";
77 import { PoincarePanel, formatEnergy } from "./poincare";
8-import { Lyapunov } from "./lyapunov";
8+import { Lyapunov } from "../common/lyapunov";
9+import { shadowOf } from "./lyapunov";
910 import { ScatterPlot } from "../common/scatter";
1011 import {
1112 energy, kinetic, normalModes, potential, stepSize, topEnergy, type Params,
@@ -498,7 +499,7 @@ function afterStrike() {
498499 poincare.clearLive();
499500 if (poincare.followMode() === "pendulum") poincare.follow(energy(P, sim.y), true);
500501 resetGhost();
501- lyap.start(sim);
502+ lyap.start(shadowOf(sim), sim.t);
502503 setRunning(true);
503504 pauseButton.disabled = false;
504505 traces.clear();
@@ -686,7 +687,7 @@ function frame(now: number) {
686687 poincare.addLive(Math.atan2(Math.sin(t2), Math.cos(t2)), q + f * (sim.y[3] - q), energy(P, sim.y));
687688 }
688689 ghost?.step(H);
689- lyap.step(sim, H);
690+ lyap.step(H, sim.t);
690691 pending -= H;
691692 traces.sample(sim.t);
692693 }
src/double-pendulum/poincare-core.tsmodified+4−3View file
@@ -10,7 +10,8 @@
1010 // with th1dot = sqrt(disc)/det > 0 is p1 = (c p2 + sqrt(disc)) / b.
1111
1212 import { energy, type Params } from "./pendulum";
13-import { Lyapunov } from "./lyapunov";
13+import { Lyapunov } from "../common/lyapunov";
14+import { shadowOf } from "./lyapunov";
1415 import { Sim } from "./sim";
1516
1617 export const CROSSINGS = 300; // crossings recorded per orbit
@@ -66,7 +67,7 @@ export class SectionOrbit {
6667 this.sim = new Sim({ ...P, tau: 0 }, 0, 0.5);
6768 const E0 = energy(this.sim.P, y);
6869 this.sim.setState(y, E0, E0);
69- this.lyap.start(this.sim);
70+ this.lyap.start(shadowOf(this.sim), this.sim.t);
7071 this.h = h;
7172 this.maxTime = maxTime;
7273 this.points = [th2, p2];
@@ -85,7 +86,7 @@ export class SectionOrbit {
8586 while (sim.t < this.maxTime) {
8687 const th1 = y[0], th2 = y[1], p2 = y[3];
8788 sim.step(h);
88- this.lyap.step(sim, h);
89+ this.lyap.step(h, sim.t);
8990 if (th1 < 0 && y[0] >= 0 && y[0] - th1 < 1) {
9091 const f = -th1 / (y[0] - th1);
9192 const d2 = Math.atan2(Math.sin(y[1] - th2), Math.cos(y[1] - th2));
src/fput-chain/README.mdmodified+1−1View file
@@ -2,7 +2,7 @@
22
33 Struck FPUT-beta chain, step 2 of "Listening for chaos" (after double-pendulum-chaos). Work in progress.
44
5-Current state: N masses (2 to 64, default 2) between two fixed walls, joined by springs with tension k d + k3 d^3. The chain sits at rest until Strike, which gives it energy E (log slider, 10 µJ to 100 J, default 2.51 J) as velocity, either on one mass j alone (default mass 1, a hammer blow at one end) or in the shape of one normal mode n. The left card animates the chain (top: to scale, rest spacing a; below: the same displacements drawn sideways as a profile, scaled to fit, with the scale given) with the time and the relative energy error. Speed: 0.1x to 20x real time (default 1x). The time and current energy show at the top left. The Ghost chip (as in the pendulum app) adds a second chain started from the same state with every mass moved by the same amount, 10^k / sqrt(N) m (total offset 10^k m, k = -12 to -1, default -10; moving only the struck mass would put the offset along the motion itself, which does not grow), run in lockstep without renormalization and drawn in translucent orange underneath; restarted on Strike, on "Restart ghost" (from the chain's current state), when switched on or when the offset changes; a readout gives the separation (root sum of squared displacement differences, m) and the time since the ghost started, and a "Separation" chart (log, last 60 s; offered only with the ghost on) plots it. A "Lyapunov" chart (last 60 s; the first 5 s after a strike left out of its y range) shows the largest Lyapunov exponent over the last 10 s and since the strike (`lyapunov.ts`, two-trajectory Benettin method as in the pendulum app: a hidden shadow chain 1e-8 away, displacements in units of d0 and momenta in units of sqrt(m k) d0, renormalized every step). Checked (N = 4, mass-1 strike, 300 s): without damping 0.1 J ~0.05 to 0.1 /s (weak chaos; N = 4 onset is near 3 k d0^2 = 75 mJ), 0.5 J ~0.5 /s, 2.51 J ~0.9 /s, 25 J ~1.7 /s; with tau 540 s, 0.1 J falls to ~0.01 /s and 2.51 J falls from ~1.0 to ~0.3 /s by 300 s. The slider is labelled "Strike energy". Parameters (m, k, k3, a, tau) fold out underneath with derived values (mode frequencies, d0, nonlinear energy scale, integrator step); Apply or changing N returns the chain to rest. The "Real time" card on the right has toggle chips for its charts. "Mode history" (off by default) is a scrolling heat map of each mode's share of the harmonic energy, log10(E_n / sum E) over 6 decades (inferno colors, mode 1 at the bottom), over the last 300 s (`LiveModeHistory` in `modehistory.ts`); the cascade of energy up the modes shows here (e.g. N = 32, mode-1 strike at 10 mJ: the top mode's share goes from ~1e-28 to ~1e-12 in 10 periods of mode 1 and ~3e-5 in 50 without damping; with tau 540 s the damping wins and the spread stalls near mode 17). "Mode bars" shows the current harmonic energy of each normal mode relative to E on a log scale (12 decades), with the number of modes sharing the energy (inverse participation, 1 / sum (E_n / sum E)^2). The other chips are charts over time (from `../common/traces.ts` and `../common/cursor.ts`, with per-group colors and a legend capped at 8 series): displacements, velocities, spring stretches, mode energies (log, last 60 s; values below 1e-12 of E left out), energies (E, T, V, W dissipated) and the energy balance error ((E + W - E0) / E0); others last 10 s; sampled every 10 ms; hover for values, click to pin. Mode bars, displacements and mode energies are on by default. Masses, springs and modes share one color per index (four fixed colors for up to 4, a blue-to-red sweep beyond).
5+Current state: N masses (2 to 64, default 2) between two fixed walls, joined by springs with tension k d + k3 d^3. The chain sits at rest until Strike, which gives it energy E (log slider, 10 µJ to 100 J, default 2.51 J) as velocity, either on one mass j alone (default mass 1, a hammer blow at one end) or in the shape of one normal mode n. The left card animates the chain (top: to scale, rest spacing a; below: the same displacements drawn sideways as a profile, scaled to fit, with the scale given) with the time and the relative energy error. Speed: 0.1x to 20x real time (default 1x). The time and current energy show at the top left. The Ghost chip (as in the pendulum app) adds a second chain started from the same state with every mass moved by the same amount, 10^k / sqrt(N) m (total offset 10^k m, k = -12 to -1, default -10; moving only the struck mass would put the offset along the motion itself, which does not grow), run in lockstep without renormalization and drawn in translucent orange underneath; restarted on Strike, on "Restart ghost" (from the chain's current state), when switched on or when the offset changes; a readout gives the separation (root sum of squared displacement differences, m) and the time since the ghost started, and a "Separation" chart (log, last 60 s; offered only with the ghost on) plots it. A "Lyapunov" chart (last 60 s; the first 5 s after a strike left out of its y range) shows the largest Lyapunov exponent over the last 10 s and since the strike (`../common/lyapunov.ts`, two-trajectory Benettin method shared with the pendulum app, with the chain's shadow in `lyapunov.ts`: a hidden shadow chain 1e-8 away, displacements in units of d0 and momenta in units of sqrt(m k) d0, renormalized every step). Checked (N = 4, mass-1 strike, 300 s): without damping 0.1 J ~0.05 to 0.1 /s (weak chaos; N = 4 onset is near 3 k d0^2 = 75 mJ), 0.5 J ~0.5 /s, 2.51 J ~0.9 /s, 25 J ~1.7 /s; with tau 540 s, 0.1 J falls to ~0.01 /s and 2.51 J falls from ~1.0 to ~0.3 /s by 300 s. The slider is labelled "Strike energy". Parameters (m, k, k3, a, tau) fold out underneath with derived values (mode frequencies, d0, nonlinear energy scale, integrator step); Apply or changing N returns the chain to rest. The "Real time" card on the right has toggle chips for its charts. "Mode history" (off by default) is a scrolling heat map of each mode's share of the harmonic energy, log10(E_n / sum E) over 6 decades (inferno colors, mode 1 at the bottom), over the last 300 s (`LiveModeHistory` in `modehistory.ts`); the cascade of energy up the modes shows here (e.g. N = 32, mode-1 strike at 10 mJ: the top mode's share goes from ~1e-28 to ~1e-12 in 10 periods of mode 1 and ~3e-5 in 50 without damping; with tau 540 s the damping wins and the spread stalls near mode 17). "Mode bars" shows the current harmonic energy of each normal mode relative to E on a log scale (12 decades), with the number of modes sharing the energy (inverse participation, 1 / sum (E_n / sum E)^2). The other chips are charts over time (from `../common/traces.ts` and `../common/cursor.ts`, with per-group colors and a legend capped at 8 series): displacements, velocities, spring stretches, mode energies (log, last 60 s; values below 1e-12 of E left out), energies (E, T, V, W dissipated) and the energy balance error ((E + W - E0) / E0); others last 10 s; sampled every 10 ms; hover for values, click to pin. Mode bars, displacements and mode energies are on by default. Masses, springs and modes share one color per index (four fixed colors for up to 4, a blue-to-red sweep beyond).
66
77 "Sped up · listen" (open by default, below the cards; otherwise as in the pendulum app): the header says how much faster it runs ("the whole decay, 321× faster: 62 min of motion in 12 s"; without damping "321× faster than the chain above"). The page starts muted (the mute button is tinted red and reads "Muted"; muted notes are still computed and shown); the audio context is created at the first Strike; closing the panel stops the note. While it is open each Strike plays a note, the same strike rendered in a Web Worker (`render-worker.ts`) sped up so the lowest linear mode sounds at the Pitch (5 to 1000 Hz, default 100; below ~20 Hz the lowest modes are below hearing but show in the spectrogram; the readout gives the lowest and highest mode pitches). A damped note is the whole decay, ending when its energy has fallen to 1e-6 of the strike's (-60 dB; predicted length (tau/2) ln 1e6 of physical time, 11.6 s of sound by default, measured 11.3 s; capped at twice that); undamped notes last 10 s. "Strike to hear the whole decay." shows until the first note; the status line gives the signal and note length ("computing…" while rendering). Replay plays the last strike again with the current settings (to compare signals). Signal (`signals.ts`): velocity (default), displacement or acceleration of mass j (default 1, the end mass, which is never at a node), force on the left wall (the first spring's tension), or sum of accelerations (the net wall force, in which even modes cancel). In one mass's motion mode n is weighted by sqrt(E_n) / omega_n, sqrt(E_n) or sqrt(E_n) omega_n. Each note is made zero-mean, scaled to a 0.9 peak (all signals equally loud; Volume sets the level) and faded out over its last 20 ms (`sound.ts`). Save WAV saves the last note as heard (16-bit mono PCM, peak at 0.9; `../common/wav.ts`; named e.g. fput-N4-mass1-2.51J-v1-100Hz.wav). Chips (as in the pendulum app, plus Modes): Spectrogram (on; `NoteSpectrogram` in `../common/spectrogram.ts`: 2048-sample Hann, zero-pad 4096, hop 512, for pitches down to 30 Hz; below, the window lengthens in proportion (up to 16384 samples at 5 Hz, zero-pad 2x, hop a quarter window) so neighbouring low modes stay resolved; dashed lines at the linear mode pitches, drawn only while at least 6 px apart, since crowded together their dashes merged into a solid on-off bar at the bottom), Modes (off by default; the note's mode-energy history, each mode's share of the energy over the whole note, from the render worker, one column per 256 samples; `ModeHistory` in `modehistory.ts`; a cascade that takes tens of mode-1 periods occupies only the first few percent of a whole-decay note, so Mode history in the Real time panel shows it better), Lyapunov (the note's largest Lyapunov exponent over the last 10 s of physical time, from a shadow run in the worker alongside the note; the first 5 s of physical time left out of the y range), Energy (log scale, so the decay is a straight line) and λ vs E (`../common/scatter.ts`: one point per 256 samples, accumulated over notes, log energy axis, colored by time into the note on a log scale, earlier notes faded, hover for values, Clear button). The first three span the whole note (strike at left, end at right) and fill in as the audio clock reaches each sample; one shared cursor (playhead while playing; hover for values, click to pin). Pause (or space) pauses the animation and the sound (audio context suspended). The note's Lyapunov shadow runs only while the Lyapunov chart or λ vs E is shown (it doubles the computation). Render time in the browser pane (with the Lyapunov shadow): ~0.6 s for N = 4, ~2 s for N = 32 (default strike). First look (N = 4, 2.51 J): λ ~1.3 /s at the strike, falling to ~0 by ~4 s of sound as the spectrogram resolves into lines; on λ vs E the note's points sit near 0 below ~0.1 J and rise to ~1.4 /s at the top.
88
src/fput-chain/lyapunov.tsmodified+33−94View file
@@ -1,102 +1,41 @@
1-// Largest Lyapunov exponent by the two-trajectory (Benettin) method, as in
2-// the double pendulum app. A shadow copy of the chain is kept at distance D0
3-// from it: after every step the separation is measured, its logarithmic
4-// growth recorded, and the shadow moved back to distance D0 along the
5-// current separation direction. The exponent is the average growth rate:
6-// over a sliding window (finite-time exponent, which follows a damped
7-// ring-down) and since the start.
1+// The chain's Lyapunov shadow (see ../common/lyapunov.ts): a second chain
2+// kept D0 from the chain.
83 //
94 // Distance: displacements in units of d0 = sqrt(k / k3) and momenta in units
105 // of sqrt(m k) d0 (the momentum of a mass swinging at sqrt(k / m) with
11-// amplitude d0), so both count comparably. The long-time exponent does not
12-// depend on this choice; finite-time values depend on it a little.
6+// amplitude d0), so both count comparably.
137
8+import { D0, type Shadow } from "../common/lyapunov";
149 import { Verlet, copyState, d0, type Params, type State } from "./chain";
1510
16-const D0 = 1e-8; // shadow distance
17-const CAP = 1 << 18; // per-step growth values kept (> window / step)
18-
19-export class Lyapunov {
20- private window: number; // s
21- private shadow: { s: State; verlet: Verlet } | null = null;
22- private uRef = 1;
23- private pRef = 1;
24- // Ring of (time, log growth) per step, for the sliding window.
25- private ts = new Float64Array(CAP);
26- private gs = new Float64Array(CAP);
27- private head = 0;
28- private tail = 0;
29- private windowSum = 0;
30- private totalSum = 0;
31- private t0 = 0; // start time
32- private t = 0; // time of the last step
33-
34- constructor(window: number) {
35- this.window = window;
36- }
37-
38- // Start (or restart) measuring from the chain's current state at time t.
39- start(P: Params, s: State, t: number): void {
40- this.uRef = d0(P);
41- this.pRef = Math.sqrt(P.m * P.k) * this.uRef;
42- const sh = copyState(s);
43- // Initial separation: equal parts in all 2N coordinates.
44- const e = D0 / Math.sqrt(2 * P.N);
45- for (let j = 0; j < P.N; j++) {
46- sh.u[j] += e * this.uRef;
47- sh.p[j] += e * this.pRef;
48- }
49- this.shadow = { s: sh, verlet: new Verlet(P, sh) };
50- this.head = this.tail = 0;
51- this.windowSum = this.totalSum = 0;
52- this.t0 = this.t = t;
53- }
54-
55- stop(): void {
56- this.shadow = null;
57- }
58-
59- // Advance the shadow by h (call right after the chain's step) and measure;
60- // t is the chain's time after the step.
61- step(s: State, h: number, t: number): void {
62- const sh = this.shadow;
63- if (!sh) return;
64- sh.verlet.step(sh.s, h);
65- const a = s, b = sh.s, ur = this.uRef, pr = this.pRef;
66- let d2 = 0;
67- for (let j = 0; j < a.u.length; j++) {
68- d2 += ((b.u[j] - a.u[j]) / ur) ** 2 + ((b.p[j] - a.p[j]) / pr) ** 2;
69- }
70- const d = Math.sqrt(d2);
71- const r = D0 / d;
72- for (let j = 0; j < a.u.length; j++) {
73- b.u[j] = a.u[j] + (b.u[j] - a.u[j]) * r;
74- b.p[j] = a.p[j] + (b.p[j] - a.p[j]) * r;
75- }
76- sh.verlet.reset(b); // forces at the moved displacements
77-
78- const g = Math.log(d / D0);
79- this.t = t;
80- this.totalSum += g;
81- this.windowSum += g;
82- this.ts[this.head] = t;
83- this.gs[this.head] = g;
84- this.head = (this.head + 1) % CAP;
85- while (this.tail !== this.head && this.ts[this.tail] <= t - this.window) {
86- this.windowSum -= this.gs[this.tail];
87- this.tail = (this.tail + 1) % CAP;
88- }
89- }
90-
91- // Exponent over the last `window` seconds (or since the start, if shorter), 1/s.
92- windowed(): number {
93- const span = Math.min(this.window, this.t - this.t0);
94- return this.shadow && span > 0 ? this.windowSum / span : NaN;
95- }
96-
97- // Exponent since the start, 1/s.
98- mean(): number {
99- const span = this.t - this.t0;
100- return this.shadow && span > 0 ? this.totalSum / span : NaN;
11+// A shadow of the chain's state s, for Lyapunov.start. s is read at every
12+// step, so it must stay the chain's state object until the next start.
13+export function shadowOf(P: Params, s: State): Shadow {
14+ const ur = d0(P), pr = Math.sqrt(P.m * P.k) * ur;
15+ const sh = copyState(s);
16+ // Initial separation: equal parts in all 2N coordinates.
17+ const e = D0 / Math.sqrt(2 * P.N);
18+ for (let j = 0; j < P.N; j++) {
19+ sh.u[j] += e * ur;
20+ sh.p[j] += e * pr;
10121 }
22+ const verlet = new Verlet(P, sh);
23+
24+ return {
25+ step: (h) => verlet.step(sh, h),
26+ separation: () => {
27+ let d2 = 0;
28+ for (let j = 0; j < s.u.length; j++) {
29+ d2 += ((sh.u[j] - s.u[j]) / ur) ** 2 + ((sh.p[j] - s.p[j]) / pr) ** 2;
30+ }
31+ return Math.sqrt(d2);
32+ },
33+ pullBack: (r) => {
34+ for (let j = 0; j < s.u.length; j++) {
35+ sh.u[j] = s.u[j] + (sh.u[j] - s.u[j]) * r;
36+ sh.p[j] = s.p[j] + (sh.p[j] - s.p[j]) * r;
37+ }
38+ verlet.reset(sh); // forces at the moved displacements
39+ },
40+ };
10241 }
src/fput-chain/main.tsmodified+4−3View file
@@ -5,7 +5,8 @@ import {
55 } from "./chain";
66 import { newCursor } from "../common/cursor";
77 import { setupHelp } from "../common/help";
8-import { Lyapunov } from "./lyapunov";
8+import { Lyapunov } from "../common/lyapunov";
9+import { shadowOf } from "./lyapunov";
910 import { LiveModeHistory, ModeHistory } from "./modehistory";
1011 import { PoincarePanel } from "./poincare";
1112 import { sectionDef } from "./poincare-core";
@@ -112,7 +113,7 @@ function strikeFrom(state: State, label: string) {
112113 t = 0;
113114 ampl = 0;
114115 resetGhost();
115- lyap.start(P, s, t);
116+ lyap.start(shadowOf(P, s), t);
116117 replayButton.disabled = false;
117118 renderNote();
118119 traces.clear();
@@ -754,7 +755,7 @@ function frame(now: number) {
754755 poincare?.addLive(x + f * (s.u[sec.plot] - x), q + f * (s.p[sec.plot] - q), energy(P, s), inside);
755756 }
756757 ghost?.verlet.step(ghost.s, H);
757- lyap.step(s, H, t + H);
758+ lyap.step(H, t + H);
758759 carry -= H;
759760 t += H;
760761 sinceSample += H;
src/fput-chain/poincare-core.tsmodified+4−3View file
@@ -22,7 +22,8 @@
2222 // when mult q^2 <= 2 m (E - vFactor v(x)) (the outline); c is the positive root.
2323
2424 import { Verlet, energy, stepFor, type Params, type State } from "./chain";
25-import { Lyapunov } from "./lyapunov";
25+import { Lyapunov } from "../common/lyapunov";
26+import { shadowOf } from "./lyapunov";
2627
2728 export const CROSSINGS = 300; // crossings recorded per orbit
2829
@@ -99,7 +100,7 @@ export class SectionOrbit {
99100 this.s = sectionState(this.P, this.D, E, x, q)!;
100101 this.verlet = new Verlet(this.P, this.s);
101102 this.h = stepFor(this.P, energy(this.P, this.s));
102- this.lyap.start(this.P, this.s, 0);
103+ this.lyap.start(shadowOf(this.P, this.s), 0);
103104 this.maxTime = maxTime;
104105 this.points = [x, q];
105106 }
@@ -119,7 +120,7 @@ export class SectionOrbit {
119120 const c = s.u[cross], x = s.u[plot], q = s.p[plot];
120121 this.verlet.step(s, h);
121122 this.t += h;
122- this.lyap.step(s, h, this.t);
123+ this.lyap.step(h, this.t);
123124 if (c < 0 && s.u[cross] >= 0) {
124125 const f = -c / (s.u[cross] - c);
125126 this.points.push(x + f * (s.u[plot] - x), q + f * (s.p[plot] - q));
src/fput-chain/render-worker.tsmodified+4−3View file
@@ -10,7 +10,8 @@
1010 // one after UNDAMPED_LENGTH seconds.
1111
1212 import { Modes, Verlet, copyState, energy, stepFor, type Params, type State } from "./chain";
13-import { Lyapunov } from "./lyapunov";
13+import { Lyapunov } from "../common/lyapunov";
14+import { shadowOf } from "./lyapunov";
1415 import { sectionDef } from "./poincare-core";
1516 import { signalByKey } from "./signals";
1617
@@ -59,7 +60,7 @@ self.onmessage = (e: MessageEvent<RenderRequest>) => {
5960 const modeE = new Float32Array(en.length * r.P.N);
6061 const lyap = new Lyapunov(LYAP_WINDOW);
6162 let t = 0; // physical time (s)
62- if (r.lyapunov) lyap.start(r.P, s, t);
63+ if (r.lyapunov) lyap.start(shadowOf(r.P, s), t);
6364 const sec = sectionDef(r.P.N);
6465 const crossings: number[] = [];
6566 let n = 0;
@@ -79,7 +80,7 @@ self.onmessage = (e: MessageEvent<RenderRequest>) => {
7980 const c = sec ? s.u[sec.cross] : 0, x = sec ? s.u[sec.plot] : 0, q = sec ? s.p[sec.plot] : 0;
8081 verlet.step(s, h);
8182 t += h;
82- if (r.lyapunov) lyap.step(s, h, t);
83+ if (r.lyapunov) lyap.step(h, t);
8384 if (sec && c < 0 && s.u[sec.cross] >= 0 && (r.P.N !== 3 || (s.u[0] === s.u[2] && s.p[0] === s.p[2]))) {
8485 const f = -c / (s.u[sec.cross] - c);
8586 crossings.push(n, x + f * (s.u[sec.plot] - x), q + f * (s.p[sec.plot] - q));