Plate: free edges (full bending energy with Poisson term, rigid motions removed, modes by symmetry blocks with a new Householder/QL eigensolver); stretching off when free for now
6 changed files+453−78
src/common/eigen.tsadded+157−0View file
@@ -0,0 +1,157 @@
1+// Eigenvalues and eigenvectors of a real symmetric matrix: Householder
2+// reduction to tridiagonal form, then the implicit QL algorithm (tred2 and
3+// tql2 of EISPACK, as in JAMA). O(n^3), a few times faster than cyclic
4+// Jacobi. A (n rows) is overwritten by the eigenvectors, eigenvector k in
5+// column k; the eigenvalues are returned in ascending order, the columns
6+// sorted with them.
7+
8+export function symmetricEigen(V: Float64Array[]): Float64Array {
9+ const n = V.length;
10+ const d = new Float64Array(n), e = new Float64Array(n);
11+ tred2(V, d, e);
12+ tql2(V, d, e);
13+ // Sort ascending (selection sort on the values, swapping columns).
14+ for (let i = 0; i < n - 1; i++) {
15+ let k = i;
16+ for (let j = i + 1; j < n; j++) if (d[j] < d[k]) k = j;
17+ if (k !== i) {
18+ const t = d[i];
19+ d[i] = d[k];
20+ d[k] = t;
21+ for (let r = 0; r < n; r++) {
22+ const x = V[r][i];
23+ V[r][i] = V[r][k];
24+ V[r][k] = x;
25+ }
26+ }
27+ }
28+ return d;
29+}
30+
31+function tred2(V: Float64Array[], d: Float64Array, e: Float64Array): void {
32+ const n = V.length;
33+ for (let j = 0; j < n; j++) d[j] = V[n - 1][j];
34+ for (let i = n - 1; i > 0; i--) {
35+ let scale = 0, h = 0;
36+ for (let k = 0; k < i; k++) scale += Math.abs(d[k]);
37+ if (scale === 0) {
38+ e[i] = d[i - 1];
39+ for (let j = 0; j < i; j++) {
40+ d[j] = V[i - 1][j];
41+ V[i][j] = 0;
42+ V[j][i] = 0;
43+ }
44+ } else {
45+ for (let k = 0; k < i; k++) {
46+ d[k] /= scale;
47+ h += d[k] * d[k];
48+ }
49+ let f = d[i - 1], g = Math.sqrt(h);
50+ if (f > 0) g = -g;
51+ e[i] = scale * g;
52+ h -= f * g;
53+ d[i - 1] = f - g;
54+ for (let j = 0; j < i; j++) e[j] = 0;
55+ for (let j = 0; j < i; j++) {
56+ f = d[j];
57+ V[j][i] = f;
58+ g = e[j] + V[j][j] * f;
59+ for (let k = j + 1; k <= i - 1; k++) {
60+ g += V[k][j] * d[k];
61+ e[k] += V[k][j] * f;
62+ }
63+ e[j] = g;
64+ }
65+ f = 0;
66+ for (let j = 0; j < i; j++) {
67+ e[j] /= h;
68+ f += e[j] * d[j];
69+ }
70+ const hh = f / (h + h);
71+ for (let j = 0; j < i; j++) e[j] -= hh * d[j];
72+ for (let j = 0; j < i; j++) {
73+ f = d[j];
74+ g = e[j];
75+ for (let k = j; k <= i - 1; k++) V[k][j] -= f * e[k] + g * d[k];
76+ d[j] = V[i - 1][j];
77+ V[i][j] = 0;
78+ }
79+ }
80+ d[i] = h;
81+ }
82+ // Accumulate the transformations.
83+ for (let i = 0; i < n - 1; i++) {
84+ V[n - 1][i] = V[i][i];
85+ V[i][i] = 1;
86+ const h = d[i + 1];
87+ if (h !== 0) {
88+ for (let k = 0; k <= i; k++) d[k] = V[k][i + 1] / h;
89+ for (let j = 0; j <= i; j++) {
90+ let g = 0;
91+ for (let k = 0; k <= i; k++) g += V[k][i + 1] * V[k][j];
92+ for (let k = 0; k <= i; k++) V[k][j] -= g * d[k];
93+ }
94+ }
95+ for (let k = 0; k <= i; k++) V[k][i + 1] = 0;
96+ }
97+ for (let j = 0; j < n; j++) {
98+ d[j] = V[n - 1][j];
99+ V[n - 1][j] = 0;
100+ }
101+ V[n - 1][n - 1] = 1;
102+ e[0] = 0;
103+}
104+
105+function tql2(V: Float64Array[], d: Float64Array, e: Float64Array): void {
106+ const n = V.length;
107+ for (let i = 1; i < n; i++) e[i - 1] = e[i];
108+ e[n - 1] = 0;
109+ let f = 0, tst1 = 0;
110+ const eps = 2 ** -52;
111+ for (let l = 0; l < n; l++) {
112+ tst1 = Math.max(tst1, Math.abs(d[l]) + Math.abs(e[l]));
113+ let m = l;
114+ while (m < n - 1 && Math.abs(e[m]) > eps * tst1) m++;
115+ if (m > l) {
116+ do {
117+ let g = d[l];
118+ let p = (d[l + 1] - g) / (2 * e[l]);
119+ let r = Math.hypot(p, 1);
120+ if (p < 0) r = -r;
121+ d[l] = e[l] / (p + r);
122+ d[l + 1] = e[l] * (p + r);
123+ const dl1 = d[l + 1];
124+ let h = g - d[l];
125+ for (let i = l + 2; i < n; i++) d[i] -= h;
126+ f += h;
127+ p = d[m];
128+ let c = 1, c2 = c, c3 = c, s = 0, s2 = 0;
129+ const el1 = e[l + 1];
130+ for (let i = m - 1; i >= l; i--) {
131+ c3 = c2;
132+ c2 = c;
133+ s2 = s;
134+ g = c * e[i];
135+ h = c * p;
136+ r = Math.hypot(p, e[i]);
137+ e[i + 1] = s * r;
138+ s = e[i] / r;
139+ c = p / r;
140+ p = c * d[i] - s * g;
141+ d[i + 1] = h + s * (c * g + s * d[i]);
142+ for (let k = 0; k < n; k++) {
143+ const row = V[k];
144+ h = row[i + 1];
145+ row[i + 1] = s * row[i] + c * h;
146+ row[i] = c * row[i] - s * h;
147+ }
148+ }
149+ p = (-s * s2 * c3 * el1 * e[l]) / dl1;
150+ e[l] = s * p;
151+ d[l] = c * p;
152+ } while (Math.abs(e[l]) > eps * tst1);
153+ }
154+ d[l] += f;
155+ e[l] = 0;
156+ }
157+}
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; 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 (free: linear for now); modes (m, n) at m² + n², pitch glides when struck hard.</div></a></li>
3838 </ul>
3939 </main>
4040 </body>
src/plate/README.mdmodified+5−1View 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.
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).
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
@@ -11,3 +11,7 @@ Integrator: velocity Verlet at ω h = 0.2 at the top of the band (as the beam),
1111 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).
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.
14+
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).
16+
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).
src/plate/index.htmlmodified+7−2View file
@@ -32,6 +32,11 @@
3232 <div class="row">
3333 <span class="muted">Masses per side</span>
3434 <input id="n" type="number" min="2" max="32" step="1" value="8" />
35+ <span class="muted" style="margin-left: 12px">Edges</span>
36+ <select id="edges">
37+ <option value="simply" selected title="held at the edge, free to turn there">simply supported</option>
38+ <option value="free" title="nothing holds the plate (rigid motion removed)">free</option>
39+ </select>
3540 <span class="muted" style="margin-left: 12px">Speed</span>
3641 <select id="speed">
3742 <option value="0.25">0.25×</option>
@@ -97,7 +102,7 @@
97102
98103 <template id="help-plate">
99104 <h2>The plate</h2>
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>
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>
101106 <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>
102107 </template>
103108 <template id="help-params">
@@ -105,7 +110,7 @@
105110 <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>
106111 <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>
107112 <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>
108- <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>
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>
109114 </template>
110115 <template id="help-sound">
111116 <h2>Sped up · listen</h2>
src/plate/main.tsmodified+63−36View file
@@ -1,8 +1,8 @@
11 import "../common/style.css";
22 import "./style.css";
33 import {
4- Modes, REF_N, Verlet, copyState, decayTimes, energy, forces, hammerShape, isDamped, newState, nonlinearScale, plateSide, refine, stepFor, strikeAt,
5- type Params, type State,
4+ 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,
66 } from "./plate";
77 import { fit } from "../common/canvas";
88 import { newCursor } from "../common/cursor";
@@ -25,6 +25,7 @@ const plateCanvas = document.getElementById("plate") as HTMLCanvasElement;
2525 const energyInput = document.getElementById("energy") as HTMLInputElement;
2626 const energyValue = document.getElementById("energy-value")!;
2727 const nInput = document.getElementById("n") as HTMLInputElement;
28+const edgesSelect = document.getElementById("edges") as HTMLSelectElement;
2829 const kindSelect = document.getElementById("strike-kind") as HTMLSelectElement;
2930 const strikeXInput = document.getElementById("strike-x") as HTMLInputElement;
3031 const strikeYInput = document.getElementById("strike-y") as HTMLInputElement;
@@ -38,7 +39,7 @@ const pauseButton = document.getElementById("pause") as HTMLButtonElement;
3839 // half an octave, 100 J is a broadband wash). base: the page's
3940 // parameters, for the plate at REF_N masses per side; P: the same plate
4041 // refined to the chosen N (see refine).
41-let base: Params = { N: REF_N, m: 1, kappa: 20, thick: 0.05, a: 0.5, tau: 540, tauTop: 200 };
42+let base: Params = { N: REF_N, edges: "simply", m: 1, kappa: 20, thick: 0.05, a: 0.5, tau: 540, tauTop: 200 };
4243 let P: Params = base;
4344 // The mallet: a physical object, not refined with N.
4445 let mallet: MalletParams = { M: 1, K: 1000 };
@@ -88,9 +89,18 @@ document.getElementById("strike")!.addEventListener("click", () => {
8889 strikeInfo.textContent = "";
8990 // The hammer's half-width is one spacing of the plate at REF_N, so the
9091 // blow is the same at any N.
92+ // A free plate is followed in its own frame: the slide and tilts a strike
93+ // gives it are removed (and reported).
94+ const free = P.edges === "free";
95+ const rigidText = (rigid: number) => (rigid > 0.005 ? `${(100 * rigid).toFixed(0)}% into sliding and tilting the plate (removed)` : "");
9196 if (kindSelect.value === "at") {
9297 strikeLabel = where;
9398 s = strikeAt(P, x, y, base.a, E);
99+ if (free) {
100+ removeRigid(P.N, s.u);
101+ const rigid = removeRigid(P.N, s.p) / (2 * P.m) / E;
102+ strikeInfo.textContent = `kick: ${rigidText(rigid) || "no rigid motion"}`;
103+ }
94104 } else {
95105 strikeLabel = `mallet${mallet.M}kg${mallet.K}-${where}`;
96106 const phi = hammerShape(P, x, y, base.a), sum = phi.reduce((a, b) => a + b, 0);
@@ -99,10 +109,16 @@ document.getElementById("strike")!.addEventListener("click", () => {
99109 mallet, E,
100110 );
101111 s = { u: r.u, p: r.p, W: 0 };
102- const share = energy(P, s) / E;
112+ const whole = energy(P, s) / E; // including any rigid motion
113+ if (free) {
114+ removeRigid(P.N, s.u);
115+ removeRigid(P.N, s.p);
116+ }
117+ const share = energy(P, s) / E, rigid = whole - share;
103118 strikeInfo.textContent =
104- `mallet: contact ${(1000 * r.contactTime).toPrecision(3)} ms, ${(100 * share).toFixed(0)}% of the energy into the plate` +
105- `${share < 0.995 ? `, ${(100 * (1 - share)).toFixed(0)}% leaves with the mallet` : ""}` +
119+ `mallet: contact ${(1000 * r.contactTime).toPrecision(3)} ms, ${(100 * share).toFixed(0)}% of the energy into the plate's vibration` +
120+ `${whole < 0.995 ? `, ${(100 * (1 - whole)).toFixed(0)}% leaves with the mallet` : ""}` +
121+ `${rigidText(rigid) ? `, ${rigidText(rigid)}` : ""}` +
106122 `; peak force ${r.peakForce.toPrecision(3)} N; shown from lift-off`;
107123 }
108124 lastStrike = copyState(s);
@@ -152,32 +168,44 @@ function showParams() {
152168 malletInputs.K.value = String(mallet.K);
153169 showDerived();
154170 }
155-// Values that follow from the parameters: the lowest modes with their (m, n)
156-// (degenerate pairs once), their ratios to the lowest against the
157-// continuum's (m^2 + n^2) / 2, and the decay times.
171+// Distinct mode frequencies (rad/s), each with how many modes share it
172+// (degenerate pairs, by the square's symmetry), lowest first.
173+function distinctModes(): { w: number; count: number; k: number }[] {
174+ const out: { w: number; count: number; k: number }[] = [];
175+ for (let k = 0; k < modes.count; k++) {
176+ const w = modes.omega[k], last = out[out.length - 1];
177+ if (last && w - last.w <= 1e-9 * w) last.count++;
178+ else out.push({ w, count: 1, k });
179+ }
180+ return out;
181+}
182+// Values that follow from the parameters: the lowest modes, with their (m,
183+// n) and the continuum's ratio (m^2 + n^2) / 2 when simply supported, and
184+// the decay times.
158185 function showDerived() {
159186 const w1 = modes.omega[0], hz = (w: number) => w / (2 * Math.PI);
160- const rows: string[] = [];
161- for (let k = 0; k < modes.count && rows.length < 8; k++) {
187+ const free = P.edges === "free";
188+ const rows = distinctModes().slice(0, 8).map(({ w, count, k }) => {
162189 const m = modes.m[k], n = modes.n[k];
163- if (m > n) continue; // (n, m) has the same frequency
164- rows.push(
165- ` (${m}, ${n})${m === n ? " " : " ×2 "}${hz(modes.omega[k]).toPrecision(3).padStart(7)} Hz ` +
166- `ratio ${(modes.omega[k] / w1).toFixed(3)} continuum ${((m * m + n * n) / 2).toFixed(1)}`,
190+ return (
191+ ` ${free ? "" : `(${m}, ${n})`}${count > 1 ? ` ×${count} ` : " "}${hz(w).toPrecision(3).padStart(7)} Hz ` +
192+ `ratio ${(w / w1).toFixed(3)}` + (free ? "" : ` continuum ${((m * m + n * n) / 2).toFixed(1)}`)
167193 );
168- }
194+ });
169195 const g = (x: number) => (x === 0 ? "0" : x.toPrecision(4));
170196 let decay = "none (no damping)";
171197 if (isDamped(P)) {
172198 const d = decayTimes(P, modes);
173- decay = `${fmtDuration(d[0])} (mode 1, 1) … ${fmtDuration(d[modes.count - 1])} (mode ${P.N}, ${P.N}), amplitude`;
199+ decay = `${fmtDuration(d[0])} (lowest mode) … ${fmtDuration(d[modes.count - 1])} (highest), amplitude`;
174200 }
175201 document.getElementById("params")!.textContent =
176202 (P.N === REF_N ? "" : `at ${P.N} × ${P.N}: m ${g(P.m)} kg, κ ${g(P.kappa)} N/m, a ${g(P.a)} m\n`) +
177203 `side ${fmtLength(plateSide(P))}, ${modes.count} modes, ${hz(w1).toPrecision(3)} … ${hz(modes.omegaMax).toPrecision(3)} Hz ` +
178204 `(top ${(modes.omegaMax / w1).toFixed(1)}× the lowest)\n` +
179- `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` +
205+ `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`) +
181209 `decay time ${decay}\n` +
182210 `integrator step ${(1000 * H).toPrecision(3)} ms (from the strike energy; ` +
183211 `${(1000 * stepFor(P, modes, 0)).toPrecision(3)} ms when weak, ` +
@@ -208,11 +236,7 @@ document.getElementById("volume")!.addEventListener("input", (e) =>
208236 const pitch = () => 10 ** Number(pitchInput.value); // Hz of the lowest mode
209237 const speedup = () => (pitch() * 2 * Math.PI) / modes.omega[0];
210238 // Distinct mode pitches (degenerate pairs once), for the spectrogram.
211-const modePitches = () => {
212- const f: number[] = [];
213- for (let k = 0; k < modes.count; k++) if (modes.m[k] <= modes.n[k]) f.push((modes.omega[k] / (2 * Math.PI)) * speedup());
214- return f;
215-};
239+const modePitches = () => distinctModes().map(({ w }) => (w / (2 * Math.PI)) * speedup());
216240 // The slowest mode's energy decays at 2 / tau1, so a drop by END_RATIO takes
217241 // about (tau1 / 2) ln(1 / END_RATIO) of physical time.
218242 const predictedNoteLength = () =>
@@ -319,7 +343,9 @@ replayButton.addEventListener("click", renderNote);
319343 // Set new parameters (for REF_N masses per side) and N: back to rest.
320344 function setParams(next: Params, N: number) {
321345 base = next;
322- P = refine(base, N);
346+ P = refine(base, next.edges === "free" ? Math.max(3, N) : N); // free: at least 3 per side
347+ edgesSelect.value = P.edges;
348+ paramInputs.thick.disabled = P.edges === "free";
323349 modes = new Modes(P);
324350 H = stepFor(P, modes, 0);
325351 s = newState(P.N);
@@ -349,6 +375,7 @@ nInput.addEventListener("change", () => {
349375 const N = Math.min(32, Math.max(2, Math.round(Number(nInput.value) || P.N)));
350376 setParams(base, N);
351377 });
378+edgesSelect.addEventListener("change", () => setParams({ ...base, edges: edgesSelect.value as Edges }, P.N));
352379 document.getElementById("param-apply")!.addEventListener("click", () => {
353380 const v = Object.fromEntries(Object.entries(paramInputs).map(([key, el]) => [key, Number(el.value)]));
354381 const mv = { M: Number(malletInputs.M.value), K: Number(malletInputs.K.value) };
@@ -387,16 +414,12 @@ function drawPlate() {
387414 const g = geometry();
388415 ctx.clearRect(0, 0, g.W, g.Ht);
389416 const { N } = P, L = plateSide(P);
390- const at = (i: number, j: number) => (i < 0 || j < 0 || i >= N || j >= N ? 0 : s.u[i + N * j]);
391417 const scale = ampl > 0 ? 1 / ampl : 0;
392- const d = heatImage.data;
418+ const d = heatImage.data, uk = (k: number) => s.u[k];
393419 for (let py = 0; py < RES; py++) {
394- const fy = ((RES - 1 - py + 0.5) / RES) * (N + 1) - 1; // node index, y up
395- const j0 = Math.floor(fy), ay = fy - j0;
420+ const y = ((RES - 1 - py + 0.5) / RES) * L; // y up
396421 for (let px = 0; px < RES; px++) {
397- const fx = ((px + 0.5) / RES) * (N + 1) - 1;
398- const i0 = Math.floor(fx), ax = fx - i0;
399- const u = (1 - ax) * (1 - ay) * at(i0, j0) + ax * (1 - ay) * at(i0 + 1, j0) + (1 - ax) * ay * at(i0, j0 + 1) + ax * ay * at(i0 + 1, j0 + 1);
422+ const u = valueAt(P, uk, ((px + 0.5) / RES) * L, y);
400423 const v = Math.max(-1, Math.min(1, u * scale)), q = 4 * (px + RES * py);
401424 // white to red (up) or blue (down)
402425 const c = v >= 0 ? [255, 255 - 200 * v, 255 - 210 * v] : [255 + 210 * v, 255 + 150 * v, 255 + 40 * v];
@@ -409,13 +432,14 @@ function drawPlate() {
409432 heatCtx.putImageData(heatImage, 0, 0);
410433 ctx.imageSmoothingEnabled = true;
411434 ctx.drawImage(heat, g.x0, g.y0, g.size, g.size);
412- // Grid lines through the masses, and the supported edge.
435+ // Grid lines through the masses, and the edge: solid when supported,
436+ // dashed when free.
413437 const px = (x: number) => g.x0 + (x / L) * g.size, py = (y: number) => g.y0 + g.size - (y / L) * g.size;
414438 ctx.strokeStyle = "rgba(0, 0, 0, 0.06)";
415439 ctx.lineWidth = 1;
416440 ctx.beginPath();
417- for (let i = 1; i <= N; i++) {
418- const z = i * P.a;
441+ for (let i = 0; i < N; i++) {
442+ const z = nodePosition(P, i);
419443 ctx.moveTo(px(z), py(0));
420444 ctx.lineTo(px(z), py(L));
421445 ctx.moveTo(px(0), py(z));
@@ -423,8 +447,11 @@ function drawPlate() {
423447 }
424448 ctx.stroke();
425449 ctx.strokeStyle = MUTED;
426- ctx.lineWidth = 2;
450+ const free = P.edges === "free";
451+ ctx.lineWidth = free ? 1 : 2;
452+ ctx.setLineDash(free ? [4, 4] : []);
427453 ctx.strokeRect(g.x0, g.y0, g.size, g.size);
454+ ctx.setLineDash([]);
428455 // Strike point and listening point.
429456 const sx = px(coord(strikeXInput)), sy = py(coord(strikeYInput));
430457 ctx.strokeStyle = TEXT;
src/plate/plate.tsmodified+220−38View file
@@ -1,8 +1,13 @@
11 // A square plate as N x N equal masses on vertical rails at the nodes of a
22 // grid, spacing a, moving up and down only, by u (m), stored row by row:
3-// node (i, j) at index i + N j, at x = (i + 1) a, y = (j + 1) a. The edge
4-// nodes (i or j = -1 or N) are supports at u = 0: the plate is simply
5-// supported (held at its edge, free to turn there). SI units. Linear for now:
3+// node (i, j) at index i + N j. SI units. Two kinds of edge (see "Free
4+// edges" below for the second):
5+// - "simply" (simply supported): node (i, j) at x = (i + 1) a, y = (j + 1)
6+// a; the edge nodes (i or j = -1 or N) are supports at u = 0 (held, free
7+// 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.
10+// Simply supported:
611 // - bending: kappa c^2 / 2 at each node, with c = (L u) the five-point
712 // Laplacian, 4 u(i, j) minus the four neighbours (kappa = D / a^2, D the
813 // flexural rigidity). A simply supported edge has no bending term: with
@@ -30,6 +35,19 @@
3035 // would need another solver. The stretching is set by the thickness h:
3136 // ks / kappa = 12 (1 - nu^2) / h^2 (from D = Eh^3 / 12 (1 - nu^2), nu = 0.3),
3237 // so it matters once the plate moves by about h.
38+// Free edges: the bending energy then needs its full form, D / 2 ((Lap w)^2
39+// - 2 (1 - nu)(w_xx w_yy - w_xy^2)), whose second (Gaussian) term only
40+// matters at a free edge and sets its natural conditions (zero moment and
41+// 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.
3351 // Optional damping, as in the beam (../beam/beam.ts): drag -gamma p on each
3452 // mass (all modes decay as exp(-t / tau)), and a frequency-dependent loss,
3553 // a dashpot alongside each grid spring, force -c (L v), which damps mode mn
@@ -37,8 +55,13 @@
3755 // State: displacements u, momenta p (kg m/s), and W, the energy (J)
3856 // dissipated so far (so E + W is conserved).
3957
58+import { symmetricEigen } from "../common/eigen";
59+
60+export type Edges = "simply" | "free";
61+
4062 export interface Params {
4163 N: number; // masses per side (N^2 in all)
64+ edges: Edges;
4265 m: number; // kg, each mass
4366 kappa: number; // N/m, bending
4467 thick: number; // m, thickness, sets the stretching (von Karman); 0 = linear
@@ -49,14 +72,16 @@ export interface Params {
4972
5073 // Resolution: the page's parameters describe the plate at REF_N masses per
5174 // side; refine() gives the same plate (same size, strike energies, decay
52-// times) with N. The spacing becomes a r, r = (REF_N + 1) / (N + 1), and
75+// times) with N. The spacing becomes a r, r the ratio of the sides in
76+// spacings (N + 1 simply supported, N free), and
5377 // each constant follows its continuum origin: m = rho h a^2 (x r^2), kappa =
5478 // D / a^2 (/ r^2), the stretching ks = Eh / a^2 (/ r^2, with the thickness
5579 // unchanged), the loss dashpots c = eta (unchanged, so a wave of given
5680 // length decays at the same rate at any N). Drag does not change.
5781 export const REF_N = 8;
5882 export function refine(base: Params, N: number): Params {
59- const r = (base.N + 1) / (N + 1);
83+ const span = (n: number) => (base.edges === "free" ? n : n + 1);
84+ const r = span(base.N) / span(N);
6085 const P = { ...base, N, m: base.m * r * r, kappa: base.kappa / (r * r), a: base.a * r };
6186 // tauTop is set at the top of the band, which moves with N: keep c.
6287 if (base.tauTop > 0) P.tauTop = base.tauTop * (laplacianMax(base) / laplacianMax(P)) * r * r;
@@ -68,23 +93,89 @@ export const isDamped = (P: Params) => P.tau > 0 || P.tauTop > 0;
6893
6994 // Eigenvalues of the one-dimensional second difference, l_m, m = 1..N.
7095 const lam1 = (N: number, m: number) => 4 * Math.sin((m * Math.PI) / (2 * (N + 1))) ** 2;
71-// Largest eigenvalue of L: 2 l_N.
72-const laplacianMax = (P: Params) => 2 * lam1(P.N, P.N);
96+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));
73100 // Loss coefficient c (kg/s); 0 when off.
74101 const lossCoefficient = (P: Params) => (P.tauTop > 0 ? (2 * P.m) / (P.tauTop * laplacianMax(P)) : 0);
75102
76-// (L x) into out, with x = 0 at the supports.
77-function laplacian(N: number, x: Float64Array, out: Float64Array): void {
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).
105+function laplacian(N: number, x: Float64Array, out: Float64Array, free = false): void {
106+ for (let j = 0, k = 0; j < N; j++) {
107+ 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];
114+ }
115+ }
116+}
117+
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;
78124 for (let j = 0, k = 0; j < N; j++) {
79125 for (let i = 0; i < N; i++, k++) {
80- let s = 4 * x[k];
81- if (i > 0) s -= x[k - 1];
82- if (i < N - 1) s -= x[k + 1];
83- if (j > 0) s -= x[k - N];
84- if (j < N - 1) s -= x[k + N];
85- out[k] = s;
126+ a0 += v[k];
127+ ax += v[k] * (i - c);
128+ ay += v[k] * (j - c);
86129 }
87130 }
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;
136+}
137+
138+// Bending energy of a free plate at u; with F given, F is set to its force.
139+function freeBending(P: Params, u: Float64Array, F?: Float64Array): number {
140+ const N = P.N, K = P.kappa;
141+ if (F) F.fill(0);
142+ let V = 0;
143+ for (let j = 0, k = 0; j < N; j++) {
144+ for (let i = 0; i < N; i++, k++) {
145+ const inX = i > 0 && i < N - 1, inY = j > 0 && j < N - 1;
146+ const xx = inX ? u[k - 1] - 2 * u[k] + u[k + 1] : 0;
147+ 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;
150+ if (F) {
151+ // dV/d(xx) = kappa (xx + nu yy), on u(i-1), u(i), u(i+1) as 1, -2, 1.
152+ if (inX) {
153+ const A = K * (xx + NU * both * yy);
154+ F[k - 1] -= A;
155+ F[k] += 2 * A;
156+ F[k + 1] -= A;
157+ }
158+ if (inY) {
159+ const B = K * (yy + NU * both * xx);
160+ F[k - N] -= B;
161+ F[k] += 2 * B;
162+ F[k + N] -= B;
163+ }
164+ }
165+ if (i < N - 1 && j < N - 1) {
166+ const tw = u[k + N + 1] - u[k + 1] - u[k + N] + u[k];
167+ V += 2 * (1 - NU) * tw * tw;
168+ if (F) {
169+ const C = 2 * (1 - NU) * K * tw;
170+ F[k + N + 1] -= C;
171+ F[k + 1] += C;
172+ F[k + N] += C;
173+ F[k] -= C;
174+ }
175+ }
176+ }
177+ }
178+ return (K * V) / 2;
88179 }
89180
90181 export interface State {
@@ -95,10 +186,11 @@ export interface State {
95186 export const newState = (N: number): State => ({ u: new Float64Array(N * N), p: new Float64Array(N * N), W: 0 });
96187 export const copyState = (s: State): State => ({ u: Float64Array.from(s.u), p: Float64Array.from(s.p), W: s.W });
97188
98-// Poisson's ratio, in ks / kappa only.
189+// Poisson's ratio: in ks / kappa, and in the bending of a free plate.
99190 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);
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);
102194
103195 // The sine transform for N: orthonormal and symmetric, S[m][i] = sqrt(2 / (N
104196 // + 1)) sin((m + 1)(i + 1) pi / (N + 1)), and 1 / (l_m + l_n)^2, the inverse
@@ -227,9 +319,13 @@ function stretching(P: Params, u: Float64Array, F?: Float64Array): number {
227319 return (ks / 8) * V;
228320 }
229321
230-// Forces: bending, F = -kappa L (L u), plus stretching; c (if given) receives
231-// L u. linear: bending only.
322+// Forces: bending, F = -kappa L (L u) (free: freeBending), plus stretching;
323+// c (if given) is scratch. linear: bending only.
232324 export function forces(P: Params, u: Float64Array, F: Float64Array, c: Float64Array = new Float64Array(u.length), linear = false): void {
325+ if (isFree(P)) {
326+ freeBending(P, u, F);
327+ return;
328+ }
233329 laplacian(P.N, u, c);
234330 laplacian(P.N, c, F);
235331 for (let k = 0; k < F.length; k++) F[k] *= -P.kappa;
@@ -243,6 +339,7 @@ export function kinetic(P: Params, s: State): number {
243339 }
244340 // Bending energy alone (the linear part) at u.
245341 export function bendingEnergy(P: Params, u: Float64Array): number {
342+ if (isFree(P)) return freeBending(P, u);
246343 const c = new Float64Array(u.length);
247344 laplacian(P.N, u, c);
248345 let V = 0;
@@ -292,6 +389,11 @@ export class Verlet {
292389 for (let k = 0; k < n; k++) p[k] += 0.5 * h * F[k];
293390 if (this.c > 0) this.loss(s, h / 2);
294391 if (gamma > 0) this.drag(s, Math.exp(-0.5 * gamma * h));
392+ 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);
396+ }
295397 }
296398 private loss(s: State, dt: number): void {
297399 const { N, m, tauTop } = this.P, p = s.p, Lp = this.tmp, n = p.length;
@@ -299,7 +401,7 @@ export class Verlet {
299401 let T0 = 0;
300402 for (let k = 0; k < n; k++) T0 += p[k] * p[k];
301403 for (let q = 0; q < sub; q++) {
302- laplacian(N, p, Lp);
404+ laplacian(N, p, Lp, isFree(this.P));
303405 for (let k = 0; k < n; k++) p[k] -= b * Lp[k];
304406 }
305407 let T1 = 0;
@@ -316,20 +418,32 @@ export class Verlet {
316418 }
317419 }
318420
319-// The normal modes, sorted by frequency (ties: by m, then n): mode numbers
320-// (m, n) and omega (rad/s).
421+// The normal modes, sorted by frequency: omega (rad/s), and for a simply
422+// 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.
321424 export class Modes {
322425 N: number;
323426 count: number;
324427 m: Int32Array;
325428 n: Int32Array;
326429 omega: Float64Array;
430+ S: Float64Array[] | null = null;
327431 constructor(P: Params) {
328- const N = P.N, w0 = Math.sqrt(P.kappa / P.m);
432+ const N = P.N;
433+ this.N = N;
434+ if (isFree(P)) {
435+ const { lambda, S } = freeModes(P);
436+ this.count = lambda.length;
437+ this.omega = Float64Array.from(lambda, (l) => Math.sqrt(Math.max(l, 0) / P.m));
438+ this.S = S;
439+ this.m = new Int32Array(this.count);
440+ this.n = new Int32Array(this.count);
441+ return;
442+ }
443+ const w0 = Math.sqrt(P.kappa / P.m);
329444 const all: [number, number, number][] = [];
330445 for (let m = 1; m <= N; m++) for (let n = 1; n <= N; n++) all.push([w0 * (lam1(N, m) + lam1(N, n)), m, n]);
331446 all.sort((x, y) => x[0] - y[0] || x[1] - y[1] || x[2] - y[2]);
332- this.N = N;
333447 this.count = all.length;
334448 this.omega = Float64Array.from(all, (x) => x[0]);
335449 this.m = Int32Array.from(all, (x) => x[1]);
@@ -340,13 +454,73 @@ export class Modes {
340454 }
341455 }
342456
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.
462+function freeModes(P: Params): { lambda: Float64Array; S: Float64Array[] } {
463+ const N = P.N, n = N * N;
464+ // One-dimensional symmetric (even) or antisymmetric (odd) unit vectors on mirror pairs.
465+ const basis1 = (odd: boolean): [number, number][][] => {
466+ const out: [number, number][][] = [];
467+ for (let k = 0; k < Math.ceil(N / 2); k++) {
468+ const k2 = N - 1 - k;
469+ if (k === k2) {
470+ if (!odd) out.push([[k, 1]]);
471+ } else out.push([[k, Math.SQRT1_2], [k2, odd ? -Math.SQRT1_2 : Math.SQRT1_2]]);
472+ }
473+ return out;
474+ };
475+ const found: { l: number; v: Float64Array }[] = [];
476+ const e = new Float64Array(n), F = new Float64Array(n);
477+ for (const ox of [false, true]) {
478+ for (const oy of [false, true]) {
479+ const bx = basis1(ox), by = basis1(oy);
480+ const vecs: [number, number][][] = [];
481+ for (const vy of by) for (const vx of bx) vecs.push(vy.flatMap(([j, wy]) => vx.map(([i, wx]): [number, number] => [i + N * j, wx * wy])));
482+ const d = vecs.length;
483+ if (d === 0) continue;
484+ const A = Array.from({ length: d }, () => new Float64Array(d));
485+ for (let q = 0; q < d; q++) {
486+ e.fill(0);
487+ for (const [k, w] of vecs[q]) e[k] = w;
488+ freeBending(P, e, F); // F = -K e
489+ for (let r = 0; r < d; r++) {
490+ let acc = 0;
491+ for (const [k, w] of vecs[r]) acc -= w * F[k];
492+ A[r][q] = acc;
493+ }
494+ }
495+ const lam = symmetricEigen(A);
496+ for (let c = 0; c < d; c++) {
497+ const v = new Float64Array(n);
498+ for (let q = 0; q < d; q++) {
499+ const a = A[q][c];
500+ for (const [k, w] of vecs[q]) v[k] += a * w;
501+ }
502+ found.push({ l: lam[c], v });
503+ }
504+ }
505+ }
506+ found.sort((x, y) => x.l - y.l);
507+ const kept = found.slice(3); // the rigid motions, lambda = 0 to round-off
508+ return { lambda: Float64Array.from(kept, (x) => x.l), S: kept.map((x) => x.v) };
509+}
510+
343511 // Amplitude decay time (s) of each mode, in Modes order (Infinity if
344-// undamped): 1 / (1 / tau + (l_m + l_n) / (lambdaMax tauTop)). Exact: L is
345-// diagonal in the modes.
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
514+// diagonal in the free modes, so this is approximate, as for the beam).
346515 export function decayTimes(P: Params, modes: Modes): Float64Array {
347- const lmax = laplacianMax(P);
516+ const lmax = laplacianMax(P), Ls = new Float64Array(P.N * P.N);
348517 return Float64Array.from(modes.omega, (_, k) => {
349- const q = lam1(P.N, modes.m[k]) + lam1(P.N, modes.n[k]);
518+ let q = 0;
519+ if (modes.S) {
520+ const v = modes.S[k];
521+ laplacian(P.N, v, Ls, true);
522+ 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]);
350524 const rate = (P.tau > 0 ? 1 / P.tau : 0) + (P.tauTop > 0 ? q / (lmax * P.tauTop) : 0);
351525 return rate > 0 ? 1 / rate : Infinity;
352526 });
@@ -356,6 +530,7 @@ export function decayTimes(P: Params, modes: Modes): Float64Array {
356530 // energy equals its bending energy (Infinity when linear). Bending goes as
357531 // the amplitude squared, stretching as its fourth power.
358532 export function nonlinearScale(P: Params): number {
533+ if (stretchCoefficient(P) === 0) return Infinity;
359534 const N = P.N, u = new Float64Array(N * N);
360535 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));
361536 const V2 = bendingEnergy(P, u), V4 = stretching(P, u);
@@ -373,11 +548,14 @@ export function stepFor(P: Params, modes: Modes, E: number): number {
373548 return STEP / Math.sqrt(modes.omegaMax ** 2 + (16 * T) / P.m);
374549 }
375550
376-// Positions (m from the corner at the supports (-1, -1)): node (i, j) is at
377-// ((i + 1) a, (j + 1) a), the plate's side is (N + 1) a; nodeIndex gives the
378-// fractional index along one side at x.
379-export const plateSide = (P: Params) => (P.N + 1) * P.a;
380-export const nodeIndex = (P: Params, x: number) => x / P.a - 1;
551+// 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);
556+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;
558+export const nodeIndex = (P: Params, x: number) => x / P.a - offset(P);
381559
382560 // The hammer's shape: a pyramid bump, (1 - |x_i - x| / w)(1 - |y_j - y| / w)
383561 // on the masses it covers, (x, y) its centre and w (m) its half-width (at
@@ -386,9 +564,9 @@ export function hammerShape(P: Params, x: number, y: number, w: number): Float64
386564 const N = P.N, we = Math.max(w, P.a), phi = new Float64Array(N * N);
387565 let any = false;
388566 for (let j = 0; j < N; j++) {
389- const fy = Math.max(0, 1 - Math.abs((j + 1) * P.a - y) / we);
567+ const fy = Math.max(0, 1 - Math.abs(nodePosition(P, j) - y) / we);
390568 for (let i = 0; i < N; i++) {
391- const v = fy * Math.max(0, 1 - Math.abs((i + 1) * P.a - x) / we);
569+ const v = fy * Math.max(0, 1 - Math.abs(nodePosition(P, i) - x) / we);
392570 phi[i + N * j] = v;
393571 if (v > 0) any = true;
394572 }
@@ -411,12 +589,16 @@ export function strikeAt(P: Params, x: number, y: number, w: number, E: number):
411589 }
412590
413591 // A per-node value v(k) at (x, y), bilinear between the four nearest nodes;
414-// the supports count as 0 (they do not move).
592+// the supports count as 0 (they do not move); past the outer masses of a
593+// free plate their values hold.
415594 export function valueAt(P: Params, v: (k: number) => number, x: number, y: number): number {
416595 const N = P.N, L = plateSide(P);
417596 const fx = nodeIndex(P, Math.min(Math.max(x, 0), L)), fy = nodeIndex(P, Math.min(Math.max(y, 0), L));
418597 const i0 = Math.floor(fx), j0 = Math.floor(fy), ax = fx - i0, ay = fy - j0;
419- const at = (i: number, j: number) => (i < 0 || j < 0 || i >= N || j >= N ? 0 : v(i + N * j));
598+ const cl = (i: number) => Math.min(N - 1, Math.max(0, i));
599+ const at = isFree(P)
600+ ? (i: number, j: number) => v(cl(i) + N * cl(j))
601+ : (i: number, j: number) => (i < 0 || j < 0 || i >= N || j >= N ? 0 : v(i + N * j));
420602 return (
421603 (1 - ax) * (1 - ay) * at(i0, j0) + ax * (1 - ay) * (ax ? at(i0 + 1, j0) : 0) +
422604 (1 - ax) * ay * (ay ? at(i0, j0 + 1) : 0) + ax * ay * (ax && ay ? at(i0 + 1, j0 + 1) : 0)