Bending stiffness κ (pinned ends) as a parameter; aliasing warning; note Lyapunov shadow only when shown; settings to try
9 changed files+111−41
src/fput-chain/README.mdmodified+12−2View file
@@ -2,14 +2,23 @@
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
7-"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). 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.
7+"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
99 Poincaré section (shown only for N = 2 and N = 3, below the sped-up panel, open; adapted from the pendulum app): its own energy slider (log, 10 µJ to 100 J, default 2.51 J), same parameters, no damping. N = 2: each time mass 1 passes its rest position moving right, mass 2's (u2, p2) is plotted. N = 3: only the symmetric family (u1 = u3, p1 = p3; only modes 1 and 3 move, since the cubic force couples odd modes only to odd modes), which is two-dimensional: each time mass 2 passes its rest position moving right, mass 1's (u1, p1) is plotted. On either section the springs are stretched by 0 or +-x, so V = vFactor (k x^2 / 2 + k3 x^4 / 4) (vFactor 2 or 4), the plotted momentum counts once (N = 2) or twice (N = 3) in the kinetic energy, the outline is where the crossing mass's momentum vanishes, and that momentum follows from E. 42 orbits seeded on a grid, 300 crossings each, in a worker (`src/poincare-core.ts`, `src/poincare-worker.ts`, `src/poincare.ts`; integrator step stepFor(E), energy error < 5e-5; ~1 s per section offline, a few s in the pane), drawn as they arrive; per-orbit colors or Color by λ (each orbit's Lyapunov exponent over its run, log scale); hover to pick out an orbit with its λ; click inside the outline to add an orbit there, and with Click to strike (on) also strike the chain from that exact state (and play it, if the sound panel is open). Follow (as in the pendulum app): "real time" makes the section's energy track the chain's as it decays (recomputed after a 3% change, at most every 2 s, and at each strike); "sped up" (which opens the sped-up panel) shows the playing note's crossings (found in the render worker, for N = 3 only while exactly symmetric, and drawn as the audio clock reaches them) as a fading trail of the last 150, with the axes fitted to the note's current energy (refit every 2%) and each dot fitted from its own energy; nothing is computed live. "Precompute sections" (shown only with Follow "sped up" until complete) builds an atlas (`src/atlas.ts`): sections at 32 energies log-spaced from 1 mJ to 100 J (a factor ~1.45 apart), 42 orbits x 150 crossings, highest first, in its own worker (~22 s in the pane for N = 3), stored in IndexedDB per N, m, k, k3 and deleted after 7 days; the sound view shows the atlas section nearest the current energy, fitted to its allowed region (exact in the linear regime, below ~25 mJ, where sections are scaled copies; an approximation across the gaps above). Without an atlas only the outline is drawn behind the trail. Setting the section's energy by hand pauses everything and clears the rings and the trail; Follow takes over again on Resume. The chain's own crossings show as dark rings when within 5% of the section's energy and, for N = 3, only while the chain is exactly symmetric (else the legend says "chain not symmetric, not shown"; strikes on mass 2 or on mode 1 or 3 are symmetric, mode shapes are made exactly mirror-symmetric for this). Findings: the two-mass chain appears regular at every energy tried (10 mJ to 100 J): nested closed curves and island chains, no chaotic sea; the grid orbits' λ falls as the run lengthens (at 10 J, median 0.045 /s over 300 crossings, 0.007 /s over 3000), as for regular orbits. A likely reason: in mode coordinates (u = S Q) the potential is (k/2) rho^2 + (k3/8) rho^4 with rho^2 = Q1^2 + 3 Q2^2, a function of one combination of the two modes. The N = 3 symmetric family is chaotic: grid orbits with λ > 0.1 /s: 0 of 42 at 25 mJ, 10 at 0.25 J (islands in a chaotic sea), 37 at 2.51 J. It is unstable to breaking the symmetry: an exactly symmetric start stays symmetric (the arithmetic is exactly mirror-symmetric), but an asymmetry of 1e-12 grows to full size within ~100 s at 0.25 to 25 J. The ghost and the Lyapunov shadow offset all masses equally, so after a symmetric strike they too stay symmetric and measure the exponent within the family.
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+
1322 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).
1423
1524 ```bash
@@ -20,6 +29,7 @@ npm run dev
2029 ## Model
2130
2231 - 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.
2333 - 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.
2434 - 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.
2535 - 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.N)!, E, GRID[0], GRID[1]);
117+ const seeds = seedGrid(P, sectionDef(P)!, 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+47−13View file
@@ -3,6 +3,13 @@
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.
613 // Optional damping: linear drag on each mass proportional to its mass,
714 // force -gamma m v = -gamma p, so every linear mode's amplitude decays as
815 // exp(-t / tau) with tau = 2 / gamma.
@@ -15,6 +22,7 @@ export interface Params {
1522 m: number; // kg, each mass
1623 k: number; // N/m, linear spring constant
1724 k3: number; // N/m^3, cubic spring coefficient
25+ kappa: number; // N/m, bending stiffness (0 = none)
1826 a: number; // m, rest spacing (for drawing only; it does not enter the equations)
1927 tau: number; // s, amplitude decay time from drag (all modes); 0 = no damping
2028 }
@@ -22,18 +30,24 @@ export interface Params {
2230 // Drag rate gamma (1/s).
2331 export const dragRate = (P: Params) => (P.tau > 0 ? 2 / P.tau : 0);
2432
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.
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.
3346 export function stepFor(P: Params, E: number): number {
3447 const d = Math.sqrt(Math.sqrt((4 * Math.max(E, 0)) / P.k3));
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));
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));
3751 }
3852
3953 export interface State {
@@ -64,6 +78,18 @@ export function forces(P: Params, u: Float64Array, F: Float64Array): void {
6478 if (j > 0) F[j - 1] = f - left;
6579 left = f;
6680 }
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+ }
6793 }
6894
6995 export function kinetic(P: Params, s: State): number {
@@ -78,6 +104,13 @@ export function potential(P: Params, s: State): number {
78104 const d2 = stretch(s.u, j) ** 2;
79105 V += (P.k * d2) / 2 + (P.k3 * d2 * d2) / 4;
80106 }
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+ }
81114 return V;
82115 }
83116
@@ -131,7 +164,8 @@ export class Verlet {
131164
132165 // Normal modes of the linear chain (k3 = 0), mode n = 1..N stored at index
133166 // n - 1: shape S[n][j] = sqrt(2 / (N + 1)) sin(pi n (j + 1) / (N + 1)),
134-// angular frequency omega_n = 2 sqrt(k / m) sin(pi n / (2 (N + 1))) (rad/s).
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))).
135169 // S is symmetric and its own inverse, so mode amplitudes are Q = S u, P = S p.
136170 export class Modes {
137171 readonly N: number;
@@ -142,12 +176,12 @@ export class Modes {
142176 const { N, m } = P;
143177 this.N = N;
144178 this.m = m;
145- const w0 = Math.sqrt(P.k / m);
146179 this.omega = new Float64Array(N);
147180 this.S = [];
148181 const c = Math.sqrt(2 / (N + 1));
149182 for (let n = 1; n <= N; n++) {
150- this.omega[n - 1] = 2 * w0 * Math.sin((Math.PI * n) / (2 * (N + 1)));
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);
151185 const row = new Float64Array(N);
152186 for (let j = 0; j < N; j++) row[j] = c * Math.sin((Math.PI * n * (j + 1)) / (N + 1));
153187 // Make the mirror symmetry exact (row[N-1-j] = (-1)^(n+1) row[j]; sin
src/fput-chain/index.htmlmodified+2−0View file
@@ -174,6 +174,7 @@
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>
177178 <label>a <input id="p-a" type="number" step="any" /> <span class="unit">m (spacing)</span></label>
178179 <label>τ <input id="p-tau" type="number" step="any" /> <span class="unit">s (decay, 0 = off)</span></label>
179180 </div>
@@ -261,6 +262,7 @@
261262 <template id="help-params">
262263 <h2>Parameters</h2>
263264 <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>
264266 </template>
265267 <template id="help-realtime">
266268 <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(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
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
1212 // depend on this choice; finite-time values depend on it a little.
1313
14-import { Verlet, copyState, d0, type Params, type State } from "./chain";
14+import { Verlet, copyState, d0, kLin, 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 * P.k) * this.uRef;
41+ this.pRef = Math.sqrt(P.m * kLin(P)) * 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+27−8View file
@@ -1,5 +1,5 @@
11 import {
2- Modes, Verlet, copyState, d0, energy, kinetic, newState, stepFor, potential, strikeMass, strikeMode, type Params, type State,
2+ Modes, Verlet, copyState, d0, energy, kLin, kinetic, newState, stepFor, potential, strikeMass, strikeMode, type Params, type State,
33 } from "./chain";
44 import { newCursor } from "./cursor";
55 import { setupHelp } from "./help";
@@ -33,7 +33,8 @@ 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-let P: Params = { N: 2, m: 1, k: 10, k3: 4000, a: 0.5, tau: 540 };
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 };
3738 const SAMPLE_DT = 0.01; // s between plotted samples
3839
3940 // The chain. It sits at rest until the first Strike; changing parameters
@@ -268,6 +269,7 @@ const paramInputs = {
268269 m: document.getElementById("p-m") as HTMLInputElement,
269270 k: document.getElementById("p-k") as HTMLInputElement,
270271 k3: document.getElementById("p-k3") as HTMLInputElement,
272+ kappa: document.getElementById("p-kappa") as HTMLInputElement,
271273 a: document.getElementById("p-a") as HTMLInputElement,
272274 tau: document.getElementById("p-tau") as HTMLInputElement,
273275 };
@@ -282,8 +284,12 @@ function showDerived() {
282284 const list = P.N <= 8 ? hz.map((f) => f.toFixed(3)).join(", ") : `${hz[0].toFixed(3)} … ${hz[P.N - 1].toFixed(3)}`;
283285 document.getElementById("params")!.textContent =
284286 `normal modes ${list} Hz\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` +
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`) +
287293 `decay time ${damped() ? `${P.tau} s (amplitude, all modes)` : "none (no damping)"}\n` +
288294 `integrator step ${(1000 * H).toPrecision(3)} ms (from the strike energy; ` +
289295 `${(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)`;
@@ -325,8 +331,13 @@ const predictedNoteLength = () =>
325331 const fmtDuration = (x: number) => (x >= 120 ? `${(x / 60).toFixed(0)} min` : `${x.toFixed(x >= 10 ? 0 : 1)} s`);
326332 function showPitch() {
327333 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;
328338 document.getElementById("pitch-value")!.textContent =
329- `lowest mode ${f[0].toFixed(0)} Hz · highest ${f[f.length - 1].toFixed(0)} Hz`;
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)` : "");
330341 // The sound runs on its own clock: the same strike, sped up.
331342 const S = speedup(), note = predictedNoteLength();
332343 soundClock.textContent = damped()
@@ -458,6 +469,9 @@ function renderNote() {
458469 const req: RenderRequest = {
459470 id: ++requestId, P, state: lastStrike, speedup: speedup(), sampleRate: ctx.sampleRate,
460471 maxLength: 2 * predictedNoteLength(), signal: sig.key, j: j - 1,
472+ // The note's Lyapunov shadow doubles the computation, so it runs only
473+ // when a display needs it (Lyapunov chart or λ vs E).
474+ lyapunov: soundTraces.shownLabels().includes("Lyapunov") || !scatterFig.hidden,
461475 };
462476 noteInfo = `${sig.label}${sig.perMass ? ` ${j}` : ""}`;
463477 noteFile =
@@ -558,8 +572,13 @@ nInput.addEventListener("change", () => {
558572 });
559573 document.getElementById("param-apply")!.addEventListener("click", () => {
560574 const v = Object.fromEntries(Object.entries(paramInputs).map(([key, el]) => [key, Number(el.value)]));
561- if (!Object.entries(v).every(([key, x]) => Number.isFinite(x) && (x > 0 || (key === "tau" && x === 0)))) {
562- paramStatus.textContent = "parameters must be positive (tau may be 0)";
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";
563582 return;
564583 }
565584 paramStatus.textContent = "";
@@ -735,7 +754,7 @@ function frame(now: number) {
735754 const dt = Math.min(0.1, (now - last) / 1000);
736755 last = now;
737756 if (running) {
738- const sec = sectionDef(P.N);
757+ const sec = sectionDef(P);
739758 carry += Number(speedSelect.value) * dt;
740759 while (carry >= H) {
741760 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+6−3View file
@@ -35,8 +35,11 @@ export interface SectionDef {
3535 pLabel: string;
3636 }
3737
38-// The section for N masses, or null if there is none.
39-export function sectionDef(N: number): SectionDef | null {
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;
4043 if (N === 2) return { cross: 0, plot: 1, vFactor: 2, mult: 1, uLabel: "u2", pLabel: "p2" };
4144 if (N === 3) return { cross: 1, plot: 0, vFactor: 4, mult: 2, uLabel: "u1 = u3", pLabel: "p1 = p3" };
4245 return null;
@@ -95,7 +98,7 @@ export class SectionOrbit {
9598
9699 constructor(P: Params, E: number, x: number, q: number, maxTime: number) {
97100 this.P = { ...P, tau: 0 };
98- this.D = sectionDef(P.N)!;
101+ this.D = sectionDef(P)!;
99102 this.s = sectionState(this.P, this.D, E, x, q)!;
100103 this.verlet = new Verlet(this.P, this.s);
101104 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.N);
157+ this.D = sectionDef(P);
158158 this.live = [];
159159 this.noteDots = [];
160160 this.outside = false;
src/fput-chain/render-worker.tsmodified+10−8View file
@@ -1,10 +1,11 @@
11 // Renders a note: the chain from a given state, sped up, sampled as audio.
22 // Each audio sample advances the physics by speedup / sampleRate seconds, in
33 // equal substeps no longer than stepFor the current energy (re-chosen every
4-// DIAG_EVERY samples, so the steps lengthen as a damped note decays). A
5-// Lyapunov shadow runs along (as in the animation: 10 s window of physical
6-// time); the energy and the windowed exponent are sent with the samples, one
7-// value each per DIAG_EVERY samples. A damped note ends
4+// DIAG_EVERY samples, so the steps lengthen as a damped note decays). On
5+// request a Lyapunov shadow runs along (as in the animation: 10 s window of
6+// physical time; it doubles the cost). The energy, the windowed exponent
7+// (NaN without the shadow) and the mode energies are sent with the samples,
8+// one value (or set) per DIAG_EVERY samples. A damped note ends
89 // once its energy has fallen to END_RATIO of the start (-60 dB); an undamped
910 // one after UNDAMPED_LENGTH seconds.
1011
@@ -27,6 +28,7 @@ export interface RenderRequest {
2728 maxLength: number; // s of sound, at most
2829 signal: string; // signal key
2930 j: number; // mass index (0-based) for per-mass signals
31+ lyapunov: boolean; // run the Lyapunov shadow
3032 }
3133
3234 export interface RenderResult {
@@ -57,15 +59,15 @@ self.onmessage = (e: MessageEvent<RenderRequest>) => {
5759 const modeE = new Float32Array(en.length * r.P.N);
5860 const lyap = new Lyapunov(LYAP_WINDOW);
5961 let t = 0; // physical time (s)
60- lyap.start(r.P, s, t);
61- const sec = sectionDef(r.P.N);
62+ if (r.lyapunov) lyap.start(r.P, s, t);
63+ const sec = sectionDef(r.P);
6264 const crossings: number[] = [];
6365 let n = 0;
6466 for (; n < nMax; n++) {
6567 if (n % DIAG_EVERY === 0) {
6668 const E = energy(r.P, s);
6769 en[n / DIAG_EVERY] = E;
68- lam[n / DIAG_EVERY] = lyap.windowed();
70+ lam[n / DIAG_EVERY] = r.lyapunov ? lyap.windowed() : NaN;
6971 modes.energies(s, Ek);
7072 modeE.set(Ek, (n / DIAG_EVERY) * r.P.N);
7173 if (r.P.tau > 0 && E <= END_RATIO * E0) break;
@@ -77,7 +79,7 @@ self.onmessage = (e: MessageEvent<RenderRequest>) => {
7779 const c = sec ? s.u[sec.cross] : 0, x = sec ? s.u[sec.plot] : 0, q = sec ? s.p[sec.plot] : 0;
7880 verlet.step(s, h);
7981 t += h;
80- lyap.step(s, h, t);
82+ if (r.lyapunov) lyap.step(s, h, t);
8183 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]))) {
8284 const f = -c / (s.u[sec.cross] - c);
8385 crossings.push(n, x + f * (s.u[sec.plot] - x), q + f * (s.p[sec.plot] - q));