Remove bending stiffness κ (moving to a separate app); keep note Lyapunov shadow only when shown
9 changed files+37−98
src/fput-chain/README.mdmodified+5−11View 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 (`src/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 `src/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 `src/traces.ts` and `src/cursor.ts` of double-pendulum-chaos, 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 (`src/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 `src/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 `src/traces.ts` and `src/cursor.ts` of double-pendulum-chaos, 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 (`src/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 (`src/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 (`src/sound.ts`). Save WAV saves the last note as heard (16-bit mono PCM, peak at 0.9; `src/wav.ts`; named e.g. fput-N4-mass1-2.51J-v1-100Hz.wav). Chips (as in the pendulum app, plus Modes): Spectrogram (on; pendulum's `NoteSpectrogram`: 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 `src/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 (`src/scatter.ts` from the pendulum app: 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
@@ -10,15 +10,6 @@ Poincaré section (shown only for N = 2 and N = 3, below the sped-up panel, open
1010
1111 Each card and panel (chain, Parameters, Real time, Sped up · listen, Poincaré section) has a small "?" button that opens a short help dialog (`src/help.ts` from the pendulum app; the texts are `<template>`s in `index.html`).
1212
13-## Settings to try
14-
15-- Cascade: Masses 32, strike mode 1 at 115 mJ, Real time Mode history on (energy climbs the modes over ~100 s); sped-up spectrogram at pitch 30 Hz.
16-- Cymbal-like beam: Parameters k = 0, κ = 20; Masses 32; strike mass 8 at 5 J; Signal velocity of mass 8; pitch ~12 Hz (a lower pitch lengthens the note without more computing; computing time is set by τ).
17-- Same beam struck gently (0.05 J): a clean inharmonic chord without the wash; compare λ (Lyapunov chip) with the 5 J strike.
18-- Stiff string: k = 10, κ = 2, Masses 32, mode 1 at 115 mJ, pitch 30 Hz, Mode history on; compare with κ = 0 (does stretching the partials slow the cascade?).
19-- Quicker beam: as the cymbal-like beam with Masses 16.
20-- Watch the pitch readout: "aliased" means the top mode is above half the sample rate; lower the pitch.
21-
2213 Live site: https://fput-chain-chaos--listening-for-chaos.vault1-sites.magland.org/ (built and deployed by `.github/workflows/pages.yml` on every push to main).
2314
2415 ```bash
@@ -26,10 +17,13 @@ npm install
2617 npm run dev
2718 ```
2819
20+## Settings to try
21+
22+- Cascade: Masses 32, strike mode 1 at 115 mJ, Real time Mode history on (energy climbs the modes over ~100 s); sped-up spectrogram at pitch 30 Hz.
23+
2924 ## Model
3025
3126 - SI units. Defaults m = 1 kg, k = 10 N/m, k3 = 4000 N/m^3, a = 0.5 m. The cubic tension equals the linear one at d0 = sqrt(k / k3) = 5 cm; how nonlinear the motion is depends only on E / (k d0^2) (25 mJ by default). The rest spacing only enters the drawing.
32-- Bending (optional, κ, default 0): energy κ c_j^2 / 2 on the curvature c_j = u_{j+1} - 2 u_j + u_{j-1} at every mass, pinned ends (no curvature at the walls); force -κ (c_{j-1} - 2 c_j + c_{j+1}). Physically this belongs to masses moving sideways (a beam; the springs then stand for tension and stretching), so with κ > 0 read the sideways view as the motion. Pinned ends keep the sine mode shapes; the frequencies become ω_n^2 = (4 k s^2 + 16 κ s^4) / m, s = sin(π n / (2 (N + 1))): k = 10, κ = 5 gives stretched partials 1, 2.01, 3.04, 4.11, 5.21; k = 0, κ = 20 gives ~n^2 (1, 3.99, 8.95, 15.8, 24.6), cymbal-like (k may be 0 when κ > 0). The reference length for the Lyapunov distance and the nonlinear scale use max(k, κ). The Poincaré section is hidden with bending (its formulas leave it out). Bending spreads the modes widely (top/bottom ~1500 at N = 32, k = 0), so the pitch readout warns when the highest mode lies above half the sample rate (aliased); at 30 Hz the top is 13 kHz. Checked: sine modes are eigenvectors of the linear forces (1e-13), forces against finite differences (4e-10), energy error <= 4e-5 from 10 mJ to 100 J for (k, κ) = (10, 0), (10, 5), (0, 20); κ = 0 gives the same steps and errors as before.
3327 - Fixed ends. Normal modes of the linear chain: shapes sqrt(2 / (N + 1)) sin(pi n j / (N + 1)), angular frequencies omega_n = 2 sqrt(k / m) sin(pi n / (2 (N + 1))); 0.311, 0.592, 0.814, 0.957 Hz for the defaults. Mode energy P_n^2 / 2m + m omega_n^2 Q_n^2 / 2.
3428 - Damping: linear drag on each mass proportional to its mass, force -gamma m v = -gamma p, so every linear mode's amplitude decays as exp(-t / tau), tau = 2 / gamma (default 540 s: ~10% amplitude per minute in the animation; 0 = off). Applied exactly in half steps around each Verlet step (p scaled by exp(-gamma h / 2)), with the kinetic energy removed added to W. Checked: modes 1 and 4 at 1e-8 J decay to 0.13539 and 0.13536 of their energy after tau (e^-2 = 0.13534); energy balance error 3e-5 over the default note (6000 s); -60 dB reached at 3635 s vs 3730 s predicted.
3529 - Integrator: velocity Verlet (order 2, symplectic; H is separable, so an explicit method suffices). Step from the energy (`stepFor`): min(0.005 sqrt(m / k), 0.02 sqrt(m / k_eff)), k_eff = k + 3 k3 d^2 with d = (4 E / k3)^(1/4) the largest stretch a spring holding all of E could reach. Set at each strike in the animation (1.58 ms when weak, 0.81 ms at 2.5 J, 0.32 ms at 100 J); re-chosen every 256 samples in the sound, so steps lengthen as a note decays. Relative energy error 1e-5 to 5e-5 from 1 mJ to 100 J (N = 4 and 32, 300 s). Render time in the pane: 0.2 to 0.3 s (N = 4), 0.6 to 0.8 s (N = 32), 2.5 to 100 J. At 100 J stretches reach ~0.55 m, more than the drawing's rest spacing, so masses pass each other in the to-scale drawing (the model has no contact).
src/fput-chain/atlas.tsmodified+1−1View file
@@ -114,7 +114,7 @@ export class Atlas {
114114 private section(P: Params, E: number, job: number): Promise<AtlasEntry | null> {
115115 if (!this.worker) this.worker = new Worker(new URL("./poincare-worker.ts", import.meta.url), { type: "module" });
116116 const worker = this.worker;
117- const seeds = seedGrid(P, sectionDef(P)!, E, GRID[0], GRID[1]);
117+ const seeds = seedGrid(P, sectionDef(P.N)!, E, GRID[0], GRID[1]);
118118 const points: number[][] = seeds.map(() => []);
119119 const lambdas: number[] = seeds.map(() => NaN);
120120 return new Promise((resolve) => {
src/fput-chain/chain.tsmodified+13−47View file
@@ -3,13 +3,6 @@
33 // SI units. The cubic term matches the linear one at the stretch
44 // d0 = sqrt(k / k3); how nonlinear the motion is depends only on the energy
55 // relative to k d0^2.
6-// Optional bending stiffness kappa (0 by default): energy kappa c_j^2 / 2 on
7-// the curvature c_j = u_{j+1} - 2 u_j + u_{j-1} at every mass, with pinned
8-// ends (no curvature at the walls). Bending only makes physical sense for
9-// masses moving sideways (a beam), so with kappa > 0 the chain is best read
10-// that way, the springs then standing for tension and stretching. Pinned
11-// ends keep the sine mode shapes; bending raises the mode frequencies
12-// unevenly (partials ~ n^2 when tension is small), making them inharmonic.
136 // Optional damping: linear drag on each mass proportional to its mass,
147 // force -gamma m v = -gamma p, so every linear mode's amplitude decays as
158 // exp(-t / tau) with tau = 2 / gamma.
@@ -22,7 +15,6 @@ export interface Params {
2215 m: number; // kg, each mass
2316 k: number; // N/m, linear spring constant
2417 k3: number; // N/m^3, cubic spring coefficient
25- kappa: number; // N/m, bending stiffness (0 = none)
2618 a: number; // m, rest spacing (for drawing only; it does not enter the equations)
2719 tau: number; // s, amplitude decay time from drag (all modes); 0 = no damping
2820 }
@@ -30,24 +22,18 @@ export interface Params {
3022 // Drag rate gamma (1/s).
3123 export const dragRate = (P: Params) => (P.tau > 0 ? 2 / P.tau : 0);
3224
33-// The linear stiffness that sets the nonlinear scale: tension, or bending if
34-// larger (so it is not 0 for a beam without tension). The cubic term matches
35-// it at d0; d0 and sqrt(m kLin) d0 are the reference length and momentum.
36-export const kLin = (P: Params) => Math.max(P.k, P.kappa);
37-export const d0 = (P: Params) => Math.sqrt(kLin(P) / P.k3);
38-
39-// Integrator step (s) for motion with energy E. Linear part: 0.01 / omega_max,
40-// with omega_max = sqrt((4 k + 16 kappa) / m) the fastest linear mode (1/200
41-// of sqrt(m / k) without bending, ~600 steps per period). Strong motion
42-// stiffens the springs: a spring holding all of E stretches by at most
43-// d = (4 E / k3)^(1/4), where its stiffness is k + 3 k3 d^2; the step is kept
44-// at 1/50 of sqrt(m) over the root of that (plus 4 kappa with bending).
45-// Relative energy error stays at a few 1e-5 from 1 mJ to 100 J.
25+export const d0 = (P: Params) => Math.sqrt(P.k / P.k3);
26+
27+// Integrator step (s) for motion with energy E. Linear part: 1/200 of
28+// sqrt(m / k), ~600 steps per period of the fastest linear mode. Strong
29+// motion stiffens the springs: a spring holding all of E stretches by at
30+// most d = (4 E / k3)^(1/4), where its stiffness is k + 3 k3 d^2; the step
31+// is kept at 1/50 of sqrt(m) over the root of that. Relative energy error
32+// stays at a few 1e-5 from 1 mJ to 100 J.
4633 export function stepFor(P: Params, E: number): number {
4734 const d = Math.sqrt(Math.sqrt((4 * Math.max(E, 0)) / P.k3));
48- const stiff = P.k + 3 * P.k3 * d * d + 4 * P.kappa;
49- const wMax = Math.sqrt((4 * P.k + 16 * P.kappa) / P.m);
50- return Math.min(0.01 / wMax, 0.02 * Math.sqrt(P.m / stiff));
35+ const stiff = P.k + 3 * P.k3 * d * d;
36+ return Math.min(0.005 * Math.sqrt(P.m / P.k), 0.02 * Math.sqrt(P.m / stiff));
5137 }
5238
5339 export interface State {
@@ -78,18 +64,6 @@ export function forces(P: Params, u: Float64Array, F: Float64Array): void {
7864 if (j > 0) F[j - 1] = f - left;
7965 left = f;
8066 }
81- if (P.kappa > 0) {
82- // Bending: F_i = -kappa (c_{i-1} - 2 c_i + c_{i+1}), curvature zero at
83- // the walls (pinned).
84- const c = (j: number) => (j < 0 || j >= N ? 0 : (j + 1 < N ? u[j + 1] : 0) - 2 * u[j] + (j > 0 ? u[j - 1] : 0));
85- let cPrev = 0, cHere = c(0);
86- for (let i = 0; i < N; i++) {
87- const cNext = c(i + 1);
88- F[i] -= P.kappa * (cPrev - 2 * cHere + cNext);
89- cPrev = cHere;
90- cHere = cNext;
91- }
92- }
9367 }
9468
9569 export function kinetic(P: Params, s: State): number {
@@ -104,13 +78,6 @@ export function potential(P: Params, s: State): number {
10478 const d2 = stretch(s.u, j) ** 2;
10579 V += (P.k * d2) / 2 + (P.k3 * d2 * d2) / 4;
10680 }
107- if (P.kappa > 0) {
108- const u = s.u, N = P.N;
109- for (let j = 0; j < N; j++) {
110- const c = (j + 1 < N ? u[j + 1] : 0) - 2 * u[j] + (j > 0 ? u[j - 1] : 0);
111- V += 0.5 * P.kappa * c * c;
112- }
113- }
11481 return V;
11582 }
11683
@@ -164,8 +131,7 @@ export class Verlet {
164131
165132 // Normal modes of the linear chain (k3 = 0), mode n = 1..N stored at index
166133 // n - 1: shape S[n][j] = sqrt(2 / (N + 1)) sin(pi n (j + 1) / (N + 1)),
167-// angular frequency omega_n = 2 sqrt(k / m) sin(pi n / (2 (N + 1))) (rad/s);
168-// with bending, omega_n^2 = (4 k s^2 + 16 kappa s^4) / m, s = sin(pi n / (2 (N + 1))).
134+// angular frequency omega_n = 2 sqrt(k / m) sin(pi n / (2 (N + 1))) (rad/s).
169135 // S is symmetric and its own inverse, so mode amplitudes are Q = S u, P = S p.
170136 export class Modes {
171137 readonly N: number;
@@ -176,12 +142,12 @@ export class Modes {
176142 const { N, m } = P;
177143 this.N = N;
178144 this.m = m;
145+ const w0 = Math.sqrt(P.k / m);
179146 this.omega = new Float64Array(N);
180147 this.S = [];
181148 const c = Math.sqrt(2 / (N + 1));
182149 for (let n = 1; n <= N; n++) {
183- const sn = Math.sin((Math.PI * n) / (2 * (N + 1)));
184- this.omega[n - 1] = Math.sqrt((4 * P.k * sn * sn + 16 * P.kappa * sn ** 4) / m);
150+ this.omega[n - 1] = 2 * w0 * Math.sin((Math.PI * n) / (2 * (N + 1)));
185151 const row = new Float64Array(N);
186152 for (let j = 0; j < N; j++) row[j] = c * Math.sin((Math.PI * n * (j + 1)) / (N + 1));
187153 // Make the mirror symmetry exact (row[N-1-j] = (-1)^(n+1) row[j]; sin
src/fput-chain/index.htmlmodified+0−2View file
@@ -174,7 +174,6 @@
174174 <label>m <input id="p-m" type="number" step="any" /> <span class="unit">kg</span></label>
175175 <label>k <input id="p-k" type="number" step="any" /> <span class="unit">N/m</span></label>
176176 <label>k3 <input id="p-k3" type="number" step="any" /> <span class="unit">N/m³</span></label>
177- <label title="Bending stiffness (pinned ends); 0 = none. Makes the partials inharmonic.">κ <input id="p-kappa" type="number" step="any" /> <span class="unit">N/m (bending)</span></label>
178177 <label>a <input id="p-a" type="number" step="any" /> <span class="unit">m (spacing)</span></label>
179178 <label>τ <input id="p-tau" type="number" step="any" /> <span class="unit">s (decay, 0 = off)</span></label>
180179 </div>
@@ -262,7 +261,6 @@
262261 <template id="help-params">
263262 <h2>Parameters</h2>
264263 <p>Mass, spring constants k and k3, the rest spacing a (drawing only), and the decay time τ of the amplitude (drag on the masses; 0 for none). Apply restarts the chain at rest. The cubic term matches the linear one at the stretch d0 listed below; the motion turns nonlinear at energies around k d0².</p>
265- <p>κ adds bending stiffness (pinned ends), which spreads the partials apart unevenly: with k = 10 and small κ they are stretched like a piano's; with k = 0 and κ = 20 they go roughly as n², and the sound turns cymbal-like. Bending belongs to masses moving sideways (a beam), so read the sideways view as the motion. The Poincaré section is not available with bending.</p>
266264 </template>
267265 <template id="help-realtime">
268266 <h2>Real time</h2>
src/fput-chain/lyapunov.tsmodified+5−5View file
@@ -6,12 +6,12 @@
66 // over a sliding window (finite-time exponent, which follows a damped
77 // ring-down) and since the start.
88 //
9-// Distance: displacements in units of d0 = sqrt(kLin / k3) and momenta in
10-// units of sqrt(m kLin) d0 (the momentum of a mass swinging at sqrt(kLin / m)
11-// with amplitude d0; kLin is k, or kappa if larger), so both count comparably. The long-time exponent does not
9+// Distance: displacements in units of d0 = sqrt(k / k3) and momenta in units
10+// 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
1212 // depend on this choice; finite-time values depend on it a little.
1313
14-import { Verlet, copyState, d0, kLin, type Params, type State } from "./chain";
14+import { Verlet, copyState, d0, type Params, type State } from "./chain";
1515
1616 const D0 = 1e-8; // shadow distance
1717 const CAP = 1 << 18; // per-step growth values kept (> window / step)
@@ -38,7 +38,7 @@ export class Lyapunov {
3838 // Start (or restart) measuring from the chain's current state at time t.
3939 start(P: Params, s: State, t: number): void {
4040 this.uRef = d0(P);
41- this.pRef = Math.sqrt(P.m * kLin(P)) * this.uRef;
41+ this.pRef = Math.sqrt(P.m * P.k) * this.uRef;
4242 const sh = copyState(s);
4343 // Initial separation: equal parts in all 2N coordinates.
4444 const e = D0 / Math.sqrt(2 * P.N);
src/fput-chain/main.tsmodified+8−24View file
@@ -1,5 +1,5 @@
11 import {
2- Modes, Verlet, copyState, d0, energy, kLin, kinetic, newState, stepFor, potential, strikeMass, strikeMode, type Params, type State,
2+ Modes, Verlet, copyState, d0, energy, kinetic, newState, stepFor, potential, strikeMass, strikeMode, type Params, type State,
33 } from "./chain";
44 import { newCursor } from "./cursor";
55 import { setupHelp } from "./help";
@@ -33,8 +33,7 @@ const modeStats = document.getElementById("mode-stats")!;
3333 // tau: amplitude decay time from drag; 540 s barely shows in the animation
3434 // (~10% amplitude per minute) but makes a struck note ring down in about
3535 // 12 s of sound. 0 = no damping.
36-// kappa: bending stiffness, 0 by default (see chain.ts).
37-let P: Params = { N: 2, m: 1, k: 10, k3: 4000, kappa: 0, a: 0.5, tau: 540 };
36+let P: Params = { N: 2, m: 1, k: 10, k3: 4000, a: 0.5, tau: 540 };
3837 const SAMPLE_DT = 0.01; // s between plotted samples
3938
4039 // The chain. It sits at rest until the first Strike; changing parameters
@@ -269,7 +268,6 @@ const paramInputs = {
269268 m: document.getElementById("p-m") as HTMLInputElement,
270269 k: document.getElementById("p-k") as HTMLInputElement,
271270 k3: document.getElementById("p-k3") as HTMLInputElement,
272- kappa: document.getElementById("p-kappa") as HTMLInputElement,
273271 a: document.getElementById("p-a") as HTMLInputElement,
274272 tau: document.getElementById("p-tau") as HTMLInputElement,
275273 };
@@ -284,12 +282,8 @@ function showDerived() {
284282 const list = P.N <= 8 ? hz.map((f) => f.toFixed(3)).join(", ") : `${hz[0].toFixed(3)} … ${hz[P.N - 1].toFixed(3)}`;
285283 document.getElementById("params")!.textContent =
286284 `normal modes ${list} Hz\n` +
287- (P.kappa > 0
288- ? `first ratios 1, ${Array.from(modes.omega.slice(1, 6), (w) => (w / modes.omega[0]).toFixed(2)).join(", ")}, … (bending)\n` +
289- `cubic = linear at stretch d0 = sqrt(max(k, κ)/k3) = ${fmtLength(d0(P))}\n` +
290- `nonlinear scale max(k, κ) d0² = ${fmtEnergy(kLin(P) * d0(P) ** 2)}\n`
291- : `cubic = linear at stretch d0 = sqrt(k/k3) = ${fmtLength(d0(P))}\n` +
292- `nonlinear scale k d0² = ${fmtEnergy(P.k * d0(P) ** 2)}\n`) +
285+ `cubic = linear at stretch d0 = sqrt(k/k3) = ${fmtLength(d0(P))}\n` +
286+ `nonlinear scale k d0² = ${fmtEnergy(P.k * d0(P) ** 2)}\n` +
293287 `decay time ${damped() ? `${P.tau} s (amplitude, all modes)` : "none (no damping)"}\n` +
294288 `integrator step ${(1000 * H).toPrecision(3)} ms (from the strike energy; ` +
295289 `${(1000 * stepFor(P, 0)).toPrecision(3)} ms when weak, ${(1000 * stepFor(P, 10 ** Number(energyInput.max))).toPrecision(3)} ms at the top of the slider)`;
@@ -331,13 +325,8 @@ const predictedNoteLength = () =>
331325 const fmtDuration = (x: number) => (x >= 120 ? `${(x / 60).toFixed(0)} min` : `${x.toFixed(x >= 10 ? 0 : 1)} s`);
332326 function showPitch() {
333327 const f = modePitches();
334- // Modes above half the sample rate cannot be represented: they fold back
335- // to lower frequencies (aliasing). Bending spreads the modes widely, so
336- // this happens at ordinary pitches; lower the pitch to avoid it.
337- const nyquist = player.sampleRate / 2;
338328 document.getElementById("pitch-value")!.textContent =
339- `lowest mode ${f[0].toFixed(0)} Hz · highest ${f[f.length - 1].toFixed(0)} Hz` +
340- (f[f.length - 1] > nyquist ? ` (above ${(nyquist / 1000).toFixed(0)} kHz: aliased, lower the pitch)` : "");
329+ `lowest mode ${f[0].toFixed(0)} Hz · highest ${f[f.length - 1].toFixed(0)} Hz`;
341330 // The sound runs on its own clock: the same strike, sped up.
342331 const S = speedup(), note = predictedNoteLength();
343332 soundClock.textContent = damped()
@@ -572,13 +561,8 @@ nInput.addEventListener("change", () => {
572561 });
573562 document.getElementById("param-apply")!.addEventListener("click", () => {
574563 const v = Object.fromEntries(Object.entries(paramInputs).map(([key, el]) => [key, Number(el.value)]));
575- const zeroOk = ["tau", "k", "kappa"];
576- if (!Object.entries(v).every(([key, x]) => Number.isFinite(x) && (x > 0 || (zeroOk.includes(key) && x === 0)))) {
577- paramStatus.textContent = "parameters must be positive (k, κ and τ may be 0)";
578- return;
579- }
580- if (v.k === 0 && v.kappa === 0) {
581- paramStatus.textContent = "k and κ cannot both be 0";
564+ if (!Object.entries(v).every(([key, x]) => Number.isFinite(x) && (x > 0 || (key === "tau" && x === 0)))) {
565+ paramStatus.textContent = "parameters must be positive (tau may be 0)";
582566 return;
583567 }
584568 paramStatus.textContent = "";
@@ -754,7 +738,7 @@ function frame(now: number) {
754738 const dt = Math.min(0.1, (now - last) / 1000);
755739 last = now;
756740 if (running) {
757- const sec = sectionDef(P);
741+ const sec = sectionDef(P.N);
758742 carry += Number(speedSelect.value) * dt;
759743 while (carry >= H) {
760744 const c = sec ? s.u[sec.cross] : 0, x = sec ? s.u[sec.plot] : 0, q = sec ? s.p[sec.plot] : 0;
src/fput-chain/poincare-core.tsmodified+3−6View file
@@ -35,11 +35,8 @@ export interface SectionDef {
3535 pLabel: string;
3636 }
3737
38-// The section for these parameters, or null if there is none (only N = 2
39-// and 3, and only without bending, which the formulas here leave out).
40-export function sectionDef(P: Params): SectionDef | null {
41- const { N } = P;
42- if (P.kappa > 0) return null;
38+// The section for N masses, or null if there is none.
39+export function sectionDef(N: number): SectionDef | null {
4340 if (N === 2) return { cross: 0, plot: 1, vFactor: 2, mult: 1, uLabel: "u2", pLabel: "p2" };
4441 if (N === 3) return { cross: 1, plot: 0, vFactor: 4, mult: 2, uLabel: "u1 = u3", pLabel: "p1 = p3" };
4542 return null;
@@ -98,7 +95,7 @@ export class SectionOrbit {
9895
9996 constructor(P: Params, E: number, x: number, q: number, maxTime: number) {
10097 this.P = { ...P, tau: 0 };
101- this.D = sectionDef(P)!;
98+ this.D = sectionDef(P.N)!;
10299 this.s = sectionState(this.P, this.D, E, x, q)!;
103100 this.verlet = new Verlet(this.P, this.s);
104101 this.h = stepFor(this.P, energy(this.P, this.s));
src/fput-chain/poincare.tsmodified+1−1View file
@@ -154,7 +154,7 @@ export class PoincarePanel {
154154 // New chain parameters (the panel shows only for N = 2 and 3).
155155 setParams(P: Params): void {
156156 this.P = P;
157- this.D = sectionDef(P);
157+ this.D = sectionDef(P.N);
158158 this.live = [];
159159 this.noteDots = [];
160160 this.outside = false;
src/fput-chain/render-worker.tsmodified+1−1View file
@@ -60,7 +60,7 @@ self.onmessage = (e: MessageEvent<RenderRequest>) => {
6060 const lyap = new Lyapunov(LYAP_WINDOW);
6161 let t = 0; // physical time (s)
6262 if (r.lyapunov) lyap.start(r.P, s, t);
63- const sec = sectionDef(r.P);
63+ const sec = sectionDef(r.P.N);
6464 const crossings: number[] = [];
6565 let n = 0;
6666 for (; n < nMax; n++) {