Plate: stretching with free edges (free in-plane edge, clamped stress function by sine solves + Woodbury); free plate moved to grid vertices (edge masses on the edge, half/quarter mass), second-order bending and stretching; mallet takes per-mass masses
7 changed files+319−126
src/common/mallet.tsmodified+5−4View file
@@ -18,7 +18,7 @@ export interface MalletParams {
1818
1919 export interface Struck {
2020 N: number;
21- m: number; // kg, each mass
21+ m: number | Float64Array; // kg, each mass (one value for all, or per mass)
2222 forces(u: Float64Array, F: Float64Array): void; // the system's own (conservative) forces
2323 weights: Float64Array; // c_j, summing to 1
2424 step: number; // s, the system's own integrator step at this energy
@@ -40,7 +40,8 @@ export function rigidContactTime(mp: MalletParams, E: number): number {
4040 }
4141
4242 export function malletStrike(sys: Struck, mp: MalletParams, E: number): MalletResult {
43- const { N, m, weights: c } = sys;
43+ const { N, weights: c } = sys;
44+ const m = typeof sys.m === "number" ? new Float64Array(N).fill(sys.m) : sys.m;
4445 const u = new Float64Array(N), p = new Float64Array(N), F = new Float64Array(N);
4546 let y = 0, vy = Math.sqrt((2 * E) / mp.M); // mallet face position and velocity
4647 // Fine step: resolve the contact on both sides, the mallet (rigid-wall
@@ -49,7 +50,7 @@ export function malletStrike(sys: Struck, mp: MalletParams, E: number): MalletRe
4950 const tr = rigidContactTime(mp, E);
5051 const dMax = ((5 * mp.M * (2 * E / mp.M)) / (4 * mp.K)) ** 0.4;
5152 let inv = 0;
52- for (let j = 0; j < N; j++) inv += (c[j] * c[j]) / m;
53+ for (let j = 0; j < N; j++) inv += (c[j] * c[j]) / m[j];
5354 const h = Math.min(sys.step, tr / 400, 0.05 * Math.sqrt(1 / inv / (1.5 * mp.K * Math.sqrt(dMax))));
5455 let t = 0, peak = 0, Fc = 0;
5556 // Total force on the masses and on the mallet at the current positions.
@@ -68,7 +69,7 @@ export function malletStrike(sys: Struck, mp: MalletParams, E: number): MalletRe
6869 for (let d = 0; t < 1000 * tr; ) {
6970 for (let j = 0; j < N; j++) p[j] += 0.5 * h * F[j];
7071 vy -= (0.5 * h * Fc) / mp.M;
71- for (let j = 0; j < N; j++) u[j] += (h * p[j]) / m;
72+ for (let j = 0; j < N; j++) u[j] += (h * p[j]) / m[j];
7273 y += h * vy;
7374 d = contact();
7475 for (let j = 0; j < N; j++) p[j] += 0.5 * h * F[j];
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 and von Kármán stretching, simply supported or free (free: linear for now); modes (m, n) at m² + n², pitch glides when struck hard.</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 or free; modes (m, n) at m² + n², pitch glides when struck hard.</div></a></li>
3838 </ul>
3939 </main>
4040 </body>
src/plate/README.mdmodified+7−3View file
@@ -1,6 +1,6 @@
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; step 2 (2026-10-09), von Kármán stretching; step 3 (2026-10-09), free edges (linear for now).
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; step 3 (2026-10-09), free edges; step 4 (2026-10-09), stretching with free edges (and the free layout moved to the grid vertices).
44
55 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
@@ -12,6 +12,10 @@ Stretching (von Kármán): with F the Airy stress function, strain energy (1/2Eh
1212
1313 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.
1414
15-Free edges (step 3, 2026-10-09; Edges select, "simply supported" default or "free"): the masses sit at cell centres (x = (i + 1/2) a, side N a, as the free-free beam), nothing holds the plate. Bending then needs the full energy D/2 ((∇²w)² − 2(1 − ν)(w_xx w_yy − w_xy²)), whose Gaussian term sets the free-edge conditions (ν = 0.3). On the grid: κ/2 times Σ xx² over nodes with both x neighbours, Σ yy² (both y neighbours), 2ν Σ xx yy (both), and 2(1 − ν) Σ t² over the cells (t the cell twist); every term is non-negative, and the energy vanishes exactly for the slide and two tilts and nothing else (checked). Those are removed from u and p after a strike and each step (`removeRigid`, the removed kinetic energy to W), as for the free-free beam; the strike readout gives their share. Modes are computed: the stiffness splits into four blocks by mirror symmetry in x and y, each eigensolved (`../common/eigen.ts`, Householder + QL, new) and the three rigid modes dropped; blocks against a full eigensolve agree to 1e-14; 0.4 s at N = 32 (Node), 4 s at N = 48. Stretching is off with free edges: it needs a free in-plane edge (F = ∂F/∂n = 0), another stress solver (next step). The frequency-dependent loss uses the free graph Laplacian, its per-mode decay times from q = S L S (approximate, as for the free beam).
15+Free edges (step 3, 2026-10-09; Edges select, "simply supported" default or "free"; layout revised in step 4): nothing holds the plate. The masses sit on the vertices of the grid, x = i a, the outer ones on the edge, each carrying its share of the area (m inside, m/2 on an edge, m/4 at a corner; `masses`), side (N − 1) a. Bending needs the full energy D/2 ((∇²w)² − 2(1 − ν)(w_xx w_yy − w_xy²)), whose Gaussian term sets the free-edge conditions (ν = 0.3). On the grid: κ/2 times the area-weighted sum over nodes of xx² + yy² + 2ν xx yy, plus 2(1 − ν) Σ t² over the cells (t the cell twist). On an edge node the curvature across the edge has no stencil; it takes the value that zeroes the bending moment there (xx = −ν yy), leaving (1 − ν²) yy²; corners have no curvature term. The energy vanishes exactly for the slide and two tilts and nothing else (checked). Those are removed from u and p (mass-weighted) after a strike and each step (`removeRigid`, the removed kinetic energy to W); the strike readout gives their share. Modes: eigenvalues of M^−1/2 K M^−1/2, split into four blocks by mirror symmetry in x and y (`../common/eigen.ts`, Householder + QL), the three rigid modes dropped; blocks against a full eigensolve agree to 1e-14. The loss uses the free graph Laplacian with half-strength springs along the edges, so m M⁻¹L has its top eigenvalue exactly 8; per-mode decay times from q = Sᵀ L S (approximate, as for the free beam).
1616
17-Check against plate theory (Leissa, free square plate, ν = 0.3, ω L² √(ρh/D)): 13.468, 19.596, 24.270, 34.801 (×2), 61.093 (×2), 63.686. Grid: 12.39, 20.12, 24.46, 32.70 (×2) at N = 8; 13.21, 19.68, 24.32, 34.28 (×2), 61.16 (×2), 62.26 at N = 32; 13.30, 19.65, 24.31, 34.46 (×2), 61.16 (×2), 62.77 at N = 48. The error falls as 1/N (first order, presumably from the edge cells, which have no curvature term), and extrapolating N = 32 and 48 gives 13.47, 19.58, 24.27, 34.81, 61.15, 63.77: within 0.1% of the table except the 61 pair (0.1%) and 63.7 (0.1%). At the default N = 8 the lowest mode is 8% low and the ratios to it are spread wider (1.63, 1.98, 2.64 against 1.455, 1.802, 2.584). A second-order edge treatment would help at small N; not done. Time domain (N = 8, kick, undamped, 1000 s): 45 of 46 distinct frequencies within 3 cents of the Verlet-corrected prediction (median 0.17 cents); energy balance 1.8e-3 (1.5e-3 with drag and loss), rigid part of u stays at round-off. Cost: the default note renders in 1.4 s in the browser at N = 8; 21 s in Node at N = 16 (twice the simply supported plate, from the longer bending stencil and the projections).
17+Check against plate theory (Leissa, free square plate, ν = 0.3, ω L² √(ρh/D)): 13.468, 19.596, 24.270, 34.801 (×2), 61.093 (×2), 63.686. Grid: 13.11, 18.38, 22.84, 32.64 (×2), 53.52 (×2), 58.53 at N = 9; 13.38, 19.27, 23.90, 34.23 (×2), 58.99 (×2), 62.29 at N = 17; 13.45, 19.51, 24.18, 34.66 (×2), 60.55 (×2), 63.33 at N = 33. Second order: the lowest mode is 2.6%, 0.67%, 0.17% low. Time domain (N = 8, kick, linear, undamped, 1000 s): all 46 distinct frequencies within 3 cents of the Verlet-corrected prediction (median 0.14 cents). The first version of free edges (c86ceb3) had the masses at cell centres, as the free-free beam; no stencil then reaches the outer half cells, which made it first order (lowest mode 8% low at N = 8) and the stretching much worse (below).
18+
19+Stretching with free edges (step 4, 2026-10-09): the edge is free in its own plane too (no normal or shear stress), F = ∂F/∂n = 0, so ∇⁴F has clamped edge conditions. s = l(u, u) at the interior nodes; φ = B⁻¹ s with B the clamped biharmonic on them, F = 0 on the edge nodes and an even ghost beyond (the standard 13-point stencil: 20 inside, 21 next to an edge, 22 at a corner): B = L² + 2C, C counting the edges a node touches. B differs from L² only on the 4(N − 3) nodes next to the edge, so B⁻¹ comes from two sine solves and a small dense capacitance matrix (Woodbury, computed once per N); against a dense solve 1e-13. Energy (ks/8) Σ s φ, force its exact gradient (checked against central differences, 1e-10), invariant under rigid motion (no net force or torque, 1e-16); a cylindrical bend (developable) has exactly zero stretching energy. Independent check: a spherical cap u = c r² (c = 0.01 /m) on a free 4 m square, against a continuum Galerkin solution of the clamped biharmonic (Legendre-based polynomials, degree 27, Python): 0.0027848 J; grid +20.3%, +5.5%, +1.4%, +0.36% at N = 9, 17, 33, 65 (second order). With the cell-centred layout the same test gave −44%, −28%, −16%, −8.7% at N = 8, 16, 32, 64 (in-plane edge on the outer masses), or +70% … +10% with it one ring outside; that is why the layout was changed. Energy balance at 100 J (stretching up to 40% of the energy): 5e-4. Nonlinear scale (lowest mode, h 5 cm, N = 8): 0.140 J. Cost per step at N = 8: 1.7 µs linear, 4.4 µs with stretching, 6.0 µs with the loss too (simply supported: 4.8 µs); the default note renders in about 6.5 s in Node, 5.5 s in the browser (simply supported: 4.3 s).
20+
21+Sound (2026-10-09, free, N = 8, h 5 cm, mallet at (1.3, 1.1), velocity at (3.1, 1.7), lowest mode at 100 Hz): at 0.01 J linear; at 1 J the energy already spreads between modes within the first second; at 10 J a broadband attack; at 100 J a wash for about the first second, the low modes emerging as it decays. Spectral centroid at 0.2 s: 880, 1210, 1270, 1290, 1500 Hz for 0.01, 1, 2.51, 10, 100 J.
src/plate/index.htmlmodified+3−3View file
@@ -102,15 +102,15 @@
102102
103103 <template id="help-plate">
104104 <h2>The plate</h2>
105- <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 edges are simply supported (held at zero height but free to turn there) or free. A free plate would also slide and tilt as a whole; that motion is removed (the plate is followed in its own frame), and the readout says how much of a strike's energy it took. Free edges have no stretching yet (it needs a free in-plane edge, next), so the free plate is linear for now. 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>
105+ <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 edges are simply supported (held at zero height but free to turn there) or free. A free plate would also slide and tilt as a whole; that motion is removed (the plate is followed in its own frame), and the readout says how much of a strike's energy it took. 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>
106106 <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>
107107 </template>
108108 <template id="help-params">
109109 <h2>Parameters</h2>
110110 <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>
111111 <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>
112- <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>
113- <p>With free edges the bending energy also includes the plate's twisting term (Poisson's ratio 0.3), which sets the edge conditions; the modes are computed, and their frequencies converge to plate theory (Leissa's tables) more slowly with N than for the simply supported plate (the lowest is 8% low at 8 × 8, 1.3% at 48 × 48). For a simply supported plate, 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>
112+ <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. A simply supported edge may move outward but not along itself; a free edge is free in its own plane too. The nonlinear scale below is the energy of the lowest mode at which its stretching energy equals its bending energy.</p>
113+ <p>With free edges the bending energy also includes the plate's twisting term (Poisson's ratio 0.3), which sets the edge conditions; the modes are computed, and their frequencies converge to plate theory (Leissa's tables) as the spacing squared (the lowest is 2.6% low at 9 × 9, 0.2% at 33 × 33). On a free plate the outer masses sit on the edge, carrying half a cell's mass (a quarter at the corners), so the free plate's side is one spacing less than the number of masses. For a simply supported plate, 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>
114114 </template>
115115 <template id="help-sound">
116116 <h2>Sped up · listen</h2>
src/plate/main.tsmodified+7−10View file
@@ -2,7 +2,7 @@ import "../common/style.css";
22 import "./style.css";
33 import {
44 Modes, REF_N, Verlet, copyState, decayTimes, energy, forces, hammerShape, isDamped, newState, nodePosition, nonlinearScale, plateSide, refine,
5- removeRigid, stepFor, strikeAt, valueAt, type Edges, type Params, type State,
5+ masses, removeRigid, stepFor, strikeAt, valueAt, type Edges, type Params, type State,
66 } from "./plate";
77 import { fit } from "../common/canvas";
88 import { newCursor } from "../common/cursor";
@@ -97,22 +97,22 @@ document.getElementById("strike")!.addEventListener("click", () => {
9797 strikeLabel = where;
9898 s = strikeAt(P, x, y, base.a, E);
9999 if (free) {
100- removeRigid(P.N, s.u);
101- const rigid = removeRigid(P.N, s.p) / (2 * P.m) / E;
100+ removeRigid(P, s.u);
101+ const rigid = removeRigid(P, s.p, true) / E;
102102 strikeInfo.textContent = `kick: ${rigidText(rigid) || "no rigid motion"}`;
103103 }
104104 } else {
105105 strikeLabel = `mallet${mallet.M}kg${mallet.K}-${where}`;
106106 const phi = hammerShape(P, x, y, base.a), sum = phi.reduce((a, b) => a + b, 0);
107107 const r = malletStrike(
108- { 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) },
108+ { N: P.N * P.N, m: masses(P), forces: (u, F) => forces(P, u, F), weights: phi.map((v) => v / sum), step: stepFor(P, modes, E) },
109109 mallet, E,
110110 );
111111 s = { u: r.u, p: r.p, W: 0 };
112112 const whole = energy(P, s) / E; // including any rigid motion
113113 if (free) {
114- removeRigid(P.N, s.u);
115- removeRigid(P.N, s.p);
114+ removeRigid(P, s.u);
115+ removeRigid(P, s.p, true);
116116 }
117117 const share = energy(P, s) / E, rigid = whole - share;
118118 strikeInfo.textContent =
@@ -203,9 +203,7 @@ function showDerived() {
203203 `side ${fmtLength(plateSide(P))}, ${modes.count} modes, ${hz(w1).toPrecision(3)} … ${hz(modes.omegaMax).toPrecision(3)} Hz ` +
204204 `(top ${(modes.omegaMax / w1).toFixed(1)}× the lowest)\n` +
205205 `lowest modes${free ? " (free edges, rigid motions removed)" : " (m, n)"}:\n${rows.join("\n")}\n` +
206- (free
207- ? "stretching off with free edges (for now)\n"
208- : `nonlinear scale ${Number.isFinite(nonlinearScale(P)) ? fmtEnergy(nonlinearScale(P)) : "none (linear)"} (mode 1, 1: stretching = bending energy)\n`) +
206+ `nonlinear scale ${Number.isFinite(nonlinearScale(P, modes)) ? fmtEnergy(nonlinearScale(P, modes)) : "none (linear)"} (${free ? "lowest mode" : "mode 1, 1"}: stretching = bending energy)\n` +
209207 `decay time ${decay}\n` +
210208 `integrator step ${(1000 * H).toPrecision(3)} ms (from the strike energy; ` +
211209 `${(1000 * stepFor(P, modes, 0)).toPrecision(3)} ms when weak, ` +
@@ -345,7 +343,6 @@ function setParams(next: Params, N: number) {
345343 base = next;
346344 P = refine(base, next.edges === "free" ? Math.max(3, N) : N); // free: at least 3 per side
347345 edgesSelect.value = P.edges;
348- paramInputs.thick.disabled = P.edges === "free";
349346 modes = new Modes(P);
350347 H = stepFor(P, modes, 0);
351348 s = newState(P.N);
src/plate/plate.tsmodified+293−102View file
@@ -5,8 +5,10 @@
55 // - "simply" (simply supported): node (i, j) at x = (i + 1) a, y = (j + 1)
66 // a; the edge nodes (i or j = -1 or N) are supports at u = 0 (held, free
77 // to turn), the plate's side (N + 1) a.
8-// - "free": nothing holds the plate; the masses sit at cell centres, x = (i
9-// + 1/2) a, as in the free-free beam, the side N a.
8+// - "free": nothing holds the plate; node (i, j) at x = i a, y = j a, the
9+// edge nodes on the plate's edge, the side (N - 1) a. Each node carries
10+// the plate's mass over its share of the surrounding cells: m inside, m /
11+// 2 on an edge, m / 4 at a corner (masses, below).
1012 // Simply supported:
1113 // - bending: kappa c^2 / 2 at each node, with c = (L u) the five-point
1214 // Laplacian, 4 u(i, j) minus the four neighbours (kappa = D / a^2, D the
@@ -39,15 +41,33 @@
3941 // - 2 (1 - nu)(w_xx w_yy - w_xy^2)), whose second (Gaussian) term only
4042 // matters at a free edge and sets its natural conditions (zero moment and
4143 // zero effective shear) (nu = 0.3, Poisson's ratio). On the grid, kappa / 2
42-// times the sum of xx^2 (at nodes with both x neighbours), yy^2 (both y
43-// neighbours), 2 nu xx yy (both), and 2 (1 - nu) t^2 over the cells, t the
44-// cell twist u(i+1, j+1) - u(i+1, j) - u(i, j+1) + u(i, j); every term is a
45-// square or (with nu < 1) a positive form, and the energy is zero exactly
46-// for the rigid motions (slide and two tilts) and nothing else. Those three
47-// motions are removed from u and p (removeRigid), as for the free-free beam:
48-// the plate is followed in its own frame. The modes are computed (Modes, by
49-// symmetry blocks); they have no (m, n). Stretching is off with free edges
50-// for now: it needs a free in-plane edge, another stress solver.
44+// times the sum over the nodes, weighted by their share of the area (1, 1/2
45+// on an edge, 1/4 at a corner), of xx^2 + yy^2 + 2 nu xx yy, plus 2 (1 - nu)
46+// t^2 over the cells, t the cell twist u(i+1, j+1) - u(i+1, j) - u(i, j+1) +
47+// u(i, j). On an edge node the curvature across the edge has no stencil; it
48+// takes the value that makes the edge's bending moment zero (xx = -nu yy on
49+// an x edge), which leaves (1 - nu^2) yy^2; a corner has no curvature term.
50+// Every term is non-negative, and the energy is zero exactly for the rigid
51+// motions (slide and two tilts) and nothing else. Those three motions are
52+// removed from u and p (removeRigid, in the mass-weighted sense), as for the
53+// free-free beam: the plate is followed in its own frame. The modes are
54+// computed (Modes, by symmetry blocks); they have no (m, n). Frequencies
55+// converge to plate theory at second order (2026-10-09: lowest mode 2.6% low
56+// at N = 9, 0.17% at N = 33). (Until 2026-10-09 the free plate had its
57+// masses at cell centres, which left the outer half cells without any
58+// stencil: first order, 8% low at N = 8, and much worse for stretching.)
59+// Stretching of a free plate: its edge is free in its own plane too (no
60+// normal or shear stress), F = dF/dn = 0, so Lap^2 F has clamped edge
61+// conditions. On the grid: s = l(u, u) at the interior nodes (both
62+// neighbours in x and y; xy the centred cross difference), and phi = B^-1 s
63+// on those nodes with B the clamped biharmonic, F = 0 on the edge nodes
64+// and an even ghost beyond (the standard 13-point stencil, 21 at the
65+// centre next to an edge): B = L^2 + 2 C, L the Dirichlet Laplacian on the
66+// interior nodes and C counting the edges a node touches (from (Lap F)^2 at
67+// the ring nodes, weight 1/2). B differs from L^2 only on the 4 (N - 3)
68+// nodes next to the edge, so B^-1 comes from sine transforms and a small
69+// dense matrix (Woodbury). Energy (ks / 8) sum s phi and force its exact
70+// gradient, as for the simply supported plate.
5171 // Optional damping, as in the beam (../beam/beam.ts): drag -gamma p on each
5272 // mass (all modes decay as exp(-t / tau)), and a frequency-dependent loss,
5373 // a dashpot alongside each grid spring, force -c (L v), which damps mode mn
@@ -80,7 +100,7 @@ export interface Params {
80100 // length decays at the same rate at any N). Drag does not change.
81101 export const REF_N = 8;
82102 export function refine(base: Params, N: number): Params {
83- const span = (n: number) => (base.edges === "free" ? n : n + 1);
103+ const span = (n: number) => (base.edges === "free" ? n - 1 : n + 1);
84104 const r = span(base.N) / span(N);
85105 const P = { ...base, N, m: base.m * r * r, kappa: base.kappa / (r * r), a: base.a * r };
86106 // tauTop is set at the top of the band, which moves with N: keep c.
@@ -94,69 +114,130 @@ export const isDamped = (P: Params) => P.tau > 0 || P.tauTop > 0;
94114 // Eigenvalues of the one-dimensional second difference, l_m, m = 1..N.
95115 const lam1 = (N: number, m: number) => 4 * Math.sin((m * Math.PI) / (2 * (N + 1))) ** 2;
96116 const isFree = (P: Params) => P.edges === "free";
97-// Largest eigenvalue of L: 2 l_N; free, twice the free-free chain's, 2 x 4
98-// sin^2((N - 1) pi / 2N).
99-const laplacianMax = (P: Params) => (isFree(P) ? 8 * Math.sin(((P.N - 1) * Math.PI) / (2 * P.N)) ** 2 : 2 * lam1(P.N, P.N));
117+// Largest eigenvalue of L: 2 l_N; free, 8 (of m M^-1 L, exactly: see laplacian).
118+const laplacianMax = (P: Params) => (isFree(P) ? 8 : 2 * lam1(P.N, P.N));
119+
120+// Share of the area at node index i along one side: 1, or 1/2 on a free edge.
121+const share = (P: Params, i: number) => (isFree(P) && (i === 0 || i === P.N - 1) ? 0.5 : 1);
122+// Mass of each node (kg): m, or for a free plate m times its share of the area.
123+const massCache = new WeakMap<Params, Float64Array>();
124+export function masses(P: Params): Float64Array {
125+ let M = massCache.get(P);
126+ if (!M) {
127+ M = new Float64Array(P.N * P.N);
128+ for (let j = 0; j < P.N; j++) for (let i = 0; i < P.N; i++) M[i + P.N * j] = P.m * share(P, i) * share(P, j);
129+ massCache.set(P, M);
130+ }
131+ return M;
132+}
100133 // Loss coefficient c (kg/s); 0 when off.
101134 const lossCoefficient = (P: Params) => (P.tauTop > 0 ? (2 * P.m) / (P.tauTop * laplacianMax(P)) : 0);
102135
103-// (L x) into out: with x = 0 at the supports, or free (the graph Laplacian
104-// of the grid springs, none beyond the edge).
136+// (L x) into out: with x = 0 at the supports, or free: the graph Laplacian
137+// of the grid springs, those along an edge at half strength (each carries
138+// half a cell). With the free masses, M^-1 L is then the sum of two
139+// one-dimensional operators with half masses at the ends, eigenvalues 4
140+// sin^2(k pi / 2 (N - 1)) / m each: largest 8 / m.
105141 function laplacian(N: number, x: Float64Array, out: Float64Array, free = false): void {
106142 for (let j = 0, k = 0; j < N; j++) {
107143 for (let i = 0; i < N; i++, k++) {
108- let s = 0, deg = 0;
109- if (i > 0) (s -= x[k - 1]), deg++;
110- if (i < N - 1) (s -= x[k + 1]), deg++;
111- if (j > 0) (s -= x[k - N]), deg++;
112- if (j < N - 1) (s -= x[k + N]), deg++;
113- out[k] = s + (free ? deg : 4) * x[k];
144+ if (!free) {
145+ let s = 4 * x[k];
146+ if (i > 0) s -= x[k - 1];
147+ if (i < N - 1) s -= x[k + 1];
148+ if (j > 0) s -= x[k - N];
149+ if (j < N - 1) s -= x[k + N];
150+ out[k] = s;
151+ continue;
152+ }
153+ const wx = j === 0 || j === N - 1 ? 0.5 : 1, wy = i === 0 || i === N - 1 ? 0.5 : 1;
154+ let s = 0;
155+ if (i > 0) s += wx * (x[k] - x[k - 1]);
156+ if (i < N - 1) s += wx * (x[k] - x[k + 1]);
157+ if (j > 0) s += wy * (x[k] - x[k - N]);
158+ if (j < N - 1) s += wy * (x[k] - x[k + N]);
159+ out[k] = s;
114160 }
115161 }
116162 }
117163
118-// Rigid motions of a free plate: orthonormal slide e0 = 1 / N and tilts ex ~
119-// i - (N - 1) / 2, ey ~ j - (N - 1) / 2. removeRigid takes their parts out
120-// of v (in place) and returns the squared norm removed.
121-export function removeRigid(N: number, v: Float64Array): number {
122- const c = (N - 1) / 2, t = Math.sqrt((N * N * (N * N - 1)) / 12);
123- let a0 = 0, ax = 0, ay = 0;
124- for (let j = 0, k = 0; j < N; j++) {
125- for (let i = 0; i < N; i++, k++) {
126- a0 += v[k];
127- ax += v[k] * (i - c);
128- ay += v[k] * (j - c);
164+// Rigid motions of a free plate: slide e0 = 1 and tilts ex = i - (N - 1) / 2,
165+// ey = j - (N - 1) / 2, orthogonal under the masses (which are symmetric).
166+// removeRigid takes their parts out of the displacements u (mass-weighted
167+// projection), or, with momentum, out of the momenta p (the rigid velocity),
168+// in place; it returns the kinetic energy removed (J, momentum only).
169+// The rigid basis, cached per Params: e_b (b = slide, x tilt, y tilt) and
170+// M e_b, each scaled by 1 / sqrt(e_b^T M e_b).
171+const rigidCache = new WeakMap<Params, { e: Float64Array[]; Me: Float64Array[] }>();
172+function rigidBasis(P: Params) {
173+ let r = rigidCache.get(P);
174+ if (!r) {
175+ const N = P.N, c = (N - 1) / 2, M = masses(P), e: Float64Array[] = [], Me: Float64Array[] = [];
176+ for (let b = 0; b < 3; b++) {
177+ const v = new Float64Array(N * N);
178+ for (let j = 0, k = 0; j < N; j++) for (let i = 0; i < N; i++, k++) v[k] = b === 0 ? 1 : b === 1 ? i - c : j - c;
179+ let d = 0;
180+ for (let k = 0; k < v.length; k++) d += M[k] * v[k] * v[k];
181+ const f = 1 / Math.sqrt(d);
182+ e.push(v.map((x) => x * f));
183+ Me.push(v.map((x, k) => x * f * M[k]));
129184 }
185+ r = { e, Me };
186+ rigidCache.set(P, r);
187+ }
188+ return r;
189+}
190+export function removeRigid(P: Params, v: Float64Array, momentum = false): number {
191+ const { e, Me } = rigidBasis(P), n = v.length;
192+ // Coefficients a_b = (M e_b)^T u (displacements) or e_b^T p (momenta);
193+ // subtract a_b e_b, or a_b M e_b.
194+ const dot = momentum ? e : Me, sub = momentum ? Me : e;
195+ let a0 = 0, a1 = 0, a2 = 0;
196+ const d0 = dot[0], d1 = dot[1], d2 = dot[2];
197+ for (let k = 0; k < n; k++) {
198+ a0 += d0[k] * v[k];
199+ a1 += d1[k] * v[k];
200+ a2 += d2[k] * v[k];
130201 }
131- a0 /= N;
132- ax /= t;
133- ay /= t;
134- for (let j = 0, k = 0; j < N; j++) for (let i = 0; i < N; i++, k++) v[k] -= a0 / N + (ax * (i - c)) / t + (ay * (j - c)) / t;
135- return a0 * a0 + ax * ax + ay * ay;
202+ const s0 = sub[0], s1 = sub[1], s2 = sub[2];
203+ for (let k = 0; k < n; k++) v[k] -= a0 * s0[k] + a1 * s1[k] + a2 * s2[k];
204+ return momentum ? (a0 * a0 + a1 * a1 + a2 * a2) / 2 : 0;
136205 }
137206
138207 // Bending energy of a free plate at u; with F given, F is set to its force.
139208 function freeBending(P: Params, u: Float64Array, F?: Float64Array): number {
140- const N = P.N, K = P.kappa;
209+ const N = P.N, K = P.kappa, edge = 1 - NU * NU, Ms = masses(P), m = P.m;
141210 if (F) F.fill(0);
142211 let V = 0;
143212 for (let j = 0, k = 0; j < N; j++) {
144213 for (let i = 0; i < N; i++, k++) {
145- const inX = i > 0 && i < N - 1, inY = j > 0 && j < N - 1;
214+ const inX = i > 0 && i < N - 1, inY = j > 0 && j < N - 1, w = Ms[k] / m; // share of the area
146215 const xx = inX ? u[k - 1] - 2 * u[k] + u[k + 1] : 0;
147216 const yy = inY ? u[k - N] - 2 * u[k] + u[k + N] : 0;
148- const both = inX && inY ? 1 : 0;
149- V += xx * xx + yy * yy + 2 * NU * both * xx * yy;
217+ // Node term and its derivatives by xx and yy (halved): inside, or on
218+ // an edge with the moment-free curvature across it.
219+ let dx = 0, dy = 0;
220+ if (inX && inY) {
221+ V += w * (xx * xx + yy * yy + 2 * NU * xx * yy);
222+ dx = xx + NU * yy;
223+ dy = yy + NU * xx;
224+ } else if (inX) {
225+ V += w * edge * xx * xx;
226+ dx = edge * xx;
227+ } else if (inY) {
228+ V += w * edge * yy * yy;
229+ dy = edge * yy;
230+ }
150231 if (F) {
151- // dV/d(xx) = kappa (xx + nu yy), on u(i-1), u(i), u(i+1) as 1, -2, 1.
232+ // On u(i-1), u(i), u(i+1) as 1, -2, 1.
152233 if (inX) {
153- const A = K * (xx + NU * both * yy);
234+ const A = K * w * dx;
154235 F[k - 1] -= A;
155236 F[k] += 2 * A;
156237 F[k + 1] -= A;
157238 }
158239 if (inY) {
159- const B = K * (yy + NU * both * xx);
240+ const B = K * w * dy;
160241 F[k - N] -= B;
161242 F[k] += 2 * B;
162243 F[k + N] -= B;
@@ -188,9 +269,8 @@ export const copyState = (s: State): State => ({ u: Float64Array.from(s.u), p: F
188269
189270 // Poisson's ratio: in ks / kappa, and in the bending of a free plate.
190271 const NU = 0.3;
191-// Stretching coefficient ks (N/m^3); 0 when linear (and, for now, with free
192-// edges).
193-export const stretchCoefficient = (P: Params) => (P.thick > 0 && !isFree(P) ? (12 * (1 - NU * NU) * P.kappa) / (P.thick * P.thick) : 0);
272+// Stretching coefficient ks (N/m^3); 0 when linear.
273+export const stretchCoefficient = (P: Params) => (P.thick > 0 ? (12 * (1 - NU * NU) * P.kappa) / (P.thick * P.thick) : 0);
194274
195275 // The sine transform for N: orthonormal and symmetric, S[m][i] = sqrt(2 / (N
196276 // + 1)) sin((m + 1)(i + 1) pi / (N + 1)), and 1 / (l_m + l_n)^2, the inverse
@@ -261,6 +341,100 @@ function inverseL2(N: number, s: Float64Array, phi: Float64Array): void {
261341 sine2(z, N, z.t2, phi);
262342 }
263343
344+// Clamped biharmonic on an M x M grid (the interior nodes of a free plate,
345+// M = N - 2): phi = B^-1 s with B = L^2 + U D U^T, U picking the nodes next
346+// to the edge and D = 2 x (edges touched). Woodbury: B^-1 = A^-1 - A^-1 U
347+// (D^-1 + U^T A^-1 U)^-1 U^T A^-1, A = L^2 (sine transforms); the r x r
348+// capacitance inverse is computed once per M.
349+interface Clamped { idx: Int32Array; capInv: Float64Array; x: Float64Array; t: Float64Array; w: Float64Array; y: Float64Array; z: Float64Array }
350+const clampedSolvers = new Map<number, Clamped>();
351+function clampedFor(M: number): Clamped {
352+ let c = clampedSolvers.get(M);
353+ if (c) return c;
354+ const n = M * M, idx: number[] = [], dg: number[] = [];
355+ for (let j = 0; j < M; j++) {
356+ for (let i = 0; i < M; i++) {
357+ const touched = (i === 0 ? 1 : 0) + (i === M - 1 ? 1 : 0) + (j === 0 ? 1 : 0) + (j === M - 1 ? 1 : 0);
358+ if (touched) {
359+ idx.push(i + M * j);
360+ dg.push(2 * touched);
361+ }
362+ }
363+ }
364+ const r = idx.length, e = new Float64Array(n), col = new Float64Array(n);
365+ const cap = Array.from({ length: r }, () => new Float64Array(r));
366+ for (let b = 0; b < r; b++) {
367+ e.fill(0);
368+ e[idx[b]] = 1;
369+ inverseL2(M, e, col);
370+ for (let a = 0; a < r; a++) cap[a][b] = col[idx[a]] + (a === b ? 1 / dg[b] : 0);
371+ }
372+ // Inverse of the (symmetric positive definite) capacitance matrix.
373+ const lam = symmetricEigen(cap), capInv = new Float64Array(r * r);
374+ for (let a = 0; a < r; a++) {
375+ for (let b = 0; b < r; b++) {
376+ let acc = 0;
377+ for (let k = 0; k < r; k++) acc += (cap[a][k] * cap[b][k]) / lam[k];
378+ capInv[a * r + b] = acc;
379+ }
380+ }
381+ c = { idx: Int32Array.from(idx), capInv, x: new Float64Array(n), t: new Float64Array(n), w: new Float64Array(n), y: new Float64Array(r), z: new Float64Array(r) };
382+ clampedSolvers.set(M, c);
383+ return c;
384+}
385+function inverseClamped(M: number, s: Float64Array, phi: Float64Array): void {
386+ const c = clampedFor(M), { idx, capInv, x, t, w, y, z } = c, r = idx.length;
387+ inverseL2(M, s, x);
388+ for (let a = 0; a < r; a++) y[a] = x[idx[a]];
389+ for (let a = 0; a < r; a++) {
390+ let acc = 0;
391+ for (let b = 0; b < r; b++) acc += capInv[a * r + b] * y[b];
392+ z[a] = acc;
393+ }
394+ t.fill(0);
395+ for (let a = 0; a < r; a++) t[idx[a]] = z[a];
396+ inverseL2(M, t, w);
397+ for (let k = 0; k < M * M; k++) phi[k] = x[k] - w[k];
398+}
399+
400+// Stretching energy of a free plate at u; with F given, its force is added.
401+function freeStretching(P: Params, ks: number, u: Float64Array, F?: Float64Array): number {
402+ const N = P.N, M = N - 2, z = stretchScratch(M);
403+ const { xx, yy, xy, s, phi } = z;
404+ for (let j = 1, q = 0; j < N - 1; j++) {
405+ for (let i = 1, k = 1 + N * j; i < N - 1; i++, k++, q++) {
406+ const cc = 2 * u[k];
407+ xx[q] = u[k - 1] - cc + u[k + 1];
408+ yy[q] = u[k - N] - cc + u[k + N];
409+ xy[q] = 0.25 * (u[k + N + 1] - u[k - N + 1] - u[k + N - 1] + u[k - N - 1]);
410+ s[q] = 2 * (xx[q] * yy[q] - xy[q] * xy[q]);
411+ }
412+ }
413+ inverseClamped(M, s, phi);
414+ let V = 0;
415+ for (let q = 0; q < M * M; q++) V += s[q] * phi[q];
416+ if (F) {
417+ // F -= (ks / 2) (Dxx^T (phi yy) + Dyy^T (phi xx) - 2 Dxy^T (phi xy)), the
418+ // differences scattered back onto the masses they read.
419+ const h = 0.5 * ks;
420+ for (let j = 1, q = 0; j < N - 1; j++) {
421+ for (let i = 1, k = 1 + N * j; i < N - 1; i++, k++, q++) {
422+ const A = h * phi[q] * yy[q], B = h * phi[q] * xx[q], C = 0.5 * h * phi[q] * xy[q];
423+ F[k - 1] -= A;
424+ F[k + 1] -= A;
425+ F[k - N] -= B;
426+ F[k + N] -= B;
427+ F[k] += 2 * (A + B);
428+ F[k + N + 1] += C;
429+ F[k - N - 1] += C;
430+ F[k - N + 1] -= C;
431+ F[k + N - 1] -= C;
432+ }
433+ }
434+ }
435+ return (ks / 8) * V;
436+}
437+
264438 // Scratch for the stretching terms, by N. Padded arrays have a ring of
265439 // zeros beyond the edge, (N + 2)^2 values, node (i, j) at (i + 1) + (N + 2)(j + 1).
266440 interface Stretch { pad: Float64Array; xx: Float64Array; yy: Float64Array; xy: Float64Array; s: Float64Array; phi: Float64Array; a: Float64Array; b: Float64Array; c: Float64Array }
@@ -278,6 +452,7 @@ function stretchScratch(N: number): Stretch {
278452 function stretching(P: Params, u: Float64Array, F?: Float64Array): number {
279453 const ks = stretchCoefficient(P);
280454 if (ks === 0) return 0;
455+ if (isFree(P)) return freeStretching(P, ks, u, F);
281456 const N = P.N, n = N * N, R = N + 2, z = stretchScratch(N);
282457 const { pad, xx, yy, xy, s, phi, a, b, c } = z;
283458 for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) pad[i + 1 + R * (j + 1)] = u[i + N * j];
@@ -324,6 +499,7 @@ function stretching(P: Params, u: Float64Array, F?: Float64Array): number {
324499 export function forces(P: Params, u: Float64Array, F: Float64Array, c: Float64Array = new Float64Array(u.length), linear = false): void {
325500 if (isFree(P)) {
326501 freeBending(P, u, F);
502+ if (!linear) stretching(P, u, F);
327503 return;
328504 }
329505 laplacian(P.N, u, c);
@@ -333,9 +509,10 @@ export function forces(P: Params, u: Float64Array, F: Float64Array, c: Float64Ar
333509 }
334510
335511 export function kinetic(P: Params, s: State): number {
512+ const M = masses(P);
336513 let T = 0;
337- for (const x of s.p) T += x * x;
338- return T / (2 * P.m);
514+ for (let k = 0; k < s.p.length; k++) T += (s.p[k] * s.p[k]) / M[k];
515+ return T / 2;
339516 }
340517 // Bending energy alone (the linear part) at u.
341518 export function bendingEnergy(P: Params, u: Float64Array): number {
@@ -365,12 +542,16 @@ const LOSS_X = 0.01;
365542 export class Verlet {
366543 private P: Params;
367544 private F: Float64Array;
545+ private M: Float64Array; // node masses
368546 private c: number;
369547 private tmp: Float64Array;
548+ private v: Float64Array;
370549 constructor(P: Params, s: State) {
371550 this.P = P;
372551 this.F = new Float64Array(P.N * P.N);
373552 this.tmp = new Float64Array(P.N * P.N);
553+ this.v = new Float64Array(P.N * P.N);
554+ this.M = masses(P);
374555 this.c = lossCoefficient(P);
375556 forces(P, s.u, this.F, this.tmp);
376557 }
@@ -379,48 +560,53 @@ export class Verlet {
379560 return this.F;
380561 }
381562 step(s: State, h: number): void {
382- const { u, p } = s, F = this.F, n = u.length, m = this.P.m;
563+ const { u, p } = s, F = this.F, n = u.length, M = this.M;
383564 const gamma = dragRate(this.P);
384565 if (gamma > 0) this.drag(s, Math.exp(-0.5 * gamma * h));
385566 if (this.c > 0) this.loss(s, h / 2);
386567 for (let k = 0; k < n; k++) p[k] += 0.5 * h * F[k];
387- for (let k = 0; k < n; k++) u[k] += (h * p[k]) / m;
568+ for (let k = 0; k < n; k++) u[k] += (h * p[k]) / M[k];
388569 forces(this.P, u, F, this.tmp);
389570 for (let k = 0; k < n; k++) p[k] += 0.5 * h * F[k];
390571 if (this.c > 0) this.loss(s, h / 2);
391572 if (gamma > 0) this.drag(s, Math.exp(-0.5 * gamma * h));
392573 if (isFree(this.P)) {
393- // Keep the plate in its own frame (the loss's edge springs, round-off).
394- removeRigid(this.P.N, u);
395- s.W += removeRigid(this.P.N, p) / (2 * m);
574+ // Keep the plate in its own frame (round-off).
575+ removeRigid(this.P, u);
576+ s.W += removeRigid(this.P, p, true);
396577 }
397578 }
579+ // dp/dt = -c L v, v = p / m per node: p -> p - c dt L v.
398580 private loss(s: State, dt: number): void {
399- const { N, m, tauTop } = this.P, p = s.p, Lp = this.tmp, n = p.length;
400- const sub = Math.ceil((2 * dt) / tauTop / LOSS_X), b = (this.c * dt) / (m * sub);
401- let T0 = 0;
402- for (let k = 0; k < n; k++) T0 += p[k] * p[k];
581+ const { N, tauTop } = this.P, p = s.p, Lv = this.tmp, v = this.v, M = this.M, n = p.length;
582+ const sub = Math.ceil((2 * dt) / tauTop / LOSS_X), b = (this.c * dt) / sub, free = isFree(this.P);
583+ let dT = 0; // twice the kinetic energy removed
403584 for (let q = 0; q < sub; q++) {
404- laplacian(N, p, Lp, isFree(this.P));
405- for (let k = 0; k < n; k++) p[k] -= b * Lp[k];
585+ for (let k = 0; k < n; k++) v[k] = p[k] / M[k];
586+ laplacian(N, v, Lv, free);
587+ for (let k = 0; k < n; k++) {
588+ const p1 = p[k] - b * Lv[k];
589+ dT += (p[k] - p1) * (p[k] + p1) / M[k];
590+ p[k] = p1;
591+ }
406592 }
407- let T1 = 0;
408- for (let k = 0; k < n; k++) T1 += p[k] * p[k];
409- s.W += (T0 - T1) / (2 * m);
593+ s.W += dT / 2;
410594 }
411595 private drag(s: State, f: number): void {
412- let p2 = 0;
413- for (let k = 0; k < s.p.length; k++) {
414- p2 += s.p[k] * s.p[k];
415- s.p[k] *= f;
596+ const p = s.p, M = this.M;
597+ let T2 = 0;
598+ for (let k = 0; k < p.length; k++) {
599+ T2 += (p[k] * p[k]) / M[k];
600+ p[k] *= f;
416601 }
417- s.W += ((1 - f * f) * p2) / (2 * this.P.m);
602+ s.W += ((1 - f * f) * T2) / 2;
418603 }
419604 }
420605
421606 // The normal modes, sorted by frequency: omega (rad/s), and for a simply
422607 // supported plate the mode numbers (m, n) (ties: by m, then n); for a free
423-// plate the shapes S (unit vectors), computed, and m = n = 0.
608+// plate the displacement shapes S (normalized as S^T M S = 1, M the node
609+// masses), computed, and m = n = 0.
424610 export class Modes {
425611 N: number;
426612 count: number;
@@ -434,7 +620,7 @@ export class Modes {
434620 if (isFree(P)) {
435621 const { lambda, S } = freeModes(P);
436622 this.count = lambda.length;
437- this.omega = Float64Array.from(lambda, (l) => Math.sqrt(Math.max(l, 0) / P.m));
623+ this.omega = Float64Array.from(lambda, (l) => Math.sqrt(Math.max(l, 0)));
438624 this.S = S;
439625 this.m = new Int32Array(this.count);
440626 this.n = new Int32Array(this.count);
@@ -454,11 +640,12 @@ export class Modes {
454640 }
455641 }
456642
457-// Modes of a free plate: eigenvalues of the bending stiffness K (lambda,
458-// kg/s^2) and unit shapes, without the three rigid motions (lambda = 0). The
459-// square is symmetric under x -> -x and y -> -y, so K splits into four
460-// blocks, even or odd in x and in y (built on mirror pairs of nodes), each
461-// a quarter the size: 1/16 of the work of one full eigensolve.
643+// Modes of a free plate: eigenvalues lambda = omega^2 of M^-1/2 K M^-1/2 (K
644+// the bending stiffness, M the node masses) and displacement shapes, without
645+// the three rigid motions (lambda = 0). The square is symmetric under x ->
646+// -x and y -> -y (masses included), so the matrix splits into four blocks,
647+// even or odd in x and in y (built on mirror pairs of nodes), each a quarter
648+// the size: 1/16 of the work of one full eigensolve.
462649 function freeModes(P: Params): { lambda: Float64Array; S: Float64Array[] } {
463650 const N = P.N, n = N * N;
464651 // One-dimensional symmetric (even) or antisymmetric (odd) unit vectors on mirror pairs.
@@ -473,7 +660,7 @@ function freeModes(P: Params): { lambda: Float64Array; S: Float64Array[] } {
473660 return out;
474661 };
475662 const found: { l: number; v: Float64Array }[] = [];
476- const e = new Float64Array(n), F = new Float64Array(n);
663+ const e = new Float64Array(n), F = new Float64Array(n), rs = Float64Array.from(masses(P), (x) => 1 / Math.sqrt(x));
477664 for (const ox of [false, true]) {
478665 for (const oy of [false, true]) {
479666 const bx = basis1(ox), by = basis1(oy);
@@ -484,11 +671,11 @@ function freeModes(P: Params): { lambda: Float64Array; S: Float64Array[] } {
484671 const A = Array.from({ length: d }, () => new Float64Array(d));
485672 for (let q = 0; q < d; q++) {
486673 e.fill(0);
487- for (const [k, w] of vecs[q]) e[k] = w;
488- freeBending(P, e, F); // F = -K e
674+ for (const [k, w] of vecs[q]) e[k] = w * rs[k];
675+ freeBending(P, e, F); // F = -K M^-1/2 b_q
489676 for (let r = 0; r < d; r++) {
490677 let acc = 0;
491- for (const [k, w] of vecs[r]) acc -= w * F[k];
678+ for (const [k, w] of vecs[r]) acc -= w * rs[k] * F[k];
492679 A[r][q] = acc;
493680 }
494681 }
@@ -497,7 +684,7 @@ function freeModes(P: Params): { lambda: Float64Array; S: Float64Array[] } {
497684 const v = new Float64Array(n);
498685 for (let q = 0; q < d; q++) {
499686 const a = A[q][c];
500- for (const [k, w] of vecs[q]) v[k] += a * w;
687+ for (const [k, w] of vecs[q]) v[k] += a * w * rs[k]; // displacement shape M^-1/2 b
501688 }
502689 found.push({ l: lam[c], v });
503690 }
@@ -509,30 +696,34 @@ function freeModes(P: Params): { lambda: Float64Array; S: Float64Array[] } {
509696 }
510697
511698 // Amplitude decay time (s) of each mode, in Modes order (Infinity if
512-// undamped): 1 / (1 / tau + q / (lambdaMax tauTop)), q = l_m + l_n (exact:
513-// L is diagonal in the modes); free, q = S L S for each shape (L is not
699+// undamped): 1 / (1 / tau + c q / 2), q = S^T L S for the mode's shape
700+// normalized by S^T M S = 1: (l_m + l_n) / m when simply supported (exact:
701+// L is diagonal in the modes); free, from the computed shape (L is not
514702 // diagonal in the free modes, so this is approximate, as for the beam).
515703 export function decayTimes(P: Params, modes: Modes): Float64Array {
516- const lmax = laplacianMax(P), Ls = new Float64Array(P.N * P.N);
704+ const c = lossCoefficient(P), Ls = new Float64Array(P.N * P.N);
517705 return Float64Array.from(modes.omega, (_, k) => {
518706 let q = 0;
519707 if (modes.S) {
520708 const v = modes.S[k];
521709 laplacian(P.N, v, Ls, true);
522710 for (let i = 0; i < v.length; i++) q += v[i] * Ls[i];
523- } else q = lam1(P.N, modes.m[k]) + lam1(P.N, modes.n[k]);
524- const rate = (P.tau > 0 ? 1 / P.tau : 0) + (P.tauTop > 0 ? q / (lmax * P.tauTop) : 0);
711+ } else q = (lam1(P.N, modes.m[k]) + lam1(P.N, modes.n[k])) / P.m;
712+ const rate = (P.tau > 0 ? 1 / P.tau : 0) + (c * q) / 2;
525713 return rate > 0 ? 1 / rate : Infinity;
526714 });
527715 }
528716
529-// Nonlinear scale (J): the energy of mode (1, 1) at which its stretching
530-// energy equals its bending energy (Infinity when linear). Bending goes as
531-// the amplitude squared, stretching as its fourth power.
532-export function nonlinearScale(P: Params): number {
717+// Nonlinear scale (J): the energy of the lowest mode ((1, 1) when simply
718+// supported) at which its stretching energy equals its bending energy
719+// (Infinity when linear). Bending goes as the amplitude squared, stretching
720+// as its fourth power.
721+export function nonlinearScale(P: Params, modes: Modes): number {
533722 if (stretchCoefficient(P) === 0) return Infinity;
534- const N = P.N, u = new Float64Array(N * N);
535- 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));
723+ const N = P.N;
724+ let u = new Float64Array(N * N);
725+ if (modes.S) u = Float64Array.from(modes.S[0]);
726+ else 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));
536727 const V2 = bendingEnergy(P, u), V4 = stretching(P, u);
537728 return V4 > 0 ? (V2 * V2) / V4 : Infinity;
538729 }
@@ -549,12 +740,12 @@ export function stepFor(P: Params, modes: Modes, E: number): number {
549740 }
550741
551742 // Positions (m from the plate's lower left corner): node (i, j) is at ((i +
552-// 1) a, (j + 1) a) and the side is (N + 1) a (simply supported), or at ((i +
553-// 1/2) a, (j + 1/2) a) and N a (free); nodeIndex gives the fractional index
554-// along one side at x.
555-const offset = (P: Params) => (isFree(P) ? 0.5 : 1);
743+// 1) a, (j + 1) a) and the side is (N + 1) a (simply supported), or at (i a,
744+// j a) and (N - 1) a (free); nodeIndex gives the fractional index along one
745+// side at x.
746+const offset = (P: Params) => (isFree(P) ? 0 : 1);
556747 export const nodePosition = (P: Params, i: number) => (i + offset(P)) * P.a;
557-export const plateSide = (P: Params) => (isFree(P) ? P.N : P.N + 1) * P.a;
748+export const plateSide = (P: Params) => (isFree(P) ? P.N - 1 : P.N + 1) * P.a;
558749 export const nodeIndex = (P: Params, x: number) => x / P.a - offset(P);
559750
560751 // The hammer's shape: a pyramid bump, (1 - |x_i - x| / w)(1 - |y_j - y| / w)
@@ -580,11 +771,11 @@ export function hammerShape(P: Params, x: number, y: number, w: number): Float64
580771
581772 // A hammer blow at (x, y): velocity in the hammer's shape, energy E.
582773 export function strikeAt(P: Params, x: number, y: number, w: number, E: number): State {
583- const s = newState(P.N), phi = hammerShape(P, x, y, w);
774+ const s = newState(P.N), phi = hammerShape(P, x, y, w), M = masses(P);
584775 let S = 0;
585- for (const v of phi) S += v * v;
586- const scale = Math.sqrt((2 * P.m * E) / S);
587- for (let k = 0; k < phi.length; k++) s.p[k] = phi[k] * scale;
776+ for (let k = 0; k < phi.length; k++) S += M[k] * phi[k] * phi[k];
777+ const scale = Math.sqrt((2 * E) / S);
778+ for (let k = 0; k < phi.length; k++) s.p[k] = M[k] * phi[k] * scale;
588779 return s;
589780 }
590781
src/plate/signals.tsmodified+3−3View file
@@ -2,7 +2,7 @@
22 // the corner), bilinear between the nearest masses, so the same listening
33 // point serves every N.
44
5-import { valueAt, type Params, type State } from "./plate";
5+import { masses, valueAt, type Params, type State } from "./plate";
66
77 export interface Signal {
88 key: string;
@@ -12,9 +12,9 @@ export interface Signal {
1212 }
1313
1414 export const SIGNALS: Signal[] = [
15- { key: "v", label: "velocity at", value: (P, s, _F, k) => s.p[k] / P.m },
15+ { key: "v", label: "velocity at", value: (P, s, _F, k) => s.p[k] / masses(P)[k] },
1616 { key: "u", label: "displacement at", value: (_P, s, _F, k) => s.u[k] },
17- { key: "a", label: "acceleration at", value: (P, _s, F, k) => F[k] / P.m },
17+ { key: "a", label: "acceleration at", value: (P, _s, F, k) => F[k] / masses(P)[k] },
1818 ];
1919
2020 export const signalByKey = (key: string) => SIGNALS.find((s) => s.key === key) ?? SIGNALS[0];