Plate: von Kármán stretching (stress function by sine transforms, exact gradient), thickness h sets it (default 5 cm), step from the strike energy, nonlinear scale
6 changed files+213−26
src/index.htmlmodified+1−1View file
@@ -34,7 +34,7 @@
3434 <li><a href="./double-pendulum/"><strong>Double pendulum</strong><div class="muted">Strike it, watch and hear the motion turn regular as the energy decays; Lyapunov exponent and Poincaré section.</div></a></li>
3535 <li><a href="./fput-chain/"><strong>FPUT chain</strong><div class="muted">A struck chain of N masses with cubic springs; energy spreading up the modes, Lyapunov exponent and Poincaré sections for N = 2 and 3.</div></a></li>
3636 <li><a href="./beam/"><strong>Beam</strong><div class="muted">A struck beam of N masses with bending stiffness and stretching: clamped, pinned or cantilever; inharmonic, cymbal-like overtones.</div></a></li>
37- <li><a href="./plate/"><strong>Plate</strong><div class="muted">A struck square plate of N × N masses with bending stiffness, simply supported; linear for now, modes (m, n) at m² + n².</div></a></li>
37+ <li><a href="./plate/"><strong>Plate</strong><div class="muted">A struck square plate of N × N masses with bending stiffness and von Kármán stretching, simply supported; modes (m, n) at m² + n², pitch glides when struck hard.</div></a></li>
3838 </ul>
3939 </main>
4040 </body>
src/plate/README.mdmodified+6−2View file
@@ -1,9 +1,13 @@
11 # plate
22
3-Struck square plate, step 4 of "Listening for chaos" (after double-pendulum, fput-chain and beam), the first two-dimensional system on the way to cymbal and gong. Work in progress: step 1 (2026-10-09), a linear, simply supported square plate with mallet strike, animation and sped-up note with spectrogram.
3+Struck square plate, step 4 of "Listening for chaos" (after double-pendulum, fput-chain and beam), the first two-dimensional system on the way to cymbal and gong. Work in progress: step 1 (2026-10-09), a linear, simply supported square plate with mallet strike, animation and sped-up note with spectrogram; step 2 (2026-10-09), von Kármán stretching.
44
5-Current state: N × N masses (N per side, 2 to 32, default 8) on vertical rails at the nodes of a square grid, spacing a, moving up and down only. Bending κ on the five-point Laplacian at each node, energy κ (L u)² / 2, force −κ L(L u). The edge nodes are supports at u = 0 with no bending term (mirror ghost), which makes the edge simply supported. Linear: no stretching yet. Damping as in the beam: drag (τ, all modes alike) and a frequency-dependent loss (dashpots alongside the grid springs, force −c L v, set by τ top at the top of the band). Parameters are given for the 8 × 8 plate (m 1 kg, κ 20 N/m, a 0.5 m, so the side is 4.5 m; τ 540 s, τ top 200 s); Masses per side refines the same plate (`refine`: m with the cell area, κ inversely, the loss coefficient c unchanged). Strike: a kick (pyramid velocity bump of half-width 0.5 m) or the shared Hertz mallet (`../common/mallet.ts`, M 1 kg, K 1000) over the same bump, at (x, y) m from the lower left corner (default (1.3, 1.1)); clicking the plate moves the strike point. The card shows the plate from above as a height map (red up, blue down, scaled to the recent largest |u|), with the strike point (cross) and the listening point (circle). The sound panel is the beam's, reduced: velocity, displacement or acceleration at a listening point (default (3.1, 1.7), bilinear between masses), pitch of the lowest mode, the whole decay as a note through the shared renderer and anti-alias filter, spectrogram with dashed lines at the distinct mode pitches, energy chart, Replay and Save WAV. Not yet: Lyapunov exponent, ghost, mode bars and histories, cross-section.
5+Current state: N × N masses (N per side, 2 to 32, default 8) on vertical rails at the nodes of a square grid, spacing a, moving up and down only. Bending κ on the five-point Laplacian at each node, energy κ (L u)² / 2, force −κ L(L u). The edge nodes are supports at u = 0 with no bending term (mirror ghost), which makes the edge simply supported. Stretching (von Kármán), below, set by the thickness h (default 5 cm; 0 = linear). Damping as in the beam: drag (τ, all modes alike) and a frequency-dependent loss (dashpots alongside the grid springs, force −c L v, set by τ top at the top of the band). Parameters are given for the 8 × 8 plate (m 1 kg, κ 20 N/m, h 0.05 m, a 0.5 m, so the side is 4.5 m; τ 540 s, τ top 200 s); Masses per side refines the same plate (`refine`: m with the cell area, κ inversely, the loss coefficient c unchanged). Strike: a kick (pyramid velocity bump of half-width 0.5 m) or the shared Hertz mallet (`../common/mallet.ts`, M 1 kg, K 1000) over the same bump, at (x, y) m from the lower left corner (default (1.3, 1.1)); clicking the plate moves the strike point. The card shows the plate from above as a height map (red up, blue down, scaled to the recent largest |u|), with the strike point (cross) and the listening point (circle). The sound panel is the beam's, reduced: velocity, displacement or acceleration at a listening point (default (3.1, 1.7), bilinear between masses), pitch of the lowest mode, the whole decay as a note through the shared renderer and anti-alias filter, spectrogram with dashed lines at the distinct mode pitches, energy chart, Replay and Save WAV. Not yet: Lyapunov exponent, ghost, mode bars and histories, cross-section.
66
77 Modes: L is diagonal in the discrete sine basis, so the modes are exact, S_mn ∝ sin(mπ(i + 1)/(N + 1)) sin(nπ(j + 1)/(N + 1)), with ω_mn = √(κ/m) (l_m + l_n), l_m = 4 sin²(mπ/(2(N + 1))) (`Modes`, sorted by frequency with (m, n)). For low modes ω_mn ∝ m² + n², the continuum plate; (m, n) and (n, m) are degenerate. Ratios to the lowest, against the continuum's 2.5 and 4.0: (1, 2) 2.309 and (2, 2) 3.618 at N = 4, 2.440 and 3.879 at N = 8, 2.483 and 3.966 at N = 16, 2.496 and 3.991 at N = 32; ω_11 stays at 1.054 to 1.089 rad/s across N. Check (2026-10-09, N = 8, kick at (1.3, 1.1), undamped, 4000 s, displacement at (3.1, 1.7)): 32 of the 33 distinct frequencies appear as spectral peaks within 0.5 cents of the Verlet-corrected prediction (2/h) asin(ωh/2) (FFT resolution 2.5 cents); the missing one, (7, 6), is barely excited there (peak 1e-4 of the largest). Energy wobble 1.2e-3 (bounded). Measured decay times match `decayTimes` (exact, L is diagonal in the modes) to the energy wobble: (1, 1) 498 s, (3, 5) 249 s, (8, 8) 146 s.
88
99 Integrator: velocity Verlet at ω h = 0.2 at the top of the band (as the beam), linear so the same step at every energy. Drag and loss in half steps around it (Strang). The loss is an explicit step p → (I − x) p, x = (c dt / m) L, whose largest eigenvalue is 2 dt / τ top (about 3e-5 by default): it never adds energy, and each mode's decay rate is off by about x / 2 relative (substeps if 2 dt / τ top passes 0.01). Unlike the beam's Crank–Nicolson step (a tridiagonal solve), a Crank–Nicolson step on the grid needs an iterative solve; with Jacobi it cost 4 to 5 times the rest of the step, for no audible difference. Cost per step: 1.2 µs at N = 8, 4.7 µs at N = 16 (Node, this laptop); the default note (4.9 s of sound, about 3400 s of plate time) renders in 1.2 s at N = 8 and 10 s at N = 16. The step shrinks as 1 / N² (ω_max / ω_11 is 32 at N = 8, 117 at N = 16, 441 at N = 32), so the cost of a note grows as N⁴; the modal integrator planned later removes the step limit.
10+
11+Stretching (von Kármán): with F the Airy stress function, strain energy (1/2Eh)|∇²F|² with ∇⁴F = −(Eh/2) l(w, w), l(f, g) = f_xx g_yy + f_yy g_xx − 2 f_xy g_xy. On the grid: s = l(u, u) from the second differences and the centred cross difference (zero beyond the edge), φ = L⁻² s by two-dimensional sine transforms (exact, with F = ∇²F = 0 at the edge), energy (ks/8) Σ s φ with ks = Eh/a² = 12(1 − ν²) κ / h² (ν = 0.3), and the force its exact gradient, −(ks/2)(Dxx(φ Dyy u) + Dyy(φ Dxx u) − 2 Dxy(φ Dxy u)). Being a gradient, Verlet conserves energy as for the linear plate. The edge conditions on F leave the edge free to move outward but not along itself (normal stresses zero, shear allowed); they were chosen because the sine transforms invert L² exactly; a fully free in-plane edge (F = ∂F/∂n = 0) would need another solver. Checks (2026-10-09): force against central differences of the energy, 1e-10 relative. Stretching over bending energy of mode (1, 1) at amplitude h: 0.257, 0.224, 0.213, 0.210, 0.2094 at N = 4, 8, 16, 32, 64, against 0.2091 from an independent continuum Galerkin sum (sine modes to 4001, Python); mode (1, 2): 0.2160 at N = 64 against 0.2150. Second order in the spacing; 7% high at the default N = 8. Energy error with a strike at 100 J (motion to 5.5 h, stretching up to 58% of the energy): 1.1e-3, the linear plate's bounded wobble. Step: STEP / √(ω_max² + 16 √(2 ks E) / m) (in-plane tension at most √(2 ks E) with all of E in one cell, twice for margin); at twice that step the energy error rises 3 to 6 times, so the margin is kept. Cost: 3.5 µs per step at N = 8, 16 µs at N = 16 (the sine transforms, 2N³ products each way using the sine matrix's mirror symmetry); the default note renders in about 5 s in Node, 3.7 s in the browser, at N = 8. The note's 4x oversampling sets the step at low energies (3.0 ms at 100 Hz, against the 5.8 ms stability needs).
12+
13+Sound (2026-10-09, N = 8, h 5 cm, mallet at (1.3, 1.1), velocity at (3.1, 1.7), lowest mode at 100 Hz): nonlinear scale 131 mJ. At 0.01 J linear. At 1 J the lowest mode starts near 140 Hz and glides down to 105 Hz within the first second; at 2.51 J (default) from about 165 Hz; at 10 J from about 200 Hz with a broadband start; at 100 J a broadband wash for the first second, the lines emerging as it decays. Spectral centroid at 0.2 s: 520, 750, 990, 1540, 2350 Hz for 0.01, 1, 2.51, 10, 100 J. The glide is downward (the stiffening relaxes as the energy decays), as in gongs that glide down.
src/plate/index.htmlmodified+5−3View file
@@ -50,6 +50,7 @@
5050 <div id="param-form">
5151 <label>m <input id="p-m" type="number" step="any" /> <span class="unit">kg</span></label>
5252 <label>κ <input id="p-kappa" type="number" step="any" /> <span class="unit">N/m (bending)</span></label>
53+ <label>h <input id="p-thick" type="number" step="any" /> <span class="unit">m (thickness, sets the stretching; 0 = linear)</span></label>
5354 <label>a <input id="p-a" type="number" step="any" /> <span class="unit">m (spacing)</span></label>
5455 <label>τ <input id="p-tau" type="number" step="any" /> <span class="unit">s (decay, 0 = off)</span></label>
5556 <label>τ top <input id="p-tautop" type="number" step="any" /> <span class="unit">s (decay of the highest modes from the loss, 0 = off)</span></label>
@@ -96,13 +97,14 @@
9697
9798 <template id="help-plate">
9899 <h2>The plate</h2>
99- <p>A square plate as a grid of equal masses on vertical rails, moving up and down only, seen from above: red is up, blue down, white at rest (scaled to fit). Each mass resists bending (κ, on how far it sits off the mean of its four neighbours). The plate is simply supported: its edge is held at zero height but free to turn. It is linear for now, so a weak and a strong strike sound the same; the nonlinearity (von Kármán stretching) comes next.</p>
100+ <p>A square plate as a grid of equal masses on vertical rails, moving up and down only, seen from above: red is up, blue down, white at rest (scaled to fit). Each mass resists bending (κ, on how far it sits off the mean of its four neighbours). The plate is simply supported: its edge is held at zero height but free to turn. A bent plate also stretches in its own plane (von Kármán), and the in-plane stress stiffens it: this is the nonlinearity. It matters once the plate moves by about its thickness h, so weak strikes ring as independent normal modes, while strong ones start sharp and glide down in pitch as they decay, and the strongest spread their energy over many modes.</p>
100101 <p>Strike gives the energy on the slider as a hammer blow (a short pyramid bump of velocity, one spacing of the 8 × 8 plate wide on each side) or with a mallet (mass M, contact stiffness K in Parameters) that presses into the plate over the same bump through a Hertz contact (force K δ^1.5) and lifts off. The point is in metres from the lower left corner; click the plate to move it there. The cross marks the strike point, the circle the listening point. The space bar pauses.</p>
101102 </template>
102103 <template id="help-params">
103104 <h2>Parameters</h2>
104- <p>Mass, bending κ, the spacing a, the decay time τ of the amplitude (drag on the masses, all modes alike; 0 for none), and τ top, the decay time a frequency-dependent loss (dashpots alongside the grid springs) gives the highest modes (at 8 × 8 masses), less for lower ones (0 for none). Apply restarts the plate at rest.</p>
105- <p>The values are for the plate with 8 masses per side. Masses per side refines the same plate: same size, strike energies and decay, with each mass scaled with the area of its cell and κ inversely, so the low modes stay put and more masses add higher ones.</p>
105+ <p>Mass, bending κ, the thickness h (it sets the stretching, as the amplitude at which it matters; 0 for a linear plate), the spacing a, the decay time τ of the amplitude (drag on the masses, all modes alike; 0 for none), and τ top, the decay time a frequency-dependent loss (dashpots alongside the grid springs) gives the highest modes (at 8 × 8 masses), less for lower ones (0 for none). Apply restarts the plate at rest.</p>
106+ <p>The values are for the plate with 8 masses per side. Masses per side refines the same plate: same size, strike energies and decay, with each mass scaled with the area of its cell and κ inversely (the thickness stays), so the low modes stay put and more masses add higher ones.</p>
107+ <p>The stretching follows the von Kármán plate equations: the in-plane stress comes from a stress function, solved on the grid at every step, held at zero with its Laplacian at the edge (the edge may move outward, not along itself). The nonlinear scale below is the energy of mode (1, 1) at which its stretching energy equals its bending energy.</p>
106108 <p>The normal modes are products of sines, numbered (m, n) by their half waves along x and y. Their frequencies go as m² + n² for the low ones (the continuum plate), so (1, 2) and (2, 1) share a pitch, and the ratios to the lowest, 2.5, 4, 5, 6.5, 8.5, ..., are not harmonic. The derived values list the lowest modes and how close the grid is to the continuum.</p>
107109 </template>
108110 <template id="help-sound">
src/plate/main.tsmodified+17−9View file
@@ -1,7 +1,7 @@
11 import "../common/style.css";
22 import "./style.css";
33 import {
4- Modes, REF_N, Verlet, copyState, decayTimes, energy, forces, hammerShape, isDamped, newState, plateSide, refine, stepFor, strikeAt,
4+ Modes, REF_N, Verlet, copyState, decayTimes, energy, forces, hammerShape, isDamped, newState, nonlinearScale, plateSide, refine, stepFor, strikeAt,
55 type Params, type State,
66 } from "./plate";
77 import { fit } from "../common/canvas";
@@ -33,10 +33,12 @@ const speedSelect = document.getElementById("speed") as HTMLSelectElement;
3333 const pauseButton = document.getElementById("pause") as HTMLButtonElement;
3434
3535 // Parameters, as the beam's defaults: bending kappa = 20 on unit masses,
36-// drag tau 540 s, the frequency-dependent loss tauTop 200 s. base: the page's
36+// drag tau 540 s, the frequency-dependent loss tauTop 200 s; thickness 5 cm
37+// for the stretching (2026-10-09: the default strike glides down by about
38+// half an octave, 100 J is a broadband wash). base: the page's
3739 // parameters, for the plate at REF_N masses per side; P: the same plate
3840 // refined to the chosen N (see refine).
39-let base: Params = { N: REF_N, m: 1, kappa: 20, a: 0.5, tau: 540, tauTop: 200 };
41+let base: Params = { N: REF_N, m: 1, kappa: 20, thick: 0.05, a: 0.5, tau: 540, tauTop: 200 };
4042 let P: Params = base;
4143 // The mallet: a physical object, not refined with N.
4244 let mallet: MalletParams = { M: 1, K: 1000 };
@@ -46,7 +48,7 @@ let mallet: MalletParams = { M: 1, K: 1000 };
4648 let modes: Modes;
4749 let s: State;
4850 let verlet: Verlet;
49-let H: number; // integrator step (s)
51+let H: number; // integrator step (s), set at each strike from its energy (see stepFor)
5052 let E0 = 0; // energy of the last strike
5153 let t = 0;
5254 let running = false;
@@ -93,7 +95,7 @@ document.getElementById("strike")!.addEventListener("click", () => {
9395 strikeLabel = `mallet${mallet.M}kg${mallet.K}-${where}`;
9496 const phi = hammerShape(P, x, y, base.a), sum = phi.reduce((a, b) => a + b, 0);
9597 const r = malletStrike(
96- { N: P.N * P.N, m: P.m, forces: (u, F) => forces(P, u, F), weights: phi.map((v) => v / sum), step: H },
98+ { N: P.N * P.N, m: P.m, forces: (u, F) => forces(P, u, F), weights: phi.map((v) => v / sum), step: stepFor(P, modes, E) },
9799 mallet, E,
98100 );
99101 s = { u: r.u, p: r.p, W: 0 };
@@ -106,6 +108,8 @@ document.getElementById("strike")!.addEventListener("click", () => {
106108 lastStrike = copyState(s);
107109 verlet = new Verlet(P, s);
108110 E0 = energy(P, s);
111+ H = stepFor(P, modes, E0);
112+ showDerived();
109113 t = 0;
110114 ampl = 0;
111115 replayButton.disabled = false;
@@ -132,6 +136,7 @@ const relErr = () => (E0 > 0 ? (energy(P, s) + s.W - E0) / E0 : 0);
132136 const paramInputs = {
133137 m: document.getElementById("p-m") as HTMLInputElement,
134138 kappa: document.getElementById("p-kappa") as HTMLInputElement,
139+ thick: document.getElementById("p-thick") as HTMLInputElement,
135140 a: document.getElementById("p-a") as HTMLInputElement,
136141 tau: document.getElementById("p-tau") as HTMLInputElement,
137142 tauTop: document.getElementById("p-tautop") as HTMLInputElement,
@@ -172,8 +177,11 @@ function showDerived() {
172177 `side ${fmtLength(plateSide(P))}, ${modes.count} modes, ${hz(w1).toPrecision(3)} … ${hz(modes.omegaMax).toPrecision(3)} Hz ` +
173178 `(top ${(modes.omegaMax / w1).toFixed(1)}× the lowest)\n` +
174179 `lowest modes (m, n):\n${rows.join("\n")}\n` +
180+ `nonlinear scale ${Number.isFinite(nonlinearScale(P)) ? fmtEnergy(nonlinearScale(P)) : "none (linear)"} (mode 1, 1: stretching = bending energy)\n` +
175181 `decay time ${decay}\n` +
176- `integrator step ${(1000 * H).toPrecision(3)} ms`;
182+ `integrator step ${(1000 * H).toPrecision(3)} ms (from the strike energy; ` +
183+ `${(1000 * stepFor(P, modes, 0)).toPrecision(3)} ms when weak, ` +
184+ `${(1000 * stepFor(P, modes, 10 ** Number(energyInput.max))).toPrecision(3)} ms at the top of the slider)`;
177185 }
178186
179187 // Sound: the last strike again, sped up so the lowest mode sounds at the
@@ -313,7 +321,7 @@ function setParams(next: Params, N: number) {
313321 base = next;
314322 P = refine(base, N);
315323 modes = new Modes(P);
316- H = stepFor(P, modes);
324+ H = stepFor(P, modes, 0);
317325 s = newState(P.N);
318326 verlet = new Verlet(P, s);
319327 E0 = 0;
@@ -344,10 +352,10 @@ nInput.addEventListener("change", () => {
344352 document.getElementById("param-apply")!.addEventListener("click", () => {
345353 const v = Object.fromEntries(Object.entries(paramInputs).map(([key, el]) => [key, Number(el.value)]));
346354 const mv = { M: Number(malletInputs.M.value), K: Number(malletInputs.K.value) };
347- const ok = Object.entries(v).every(([key, x]) => Number.isFinite(x) && (x > 0 || (["tau", "tauTop"].includes(key) && x === 0))) &&
355+ const ok = Object.entries(v).every(([key, x]) => Number.isFinite(x) && (x > 0 || (["tau", "tauTop", "thick"].includes(key) && x === 0))) &&
348356 Object.values(mv).every((x) => Number.isFinite(x) && x > 0);
349357 if (!ok) {
350- paramStatus.textContent = "parameters must be positive (τ, τ top may be 0)";
358+ paramStatus.textContent = "parameters must be positive (h, τ, τ top may be 0)";
351359 return;
352360 }
353361 paramStatus.textContent = "";
src/plate/plate.tsmodified+181−9View file
@@ -13,6 +13,23 @@
1313 // 1)) sin(n pi (j + 1) / (N + 1)), m, n = 1..N, with omega_mn = sqrt(kappa /
1414 // m) (l_m + l_n), l_m = 4 sin^2(m pi / (2 (N + 1))). For low modes l_m ~ (m
1515 // pi / (N + 1))^2, so omega_mn ~ m^2 + n^2: the continuum plate's spectrum.
16+// Stretching, the nonlinearity (von Karman): a bent plate stretches in its
17+// own plane, and the in-plane stress stiffens it. With F the Airy stress
18+// function, the strain energy is (1 / 2Eh) |Lap F|^2 with Lap^2 F = -(Eh /
19+// 2) l(w, w), l(f, g) = f_xx g_yy + f_yy g_xx - 2 f_xy g_xy. On the grid:
20+// s = l(u, u) at each node from the second differences Dxx, Dyy and the
21+// centred cross difference Dxy (zero beyond the edge), phi = L^-2 s (L^2
22+// inverted exactly by sine transforms, with F = Lap F = 0 at the edge), and
23+// the stretching energy is (ks / 8) sum s phi, ks = Eh / a^2. The force is
24+// its exact gradient, -(ks / 2) (Dxx (phi Dyy u) + Dyy (phi Dxx u) - 2 Dxy
25+// (phi Dxy u)) (the continuum's l(w, F)), so the same Verlet conserves
26+// energy as for the linear plate. The stress function's edge conditions (F
27+// = Lap F = 0) leave the edge free to move outward but not along itself
28+// (no normal stresses there, shear allowed), chosen because the sine
29+// transforms invert L^2 exactly; a fully free in-plane edge (F = dF/dn = 0)
30+// would need another solver. The stretching is set by the thickness h:
31+// ks / kappa = 12 (1 - nu^2) / h^2 (from D = Eh^3 / 12 (1 - nu^2), nu = 0.3),
32+// so it matters once the plate moves by about h.
1633 // Optional damping, as in the beam (../beam/beam.ts): drag -gamma p on each
1734 // mass (all modes decay as exp(-t / tau)), and a frequency-dependent loss,
1835 // a dashpot alongside each grid spring, force -c (L v), which damps mode mn
@@ -24,6 +41,7 @@ export interface Params {
2441 N: number; // masses per side (N^2 in all)
2542 m: number; // kg, each mass
2643 kappa: number; // N/m, bending
44+ thick: number; // m, thickness, sets the stretching (von Karman); 0 = linear
2745 a: number; // m, rest spacing (for positions and drawing; it does not enter the equations)
2846 tau: number; // s, amplitude decay time from drag (all modes); 0 = no damping
2947 tauTop: number; // s, amplitude decay time of the top of the band from the loss alone; 0 = off
@@ -33,7 +51,8 @@ export interface Params {
3351 // side; refine() gives the same plate (same size, strike energies, decay
3452 // times) with N. The spacing becomes a r, r = (REF_N + 1) / (N + 1), and
3553 // each constant follows its continuum origin: m = rho h a^2 (x r^2), kappa =
36-// D / a^2 (/ r^2), the loss dashpots c = eta (unchanged, so a wave of given
54+// D / a^2 (/ r^2), the stretching ks = Eh / a^2 (/ r^2, with the thickness
55+// unchanged), the loss dashpots c = eta (unchanged, so a wave of given
3756 // length decays at the same rate at any N). Drag does not change.
3857 export const REF_N = 8;
3958 export function refine(base: Params, N: number): Params {
@@ -76,11 +95,145 @@ export interface State {
7695 export const newState = (N: number): State => ({ u: new Float64Array(N * N), p: new Float64Array(N * N), W: 0 });
7796 export const copyState = (s: State): State => ({ u: Float64Array.from(s.u), p: Float64Array.from(s.p), W: s.W });
7897
79-// Forces F = -kappa L (L u); c (if given) receives L u.
80-export function forces(P: Params, u: Float64Array, F: Float64Array, c: Float64Array = new Float64Array(u.length)): void {
98+// Poisson's ratio, in ks / kappa only.
99+const NU = 0.3;
100+// Stretching coefficient ks (N/m^3); 0 when linear.
101+export const stretchCoefficient = (P: Params) => (P.thick > 0 ? (12 * (1 - NU * NU) * P.kappa) / (P.thick * P.thick) : 0);
102+
103+// The sine transform for N: orthonormal and symmetric, S[m][i] = sqrt(2 / (N
104+// + 1)) sin((m + 1)(i + 1) pi / (N + 1)), and 1 / (l_m + l_n)^2, the inverse
105+// of L^2 in the sine basis.
106+interface Sine { S: Float64Array; inv2: Float64Array; t: Float64Array; t2: Float64Array; e: Float64Array; o: Float64Array }
107+const sines = new Map<number, Sine>();
108+function sineFor(N: number): Sine {
109+ let z = sines.get(N);
110+ if (!z) {
111+ const S = new Float64Array(N * N), inv2 = new Float64Array(N * N);
112+ for (let m = 0; m < N; m++) for (let i = 0; i < N; i++) S[m * N + i] = Math.sqrt(2 / (N + 1)) * Math.sin(((m + 1) * (i + 1) * Math.PI) / (N + 1));
113+ for (let n = 0; n < N; n++) for (let m = 0; m < N; m++) inv2[m + N * n] = 1 / (lam1(N, m + 1) + lam1(N, n + 1)) ** 2;
114+ z = { S, inv2, t: new Float64Array(N * N), t2: new Float64Array(N * N), e: new Float64Array(N), o: new Float64Array(N) };
115+ sines.set(N, z);
116+ }
117+ return z;
118+}
119+// out = S X S for X stored row by row (x fastest); out may not be X. Each
120+// one-dimensional transform uses the mirror symmetry S[m][N - 1 - i] = (-1)^m
121+// S[m][i]: even m (0-based) see x_i + x_(N-1-i), odd m x_i - x_(N-1-i), half
122+// the products.
123+function sine2(z: Sine, N: number, X: Float64Array, out: Float64Array): void {
124+ const { S, t, e, o } = z, half = N >> 1, mid = N & 1;
125+ // t = transform of each row j (along x), stored transposed: t[j + N m].
126+ for (let j = 0; j < N; j++) {
127+ const r = N * j;
128+ for (let i = 0; i < half; i++) {
129+ e[i] = X[r + i] + X[r + N - 1 - i];
130+ o[i] = X[r + i] - X[r + N - 1 - i];
131+ }
132+ const xm = mid ? X[r + half] : 0;
133+ for (let m = 0; m < N; m++) {
134+ const Sm = m * N;
135+ let acc = 0;
136+ if (m & 1) for (let i = 0; i < half; i++) acc += o[i] * S[Sm + i];
137+ else {
138+ for (let i = 0; i < half; i++) acc += e[i] * S[Sm + i];
139+ if (mid) acc += xm * S[Sm + half];
140+ }
141+ t[j + N * m] = acc;
142+ }
143+ }
144+ // out[m + N n] = transform of column m of the result (along y).
145+ for (let m = 0; m < N; m++) {
146+ const r = N * m;
147+ for (let j = 0; j < half; j++) {
148+ e[j] = t[r + j] + t[r + N - 1 - j];
149+ o[j] = t[r + j] - t[r + N - 1 - j];
150+ }
151+ const xm = mid ? t[r + half] : 0;
152+ for (let n = 0; n < N; n++) {
153+ const Sn = n * N;
154+ let acc = 0;
155+ if (n & 1) for (let j = 0; j < half; j++) acc += o[j] * S[Sn + j];
156+ else {
157+ for (let j = 0; j < half; j++) acc += e[j] * S[Sn + j];
158+ if (mid) acc += xm * S[Sn + half];
159+ }
160+ out[m + N * n] = acc;
161+ }
162+ }
163+}
164+// phi = L^-2 s (in place in phi; s unchanged).
165+function inverseL2(N: number, s: Float64Array, phi: Float64Array): void {
166+ const z = sineFor(N);
167+ sine2(z, N, s, z.t2);
168+ for (let k = 0; k < N * N; k++) z.t2[k] *= z.inv2[k];
169+ sine2(z, N, z.t2, phi);
170+}
171+
172+// Scratch for the stretching terms, by N. Padded arrays have a ring of
173+// zeros beyond the edge, (N + 2)^2 values, node (i, j) at (i + 1) + (N + 2)(j + 1).
174+interface Stretch { pad: Float64Array; xx: Float64Array; yy: Float64Array; xy: Float64Array; s: Float64Array; phi: Float64Array; a: Float64Array; b: Float64Array; c: Float64Array }
175+const scratch = new Map<number, Stretch>();
176+function stretchScratch(N: number): Stretch {
177+ let z = scratch.get(N);
178+ if (!z) {
179+ const f = () => new Float64Array(N * N), g = () => new Float64Array((N + 2) * (N + 2));
180+ z = { pad: g(), xx: f(), yy: f(), xy: f(), s: f(), phi: f(), a: g(), b: g(), c: g() };
181+ scratch.set(N, z);
182+ }
183+ return z;
184+}
185+// Stretching energy (J) at u; with F given, its force is added to F.
186+function stretching(P: Params, u: Float64Array, F?: Float64Array): number {
187+ const ks = stretchCoefficient(P);
188+ if (ks === 0) return 0;
189+ const N = P.N, n = N * N, R = N + 2, z = stretchScratch(N);
190+ const { pad, xx, yy, xy, s, phi, a, b, c } = z;
191+ for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) pad[i + 1 + R * (j + 1)] = u[i + N * j];
192+ // Second differences at every node (zero beyond the edge): xx, yy and the
193+ // centred cross difference xy; s = l(u, u).
194+ for (let j = 0, k = 0; j < N; j++) {
195+ for (let i = 0, K = 1 + R * (j + 1); i < N; i++, k++, K++) {
196+ const cc = 2 * pad[K];
197+ xx[k] = pad[K - 1] - cc + pad[K + 1];
198+ yy[k] = pad[K - R] - cc + pad[K + R];
199+ xy[k] = 0.25 * (pad[K + R + 1] - pad[K - R + 1] - pad[K + R - 1] + pad[K - R - 1]);
200+ s[k] = 2 * (xx[k] * yy[k] - xy[k] * xy[k]);
201+ }
202+ }
203+ inverseL2(N, s, phi);
204+ let V = 0;
205+ for (let k = 0; k < n; k++) V += s[k] * phi[k];
206+ if (F) {
207+ // F -= (ks / 2) (Dxx (phi yy) + Dyy (phi xx) - 2 Dxy (phi xy)); the
208+ // difference operators are symmetric, so each is its own transpose.
209+ for (let j = 0, k = 0; j < N; j++) {
210+ for (let i = 0, K = 1 + R * (j + 1); i < N; i++, k++, K++) {
211+ a[K] = phi[k] * yy[k];
212+ b[K] = phi[k] * xx[k];
213+ c[K] = phi[k] * xy[k];
214+ }
215+ }
216+ const h = 0.5 * ks;
217+ for (let j = 0, k = 0; j < N; j++) {
218+ for (let i = 0, K = 1 + R * (j + 1); i < N; i++, k++, K++) {
219+ const g =
220+ a[K - 1] - 2 * a[K] + a[K + 1] +
221+ b[K - R] - 2 * b[K] + b[K + R] -
222+ 0.5 * (c[K + R + 1] - c[K - R + 1] - c[K + R - 1] + c[K - R - 1]);
223+ F[k] -= h * g;
224+ }
225+ }
226+ }
227+ return (ks / 8) * V;
228+}
229+
230+// Forces: bending, F = -kappa L (L u), plus stretching; c (if given) receives
231+// L u. linear: bending only.
232+export function forces(P: Params, u: Float64Array, F: Float64Array, c: Float64Array = new Float64Array(u.length), linear = false): void {
81233 laplacian(P.N, u, c);
82234 laplacian(P.N, c, F);
83235 for (let k = 0; k < F.length; k++) F[k] *= -P.kappa;
236+ if (!linear) stretching(P, u, F);
84237 }
85238
86239 export function kinetic(P: Params, s: State): number {
@@ -88,13 +241,16 @@ export function kinetic(P: Params, s: State): number {
88241 for (const x of s.p) T += x * x;
89242 return T / (2 * P.m);
90243 }
91-export function potential(P: Params, s: State): number {
92- const c = new Float64Array(s.u.length);
93- laplacian(P.N, s.u, c);
244+// Bending energy alone (the linear part) at u.
245+export function bendingEnergy(P: Params, u: Float64Array): number {
246+ const c = new Float64Array(u.length);
247+ laplacian(P.N, u, c);
94248 let V = 0;
95249 for (const x of c) V += x * x;
96250 return (P.kappa * V) / 2;
97251 }
252+export const stretchEnergy = (P: Params, u: Float64Array) => stretching(P, u);
253+export const potential = (P: Params, s: State) => bendingEnergy(P, s.u) + stretching(P, s.u);
98254 export const energy = (P: Params, s: State) => kinetic(P, s) + potential(P, s);
99255
100256 // Velocity Verlet, with drag and loss in half steps around it (Strang), as
@@ -196,10 +352,26 @@ export function decayTimes(P: Params, modes: Modes): Float64Array {
196352 });
197353 }
198354
199-// Integrator step (s): STEP / omega at the top of the band, as in the beam
200-// (linear, so the same at every energy).
355+// Nonlinear scale (J): the energy of mode (1, 1) at which its stretching
356+// energy equals its bending energy (Infinity when linear). Bending goes as
357+// the amplitude squared, stretching as its fourth power.
358+export function nonlinearScale(P: Params): number {
359+ const N = P.N, u = new Float64Array(N * N);
360+ for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) u[i + N * j] = Math.sin((Math.PI * (i + 1)) / (N + 1)) * Math.sin((Math.PI * (j + 1)) / (N + 1));
361+ const V2 = bendingEnergy(P, u), V4 = stretching(P, u);
362+ return V4 > 0 ? (V2 * V2) / V4 : Infinity;
363+}
364+
365+// Integrator step (s) for motion with energy E: STEP / omega, omega at the
366+// top of the band, as in the beam. Stretching adds an in-plane tension T
367+// (N/m), which raises the top mode's omega^2 by about 8 T / m; T is at most
368+// about sqrt(2 ks E) (all of E as stretching energy in one cell), taken
369+// twice for margin.
201370 const STEP = 0.2;
202-export const stepFor = (P: Params, modes: Modes) => STEP / modes.omegaMax;
371+export function stepFor(P: Params, modes: Modes, E: number): number {
372+ const T = Math.sqrt(2 * stretchCoefficient(P) * Math.max(E, 0));
373+ return STEP / Math.sqrt(modes.omegaMax ** 2 + (16 * T) / P.m);
374+}
203375
204376 // Positions (m from the corner at the supports (-1, -1)): node (i, j) is at
205377 // ((i + 1) a, (j + 1) a), the plate's side is (N + 1) a; nodeIndex gives the
src/plate/render-worker.tsmodified+3−2View file
@@ -1,5 +1,6 @@
11 // Renders a note: the plate from a given state, sped up, sampled as audio
2-// (the loop and energy: ../common/note-render.ts). No Lyapunov shadow yet.
2+// (the loop and energy: ../common/note-render.ts), the step re-chosen as the
3+// energy falls. No Lyapunov shadow yet.
34
45 import { Modes, Verlet, copyState, energy, isDamped, stepFor, type Params, type State } from "./plate";
56 import { renderNote, type Note } from "../common/note-render";
@@ -31,7 +32,7 @@ self.onmessage = (e: MessageEvent<RenderRequest>) => {
3132 const note = renderNote(
3233 {
3334 energy: () => energy(P, s),
34- stepFor: () => stepFor(P, modes),
35+ stepFor: (E) => stepFor(P, modes, E),
3536 step: (h) => verlet.step(s, h),
3637 signal: () => signalAt(signal, P, s, verlet.force, r.x, r.y),
3738 shadow: () => {