Driven plate: new app, rectangular simply supported von Kármán plate driven in one mode (modal force in phase with the mode's velocity, energy-target controller on a log ramp), height map, energy chart, perturbation option
8 changed files+860−0
src/driven-plate/README.mdadded+13−0View file
@@ -0,0 +1,13 @@
1+# driven-plate
2+
3+A rectangular plate driven in one mode, without damping, its energy raised slowly, to watch the energy move from the driven mode into the others. Step 1 (2026-10-09): simply supported Nx × Ny plate with von Kármán stretching, modal drive with an energy-target controller, log ramp, height map, energy chart.
4+
5+Plate (`plate.ts`): the square plate of `../plate` (simply supported, flat) generalized to a rectangle, with spacings ax and ay. Nx follows from the masses across the short side and the aspect so the cells are nearly square (default 8 across, aspect 1.3: 11 × 8). Units scaled away (D = ρh = h = 1, short side 1): with energy given against the nonlinear scale, displacement against h and time in driven periods, only the aspect and ν = 0.3 remain. Bending on the five-point Laplacian, exact sine modes ω_mn = l_m / ax² + l_n / ay²; stretching by the stress function, L⁻² by two sine transforms (plain matrix products, no mirror trick yet), force the exact gradient.
6+
7+Drive (`drive.ts`): force in the shape of the driven mode (unit shape S, node force Q √M S on the modal coordinate q = √M S·u), so in the linear plate only that mode moves. Q = G q̇, in phase with the mode's velocity; G = (E_target − E) / (τc Ed) every 1/32 period (τc = 2 periods, Ed the driven mode's linear energy, |G| ≤ 0.2 ω). Q held over each step, its work Q Δq summed; energy error = (E − E_start − W) / E. Start: the driven mode at rest, displaced, at 10⁻³ E_nl; Perturbation adds 10⁻⁶ of that in every other mode (equal energies, random phases, fixed seed). E_nl: energy of the driven mode's shape at which stretching equals bending, Vb(1)² / Vs(1). Step: as `../plate` (0.2 / ω at the top of the band, raised for the in-plane force).
8+
9+Page (`main.ts`): runs on the main thread, advancing speed × frame time in driven periods, at most 14 ms of compute per frame (the actual speed is shown). Ramp: log10 of the target rises by 1 per "periods per decade" (default 400) up to 10². Energy chart: total, target, driven mode's linear energy, sum of the other modes' linear energies (from the sine transforms of u and p), log scale in E_nl.
10+
11+Checks (2026-10-09, Node/bun scratch script): force against central differences of the energy, 5e-10 relative (Nx 9, Ny 6); stretching over bending of the (1, 1) shape at amplitude h on the square, 0.2243, 0.2134, 0.2094 at N = 8, 16, 64 (`../plate`: 0.224, 0.213, 0.2094; continuum 0.2091); linear period of (1, 1) against the Verlet-corrected prediction, 7e-9; undriven energy error over 200 periods at 0.1, 1, 10, 100 E_nl: 4e-6, 1e-4, 6e-4, 9e-4 (bounded); controller on the default ramp (10⁻³ to 10² over 2000 periods) tracks to 1%, energy error 6e-4. Cost 4.3 µs per step (11 × 8); steps per driven period 1700 at 0.1 E_nl, 6000 at 10, 18000 at 100 (the plate stiffens: the driven mode's period falls to 0.1 of its linear value at 100 E_nl), so the page keeps up with 20 periods/s only to a few E_nl.
12+
13+First look (default settings): the other modes' energy rises steadily with the drive as the driven mode's own nonlinear response (∝ E³), then jumps by about two decades near 0.5 E_nl, after which they hold a third or more of the energy; the driven mode's pitch rises (1.4 × linear at 2 E_nl). Without the perturbation the route differs (42% left in the driven mode at 1.7 E_nl against 70% at 2 E_nl with it).
src/driven-plate/drive.tsadded+199−0View file
@@ -0,0 +1,199 @@
1+// The driven plate: the plate of plate.ts pushed by a force in the shape of
2+// one mode (the driven mode), with a controller that steers the total energy
3+// toward a target, and no damping.
4+// Modal force: with the driven mode's unit shape S (sum of squares 1) and
5+// the node mass M, the modal coordinate is q = sqrt(M) S.u and its velocity
6+// qd = S.p / sqrt(M); a generalized force Q on q is the node force Q sqrt(M)
7+// S. In the linear plate it moves the driven mode and nothing else, so any
8+// other mode that moves got its energy through the stretching.
9+// Controller: Q = G qd, a push in phase with the mode's own velocity (a
10+// negative damping), so it keeps adding energy whatever the mode's
11+// frequency, which rises as the plate stiffens with amplitude. Over a cycle
12+// the mean of qd^2 is the driven mode's energy Ed = (qd^2 + omega^2 q^2) / 2
13+// (linear), so G = (Et - E) / (TAU_C Ed) brings E to the target Et at rate
14+// 1 / TAU_C; G is capped at G_MAX omega (and Ed floored) so the push stays
15+// gentle when the driven mode holds little of the energy. G is set every
16+// control interval from E; Q is held fixed over each step, and its work Q dq
17+// is added to Wd, so E - Wd stays constant up to the integrator's error.
18+// Time: velocity Verlet at the step of plate.ts, recomputed from the energy.
19+
20+import { Modes, Plate, gridFor, modeShape, nonlinearScale, omegaOf, stepFor, type Grid, type Params } from "./plate";
21+
22+const TAU_C = 2; // controller time constant, in driven periods
23+const G_MAX = 0.2; // cap on G, times the driven mode's omega
24+const PERTURB = 1e-6; // energy of the perturbation, relative to the starting energy
25+
26+export interface DriveOptions {
27+ m: number; // driven mode (m, n)
28+ n: number;
29+ E0: number; // starting energy, in units of the nonlinear scale
30+ perturb: boolean; // add a tiny random motion in every other mode
31+}
32+
33+export class DrivenPlate {
34+ P: Params;
35+ g: Grid;
36+ plate: Plate;
37+ modes: Modes;
38+ opts: DriveOptions;
39+ omega: number; // the driven mode's (linear) omega
40+ period: number; // its period
41+ Enl: number; // its nonlinear scale
42+ u: Float64Array;
43+ p: Float64Array;
44+ t = 0; // time
45+ E = 0; // energy now (kinetic + potential)
46+ Et: number; // target energy (absolute units)
47+ Wd = 0; // work done by the drive
48+ Estart: number; // energy at the start, for the energy error
49+ G = 0; // controller gain
50+ driveOn = true; // false: no drive (G = 0)
51+ q = 0; // driven mode's coordinate and velocity
52+ qd = 0;
53+ private S: Float64Array; // driven mode's unit shape
54+ private F: Float64Array;
55+ private V = 0; // potential energy at u
56+ private h = 0;
57+ private sqM: number;
58+ // Measured period of the driven mode: upward zero crossings of q.
59+ private lastCross = NaN;
60+ measuredPeriod = NaN;
61+
62+ constructor(P: Params, opts: DriveOptions) {
63+ this.P = P;
64+ this.g = gridFor(P);
65+ this.plate = new Plate(this.g);
66+ this.modes = new Modes(this.g);
67+ this.opts = opts;
68+ const { m, n } = opts;
69+ this.omega = omegaOf(this.g, m, n);
70+ this.period = (2 * Math.PI) / this.omega;
71+ this.Enl = nonlinearScale(this.plate, m, n);
72+ this.S = modeShape(this.g, m, n);
73+ this.sqM = Math.sqrt(this.g.mass);
74+ this.u = new Float64Array(this.g.n);
75+ this.p = new Float64Array(this.g.n);
76+ this.F = new Float64Array(this.g.n);
77+ // Start in the driven mode at E0 (at rest, displaced), plus the perturbation.
78+ const E0 = opts.E0 * this.Enl, q0 = Math.sqrt(2 * E0) / this.omega;
79+ for (let k = 0; k < this.g.n; k++) this.u[k] = (q0 * this.S[k]) / this.sqM;
80+ if (opts.perturb) this.addPerturbation(PERTURB * E0);
81+ this.Et = E0;
82+ this.V = this.plate.forces(this.u, this.F);
83+ this.E = this.kinetic() + this.V;
84+ this.Estart = this.E;
85+ this.readMode();
86+ this.setStep();
87+ }
88+
89+ // Random phases, equal energy e / (count - 1) in every mode but the driven
90+ // one (linear energies; tiny, so the stretching is negligible). A fixed
91+ // seed, so runs repeat.
92+ private addPerturbation(e: number): void {
93+ const g = this.g, M = this.modes, each = e / (M.count - 1);
94+ let seed = 12345;
95+ const rand = () => ((seed = (seed * 1103515245 + 12345) >>> 0) / 2 ** 32);
96+ for (let k = 0; k < M.count; k++) {
97+ if (M.m[k] === this.opts.m && M.n[k] === this.opts.n) continue;
98+ const S = modeShape(g, M.m[k], M.n[k]), w = M.omega[k], a = Math.sqrt(2 * each), ph = 2 * Math.PI * rand();
99+ const q = (a / w) * Math.cos(ph), qd = a * Math.sin(ph);
100+ for (let i = 0; i < g.n; i++) {
101+ this.u[i] += (q * S[i]) / this.sqM;
102+ this.p[i] += qd * S[i] * this.sqM;
103+ }
104+ }
105+ }
106+
107+ kinetic(): number {
108+ let T = 0;
109+ for (let k = 0; k < this.g.n; k++) T += this.p[k] * this.p[k];
110+ return T / (2 * this.g.mass);
111+ }
112+ // Driven mode's linear energy.
113+ get Ed(): number {
114+ return (this.qd * this.qd + this.omega * this.omega * this.q * this.q) / 2;
115+ }
116+ private readMode(): void {
117+ let a = 0, b = 0;
118+ for (let k = 0; k < this.g.n; k++) {
119+ a += this.S[k] * this.u[k];
120+ b += this.S[k] * this.p[k];
121+ }
122+ this.q = a * this.sqM;
123+ this.qd = b / this.sqM;
124+ }
125+ private setStep(): void {
126+ this.h = stepFor(this.g, this.modes, Math.max(this.E, this.Et));
127+ }
128+ // Controller gain from the energy now and the target.
129+ private setGain(): void {
130+ const Ed = Math.max(this.Ed, 1e-3 * this.E, 1e-300);
131+ const G = (this.Et - this.E) / (TAU_C * this.period * Ed), cap = G_MAX * this.omega;
132+ this.G = this.driveOn ? Math.max(-cap, Math.min(cap, G)) : 0;
133+ }
134+
135+ // Advance by about dt (whole steps); returns the steps taken.
136+ advance(dt: number): number {
137+ this.setStep();
138+ const steps = Math.max(1, Math.round(dt / this.h)), h = dt / steps;
139+ const { u, p, F, S } = this, n = this.g.n, M = this.g.mass, sqM = this.sqM;
140+ // Gain update every control interval (a few per driven period).
141+ const every = Math.max(1, Math.round(this.period / 32 / h));
142+ for (let s = 0; s < steps; s++) {
143+ if (s % every === 0) {
144+ this.E = this.kinetic() + this.V;
145+ this.setGain();
146+ }
147+ const Q = this.G * this.qd, f = Q * sqM, qBefore = this.q;
148+ for (let k = 0; k < n; k++) p[k] += 0.5 * h * (F[k] + f * S[k]);
149+ for (let k = 0; k < n; k++) u[k] += (h * p[k]) / M;
150+ this.V = this.plate.forces(u, F);
151+ for (let k = 0; k < n; k++) p[k] += 0.5 * h * (F[k] + f * S[k]);
152+ this.readMode();
153+ this.Wd += Q * (this.q - qBefore);
154+ if (qBefore < 0 && this.q >= 0) {
155+ const tc = this.t + (h * -qBefore) / (this.q - qBefore);
156+ if (Number.isFinite(this.lastCross)) this.measuredPeriod = tc - this.lastCross;
157+ this.lastCross = tc;
158+ }
159+ this.t += h;
160+ }
161+ this.E = this.kinetic() + this.V;
162+ return steps;
163+ }
164+
165+ // Relative energy error: (E - Estart - Wd) / E.
166+ get energyError(): number {
167+ return (this.E - this.Estart - this.Wd) / this.E;
168+ }
169+ get step(): number {
170+ return this.h;
171+ }
172+ // Linear energy of every mode, (qd^2 + omega^2 q^2) / 2, into out at index
173+ // (m - 1) + Nx (n - 1), from the sine transforms of u and p.
174+ private Uhat: Float64Array | null = null;
175+ private Phat: Float64Array | null = null;
176+ private omegaGrid: Float64Array | null = null;
177+ modalEnergies(out: Float64Array): void {
178+ const g = this.g, n = g.n;
179+ const U = (this.Uhat ??= new Float64Array(n)), Pp = (this.Phat ??= new Float64Array(n));
180+ if (!this.omegaGrid) {
181+ this.omegaGrid = new Float64Array(n);
182+ for (let b = 0; b < g.Ny; b++) for (let a = 0; a < g.Nx; a++) this.omegaGrid[a + g.Nx * b] = omegaOf(g, a + 1, b + 1);
183+ }
184+ this.plate.sine2(this.u, U);
185+ this.plate.sine2(this.p, Pp);
186+ const W = this.omegaGrid, M = g.mass;
187+ for (let k = 0; k < n; k++) out[k] = (Pp[k] * Pp[k]) / (2 * M) + (W[k] * W[k] * M * U[k] * U[k]) / 2;
188+ }
189+ // Index of the driven mode in modalEnergies' output.
190+ get drivenIndex(): number {
191+ return this.opts.m - 1 + this.g.Nx * (this.opts.n - 1);
192+ }
193+ // Largest |u| (in thicknesses).
194+ maxAbsU(): number {
195+ let a = 0;
196+ for (let k = 0; k < this.g.n; k++) a = Math.max(a, Math.abs(this.u[k]));
197+ return a;
198+ }
199+}
src/driven-plate/index.htmladded+82−0View file
@@ -0,0 +1,82 @@
1+<!doctype html>
2+<html lang="en">
3+ <head>
4+ <meta charset="UTF-8" />
5+ <meta name="viewport" content="width=device-width, initial-scale=1.0" />
6+ <title>Driven plate</title>
7+ </head>
8+ <body>
9+ <main>
10+ <header><h1>Driven plate</h1></header>
11+ <div class="grid">
12+ <section class="card">
13+ <canvas id="plate"></canvas>
14+ <div id="readout" class="num"></div>
15+ <div id="controls">
16+ <div class="row">
17+ <button id="run" class="primary">Pause</button>
18+ <button id="reset">Reset</button>
19+ <label class="chip" title="Raise the target energy steadily (log scale) up to the top of the slider"><input id="ramp" type="checkbox" checked /><span>Ramp</span></label>
20+ <span class="muted"><input id="per-decade" type="number" min="10" step="10" value="400" /> periods per decade</span>
21+ <button type="button" class="help" data-help="drive" title="About the driven plate" aria-label="About the driven plate">?</button>
22+ </div>
23+ <div class="row">
24+ <span class="muted label-w">Target E / E<sub>nl</sub></span>
25+ <input id="target" type="range" min="-4" max="2" step="0.01" value="-3" />
26+ <span id="target-value" class="num value-w"></span>
27+ </div>
28+ <div class="row">
29+ <span class="muted label-w">Speed</span>
30+ <input id="speed" type="range" min="0" max="2.5" step="0.01" value="1.3" />
31+ <span id="speed-value" class="num value-w"></span>
32+ </div>
33+ </div>
34+ </section>
35+ <div class="column">
36+ <section class="card">
37+ <div class="card-head">Plate and drive <button type="button" class="help" data-help="params" title="About the parameters" aria-label="About the parameters">?</button></div>
38+ <div class="row">
39+ <span class="muted">Masses across</span>
40+ <input id="ny" type="number" min="2" max="32" step="1" value="8" />
41+ <span class="muted" style="margin-left: 8px">Aspect</span>
42+ <input id="aspect" type="number" min="1" max="3" step="0.05" value="1.3" />
43+ <span id="grid-info" class="muted"></span>
44+ </div>
45+ <div class="row">
46+ <span class="muted">Driven mode</span>
47+ <select id="mode"></select>
48+ <label class="chip" title="Start with a tiny random motion (a millionth of the starting energy) in every other mode; off, the plate starts in the driven mode exactly"><input id="perturb" type="checkbox" checked /><span>Perturbation</span></label>
49+ </div>
50+ <div id="params"></div>
51+ </section>
52+ <section class="card">
53+ <div class="card-head">Energy <button type="button" class="help" data-help="energy" title="About the energy chart" aria-label="About the energy chart">?</button></div>
54+ <canvas id="chart"></canvas>
55+ </section>
56+ </div>
57+ </div>
58+ </main>
59+ <dialog id="help">
60+ <div id="help-body"></div>
61+ <div id="help-actions"><button id="help-close" type="button">Close</button></div>
62+ </dialog>
63+
64+ <template id="help-drive">
65+ <h2>The driven plate</h2>
66+ <p>A thin rectangular plate, simply supported at its edges, seen from above: red is up, blue down, white at rest (scaled to the recent largest displacement). It is a grid of equal masses with bending stiffness and von Kármán stretching: a bent plate stretches in its own plane, which stiffens it once it moves by about its thickness h. There is no damping.</p>
67+ <p>A force in the shape of one normal mode (the driven mode) pushes the plate. In a linear plate such a force moves that mode and nothing else, so any other mode that moves got its energy through the stretching. The force is always in phase with the driven mode's own velocity, so it adds energy whatever the mode's frequency (which rises as the plate stiffens), and a controller sets its strength so that the plate's total energy follows the target.</p>
68+ <p>The target is given against E<sub>nl</sub>, the energy of the driven mode at which its stretching energy equals its bending energy: well below 1 the plate is nearly linear. Ramp raises the target steadily, the same number of periods for each factor of ten, up to the top of the slider; untick it to hold the energy where it is. Reset starts again in the driven mode at 10<sup>−3</sup> E<sub>nl</sub>. Time is counted in periods of the driven mode (linear), and Speed sets how many pass per second; at high energies the plate needs many more steps per period, and the actual speed (shown) can fall below the setting. The space bar pauses.</p>
69+ </template>
70+ <template id="help-params">
71+ <h2>Plate and drive</h2>
72+ <p>Masses across is the number of masses across the short side; the long side gets as many as keep the cells nearly square. Aspect is the ratio of the sides. A square plate has pairs of modes at the same pitch ((1, 2) and (2, 1)), which trade energy easily; a non-square plate splits them. The default 1.3 keeps the lowest modes at least 1% apart in pitch.</p>
73+ <p>The modes are products of sines, numbered (m, n) by their half waves along x and y. The list gives each one's pitch relative to the lowest. Perturbation starts the plate with a tiny random motion in every other mode, a millionth of the starting energy in all; without it the plate starts in the driven mode exactly, and by symmetry only modes with the driven mode's symmetry can then grow (for (1, 1): odd m and odd n), at first; round-off seeds the rest. Any change here starts again.</p>
74+ <p>Units are scaled away: for the von Kármán plate only the aspect ratio and Poisson's ratio (0.3) matter once energy is given against E<sub>nl</sub>, displacement against the thickness h and time in periods.</p>
75+ </template>
76+ <template id="help-energy">
77+ <h2>Energy</h2>
78+ <p>Against time in driven periods, on a log scale in units of E<sub>nl</sub>: the total energy (black), the target (dashed), the driven mode's energy (blue) and the sum over all other modes (orange). The mode energies are the linear ones, (velocity² + ω² displacement²) / 2 per mode; once stretching matters they no longer add up exactly to the total.</p>
79+ </template>
80+ <script type="module" src="./main.ts"></script>
81+ </body>
82+</html>
src/driven-plate/main.tsadded+298−0View file
@@ -0,0 +1,298 @@
1+import "../common/style.css";
2+import "./style.css";
3+import { fit } from "../common/canvas";
4+import { setupHelp } from "../common/help";
5+import { DrivenPlate } from "./drive";
6+import { Modes, gridFor, type Params } from "./plate";
7+
8+setupHelp();
9+
10+const LOG_START = -3; // log10 of the starting energy, in E_nl
11+const BUDGET_MS = 14; // compute per frame at most
12+
13+const $ = <T extends HTMLElement>(id: string) => document.getElementById(id) as T;
14+const plateCanvas = $<HTMLCanvasElement>("plate");
15+const chartCanvas = $<HTMLCanvasElement>("chart");
16+const readout = $("readout");
17+const runButton = $<HTMLButtonElement>("run");
18+const rampBox = $<HTMLInputElement>("ramp");
19+const perDecadeInput = $<HTMLInputElement>("per-decade");
20+const targetSlider = $<HTMLInputElement>("target");
21+const targetValue = $("target-value");
22+const speedSlider = $<HTMLInputElement>("speed");
23+const speedValue = $("speed-value");
24+const nyInput = $<HTMLInputElement>("ny");
25+const aspectInput = $<HTMLInputElement>("aspect");
26+const gridInfo = $("grid-info");
27+const modeSelect = $<HTMLSelectElement>("mode");
28+const perturbBox = $<HTMLInputElement>("perturb");
29+const paramsDiv = $("params");
30+
31+let P: Params = { Ny: 8, aspect: 1.3 };
32+let sim: DrivenPlate;
33+let running = true;
34+let logTarget = LOG_START;
35+let actualSpeed = 0; // periods per second, smoothed
36+let ampl = 0; // colour scale: recent largest |u|
37+let modal = new Float64Array(0);
38+// History for the chart: time (periods) and log10 of E, target, driven, others (in E_nl).
39+let hist: { t: number[]; E: number[]; Et: number[]; Ed: number[]; Eo: number[] };
40+
41+const fmt = (x: number, d = 3) => (x === 0 ? "0" : Math.abs(x) >= 1e-2 && Math.abs(x) < 1e4 ? x.toPrecision(d) : x.toExponential(d - 1));
42+const speed = () => 10 ** Number(speedSlider.value);
43+
44+// The mode list for the current grid: lowest modes by pitch.
45+function fillModes(): void {
46+ const g = gridFor(P), M = new Modes(g), keep = modeSelect.value || "1,1";
47+ modeSelect.replaceChildren();
48+ for (let k = 0; k < Math.min(M.count, 20); k++) {
49+ const o = document.createElement("option");
50+ o.value = `${M.m[k]},${M.n[k]}`;
51+ o.textContent = `(${M.m[k]}, ${M.n[k]}) ${(M.omega[k] / M.omega[0]).toFixed(2)}`;
52+ modeSelect.append(o);
53+ }
54+ modeSelect.value = keep;
55+ if (modeSelect.selectedIndex < 0) modeSelect.value = "1,1";
56+ gridInfo.textContent = `${g.Nx} × ${g.Ny} masses`;
57+}
58+
59+function reset(): void {
60+ const [m, n] = modeSelect.value.split(",").map(Number);
61+ sim = new DrivenPlate(P, { m, n, E0: 10 ** LOG_START, perturb: perturbBox.checked });
62+ logTarget = LOG_START;
63+ targetSlider.value = String(LOG_START);
64+ modal = new Float64Array(sim.g.n);
65+ hist = { t: [], E: [], Et: [], Ed: [], Eo: [] };
66+ ampl = 0;
67+ record();
68+ showParams();
69+}
70+
71+function showParams(): void {
72+ const g = sim.g, M = sim.modes;
73+ const lines = [
74+ `grid ${g.Nx} × ${g.Ny}, sides ${g.Lx.toFixed(2)} × ${g.Ly.toFixed(2)} (short side 1), ν 0.3`,
75+ `driven (${sim.opts.m}, ${sim.opts.n}), pitch ${(sim.omega / M.omega[0]).toFixed(3)} × lowest`,
76+ `top mode ${(M.omegaMax / sim.omega).toFixed(1)} × driven`,
77+ `lowest: ${Array.from(M.omega.slice(0, 8), (w, k) => `(${M.m[k]},${M.n[k]}) ${(w / M.omega[0]).toFixed(2)}`).join(" ")}`,
78+ ];
79+ paramsDiv.textContent = lines.join("\n");
80+}
81+
82+function record(): void {
83+ sim.modalEnergies(modal);
84+ let others = 0;
85+ const d = sim.drivenIndex;
86+ for (let k = 0; k < modal.length; k++) if (k !== d) others += modal[k];
87+ const L = (x: number) => Math.log10(Math.max(x / sim.Enl, 1e-30));
88+ hist.t.push(sim.t / sim.period);
89+ hist.E.push(L(sim.E));
90+ hist.Et.push(L(sim.Et));
91+ hist.Ed.push(L(modal[d]));
92+ hist.Eo.push(L(others));
93+ // Keep the history short: drop every other point when it grows long.
94+ if (hist.t.length > 4000) for (const key of ["t", "E", "Et", "Ed", "Eo"] as const) hist[key] = hist[key].filter((_, i) => i % 2 === 0);
95+}
96+
97+function setTarget(): void {
98+ sim.Et = 10 ** logTarget * sim.Enl;
99+ targetValue.textContent = fmt(10 ** logTarget);
100+}
101+
102+// One frame: advance by speed x frame time (in driven periods), within the
103+// compute budget.
104+let last = performance.now();
105+let msPerPeriod = 1; // compute cost, measured
106+function frame(now: number): void {
107+ const wall = Math.min(0.1, (now - last) / 1000);
108+ last = now;
109+ if (running) {
110+ let periods = speed() * wall;
111+ periods = Math.min(periods, BUDGET_MS / msPerPeriod);
112+ if (rampBox.checked) {
113+ const top = Number(targetSlider.max);
114+ logTarget = Math.min(top, logTarget + periods / Number(perDecadeInput.value || 400));
115+ targetSlider.value = String(logTarget);
116+ }
117+ setTarget();
118+ if (periods > 0) {
119+ const t0 = performance.now();
120+ sim.advance(periods * sim.period);
121+ const ms = performance.now() - t0;
122+ msPerPeriod = 0.8 * msPerPeriod + 0.2 * Math.max(ms / periods, 1e-3);
123+ }
124+ actualSpeed = 0.9 * actualSpeed + 0.1 * (wall > 0 ? periods / wall : 0);
125+ record();
126+ } else actualSpeed = 0;
127+ draw();
128+ requestAnimationFrame(frame);
129+}
130+
131+function draw(): void {
132+ drawPlate();
133+ drawChart();
134+ speedValue.textContent = `${fmt(speed(), 2)} /s`;
135+ const others = 10 ** hist.Eo[hist.Eo.length - 1], driven = 10 ** hist.Ed[hist.Ed.length - 1], E = sim.E / sim.Enl;
136+ const T = sim.measuredPeriod / sim.period;
137+ readout.textContent =
138+ `t ${(sim.t / sim.period).toFixed(0)} periods · E ${fmt(E)} E_nl (target ${fmt(sim.Et / sim.Enl)}) · driven mode ${fmt((100 * driven) / E, 3)}% · other modes ${fmt((100 * others) / E, 2)}%\n` +
139+ `driven pitch ${Number.isFinite(T) ? fmt(1 / T) : "–"} × linear · max |w| ${fmt(sim.maxAbsU(), 2)} h · ` +
140+ `speed ${running ? fmt(actualSpeed, 2) : "0"} periods/s · energy error ${Math.abs(sim.energyError).toExponential(1)}`;
141+}
142+
143+// Height map, as in ../plate: u bilinear between the nodes (0 at the
144+// supports), red up, blue down, scaled to the recent largest |u|.
145+const RES = 200;
146+const heat = document.createElement("canvas");
147+const heatCtx = heat.getContext("2d")!;
148+let heatImage: ImageData | null = null;
149+function drawPlate(): void {
150+ const ctx = fit(plateCanvas), W = plateCanvas.clientWidth, H = plateCanvas.clientHeight;
151+ ctx.clearRect(0, 0, W, H);
152+ const g = sim.g, u = sim.u;
153+ const mx = sim.maxAbsU();
154+ ampl = Math.max(mx, ampl * 0.97);
155+ const scale = ampl > 0 ? 1 / ampl : 0;
156+ const rw = RES, rh = Math.max(2, Math.round((RES * g.Ly) / g.Lx));
157+ if (!heatImage || heatImage.width !== rw || heatImage.height !== rh) {
158+ heat.width = rw;
159+ heat.height = rh;
160+ heatImage = heatCtx.createImageData(rw, rh);
161+ }
162+ const d = heatImage.data;
163+ const at = (i: number, j: number) => (i < 0 || j < 0 || i >= g.Nx || j >= g.Ny ? 0 : u[i + g.Nx * j]);
164+ for (let py = 0; py < rh; py++) {
165+ const fy = (((rh - 1 - py + 0.5) / rh) * g.Ly) / g.ay - 1; // y up
166+ const j0 = Math.floor(fy), by = fy - j0;
167+ for (let px = 0; px < rw; px++) {
168+ const fx = (((px + 0.5) / rw) * g.Lx) / g.ax - 1, i0 = Math.floor(fx), bx = fx - i0;
169+ const v0 = (1 - bx) * at(i0, j0) + bx * at(i0 + 1, j0), v1 = (1 - bx) * at(i0, j0 + 1) + bx * at(i0 + 1, j0 + 1);
170+ const v = Math.max(-1, Math.min(1, ((1 - by) * v0 + by * v1) * scale)), q = 4 * (px + rw * py);
171+ if (v >= 0) {
172+ d[q] = 255;
173+ d[q + 1] = 255 - 200 * v;
174+ d[q + 2] = 255 - 210 * v;
175+ } else {
176+ d[q] = 255 + 210 * v;
177+ d[q + 1] = 255 + 150 * v;
178+ d[q + 2] = 255 + 40 * v;
179+ }
180+ d[q + 3] = 255;
181+ }
182+ }
183+ heatCtx.putImageData(heatImage, 0, 0);
184+ const pad = 12, top = 8, bottom = 22;
185+ const s = Math.min((W - 2 * pad) / g.Lx, (H - top - bottom) / g.Ly);
186+ const w = s * g.Lx, h = s * g.Ly, x0 = (W - w) / 2, y0 = top;
187+ ctx.imageSmoothingEnabled = true;
188+ ctx.drawImage(heat, x0, y0, w, h);
189+ ctx.strokeStyle = "#7a7d82";
190+ ctx.lineWidth = 2;
191+ ctx.strokeRect(x0, y0, w, h);
192+ ctx.fillStyle = "#7a7d82";
193+ ctx.font = "11px system-ui, sans-serif";
194+ ctx.textAlign = "right";
195+ ctx.fillText(ampl > 0 ? `colour full at ±${fmt(ampl, 2)} h` : "", x0 + w, H - 6);
196+ ctx.textAlign = "left";
197+}
198+
199+// Energy chart: log10(E / E_nl) against time in periods.
200+const SERIES = [
201+ { key: "E", color: "#23262b", dash: [] as number[], label: "total" },
202+ { key: "Et", color: "#7a7d82", dash: [4, 4], label: "target" },
203+ { key: "Ed", color: "#2f5d8a", dash: [], label: "driven mode" },
204+ { key: "Eo", color: "#c8742c", dash: [], label: "other modes" },
205+] as const;
206+function drawChart(): void {
207+ const ctx = fit(chartCanvas), W = chartCanvas.clientWidth, H = chartCanvas.clientHeight;
208+ ctx.clearRect(0, 0, W, H);
209+ const n = hist.t.length;
210+ if (n < 2) return;
211+ const l = 40, r = 8, t = 22, b = 22, pw = W - l - r, ph = H - t - b;
212+ const tMax = Math.max(hist.t[n - 1], 1);
213+ let lo = Infinity, hi = -Infinity;
214+ for (const s of SERIES) for (const v of hist[s.key]) if (v > -30) {
215+ lo = Math.min(lo, v);
216+ hi = Math.max(hi, v);
217+ }
218+ lo = Math.max(Math.floor(lo), hi - 14);
219+ hi = Math.ceil(hi + 0.01);
220+ const X = (x: number) => l + (x / tMax) * pw, Y = (y: number) => t + ((hi - Math.max(lo, y)) / (hi - lo)) * ph;
221+ ctx.font = "11px system-ui, sans-serif";
222+ ctx.fillStyle = "#7a7d82";
223+ ctx.strokeStyle = "#e3e2dc";
224+ ctx.lineWidth = 1;
225+ const stepY = hi - lo > 8 ? 2 : 1;
226+ ctx.textAlign = "right";
227+ for (let y = Math.ceil(lo); y <= hi; y += stepY) {
228+ ctx.beginPath();
229+ ctx.moveTo(l, Y(y));
230+ ctx.lineTo(l + pw, Y(y));
231+ ctx.stroke();
232+ ctx.fillText(`1e${y}`, l - 4, Y(y) + 4);
233+ }
234+ ctx.textAlign = "center";
235+ ctx.fillText(`time (driven periods), 0 to ${tMax.toFixed(0)}`, l + pw / 2, H - 6);
236+ for (const s of SERIES) {
237+ const v = hist[s.key];
238+ ctx.strokeStyle = s.color;
239+ ctx.setLineDash(s.dash);
240+ ctx.lineWidth = 1.5;
241+ ctx.beginPath();
242+ for (let i = 0; i < n; i++) (i ? ctx.lineTo : ctx.moveTo).call(ctx, X(hist.t[i]), Y(v[i]));
243+ ctx.stroke();
244+ }
245+ ctx.setLineDash([]);
246+ // Legend.
247+ ctx.textAlign = "left";
248+ let x = l;
249+ for (const s of SERIES) {
250+ ctx.strokeStyle = s.color;
251+ ctx.setLineDash(s.dash);
252+ ctx.beginPath();
253+ ctx.moveTo(x, 10);
254+ ctx.lineTo(x + 16, 10);
255+ ctx.stroke();
256+ ctx.setLineDash([]);
257+ ctx.fillStyle = "#23262b";
258+ ctx.fillText(s.label, x + 20, 14);
259+ x += 28 + ctx.measureText(s.label).width;
260+ }
261+ ctx.fillStyle = "#7a7d82";
262+ ctx.textAlign = "right";
263+ ctx.fillText("E / E_nl", W - r, 14);
264+}
265+
266+// Controls.
267+function setRunning(on: boolean): void {
268+ running = on;
269+ runButton.textContent = on ? "Pause" : "Run";
270+}
271+runButton.addEventListener("click", () => setRunning(!running));
272+$("reset").addEventListener("click", reset);
273+document.addEventListener("keydown", (e) => {
274+ if (e.code === "Space" && !(e.target instanceof HTMLInputElement || e.target instanceof HTMLSelectElement || e.target instanceof HTMLButtonElement)) {
275+ e.preventDefault();
276+ setRunning(!running);
277+ }
278+});
279+targetSlider.addEventListener("input", () => {
280+ logTarget = Number(targetSlider.value);
281+ setTarget();
282+});
283+function applyGrid(): void {
284+ const Ny = Math.round(Number(nyInput.value)), aspect = Number(aspectInput.value);
285+ if (!(Ny >= 2 && Ny <= 32 && aspect >= 1 && aspect <= 3)) return;
286+ P = { Ny, aspect };
287+ fillModes();
288+ reset();
289+}
290+nyInput.addEventListener("change", applyGrid);
291+aspectInput.addEventListener("change", applyGrid);
292+modeSelect.addEventListener("change", reset);
293+perturbBox.addEventListener("change", reset);
294+
295+fillModes();
296+reset();
297+setTarget();
298+requestAnimationFrame(frame);
src/driven-plate/plate.tsadded+252−0View file
@@ -0,0 +1,252 @@
1+// A rectangular plate, simply supported, as Nx x Ny equal masses on vertical
2+// rails at the interior nodes of a grid, moving up and down only, with von
3+// Karman stretching. Node (i, j) at x = (i + 1) ax, y = (j + 1) ay, index i +
4+// Nx j; the edge nodes (i = -1 or Nx, j = -1 or Ny) are supports at u = 0.
5+// The sides are Lx = (Nx + 1) ax and Ly = (Ny + 1) ay.
6+// Units: the von Karman plate, scaled, has only the aspect ratio and Poisson's
7+// ratio left as parameters (lengths in the plane by Ly, displacement by the
8+// thickness h, time by Ly^2 sqrt(rho h / D)). So D = 1, rho h = 1, h = 1, Ly
9+// = 1 and Lx = aspect; energies are reported against the nonlinear scale and
10+// times in periods of the driven mode, so these units never show.
11+// Bending: D / 2 times the area-weighted sum of (Lap u)^2, Lap the five-point
12+// Laplacian with spacings ax, ay (L u = -Lap u below, positive); a simply
13+// supported edge has no term (mirror ghost). Force -D ax ay L (L u). The modes
14+// are exact products of sines, S_mn(i, j) = sx_m(i) sy_n(j) with sx_m(i) =
15+// sqrt(2 / (Nx + 1)) sin(m pi (i + 1) / (Nx + 1)), and L S_mn = Lam_mn S_mn,
16+// Lam_mn = lx_m / ax^2 + ly_n / ay^2, l_m = 4 sin^2(m pi / 2 (N + 1)), so
17+// omega_mn = Lam_mn sqrt(D / rho h): for low modes (m pi / Lx)^2 + (n pi /
18+// Ly)^2, the continuum plate.
19+// Stretching (von Karman), as in ../plate/plate.ts with two spacings: s =
20+// l(u, u) = 2 (X Y - Z^2) at each node, X = Dxx u / ax^2, Y = Dyy u / ay^2, Z
21+// = Dxy u / (ax ay) (second differences, zero beyond the edge, Dxy the centred
22+// cross difference), phi = L^-2 s by sine transforms (F = Lap F = 0 at the
23+// edge), energy (Eh / 8) ax ay sum s phi, Eh = 12 (1 - nu^2) D / h^2. The
24+// force is its exact gradient, -(Eh / 2) ax ay (Dxx (phi Y) / ax^2 + Dyy (phi
25+// X) / ay^2 - 2 Dxy (phi Z) / (ax ay)), so velocity Verlet conserves energy.
26+// Loss (none yet) and the drive live in drive.ts.
27+
28+export const NU = 0.3;
29+const D = 1, RHO_H = 1, THICK = 1;
30+const EH = (12 * (1 - NU * NU) * D) / (THICK * THICK);
31+
32+export interface Params {
33+ Ny: number; // masses across the short side (y)
34+ aspect: number; // Lx / Ly, at least 1
35+}
36+
37+// The grid for the given parameters: Nx chosen so the cells are as near
38+// square as possible.
39+export interface Grid {
40+ Nx: number;
41+ Ny: number;
42+ n: number; // Nx Ny
43+ Lx: number;
44+ Ly: number;
45+ ax: number;
46+ ay: number;
47+ mass: number; // of each node, rho h ax ay
48+}
49+export function gridFor(P: Params): Grid {
50+ const Ly = 1, Lx = P.aspect, Ny = P.Ny;
51+ const Nx = Math.max(1, Math.round(Lx * (Ny + 1)) - 1);
52+ const ax = Lx / (Nx + 1), ay = Ly / (Ny + 1);
53+ return { Nx, Ny, n: Nx * Ny, Lx, Ly, ax, ay, mass: RHO_H * ax * ay };
54+}
55+
56+const lam1 = (N: number, m: number) => 4 * Math.sin((m * Math.PI) / (2 * (N + 1))) ** 2;
57+// Orthonormal sine matrix for N points, S[m * N + i], m = 0..N-1 (mode m + 1).
58+function sineMatrix(N: number): Float64Array {
59+ const S = new Float64Array(N * N);
60+ for (let m = 0; m < N; m++) for (let i = 0; i < N; i++) S[m * N + i] = Math.sqrt(2 / (N + 1)) * Math.sin(((m + 1) * (i + 1) * Math.PI) / (N + 1));
61+ return S;
62+}
63+
64+// The normal modes (exact), sorted by frequency: omega (rad / time unit) and
65+// (m, n), ties by m then n. Mode (m, n) has Lam = omega (D = rho h = 1).
66+export class Modes {
67+ count: number;
68+ m: Int32Array;
69+ n: Int32Array;
70+ omega: Float64Array;
71+ constructor(g: Grid) {
72+ const all: [number, number, number][] = [];
73+ for (let m = 1; m <= g.Nx; m++) for (let n = 1; n <= g.Ny; n++) all.push([omegaOf(g, m, n), m, n]);
74+ all.sort((x, y) => x[0] - y[0] || x[1] - y[1] || x[2] - y[2]);
75+ this.count = all.length;
76+ this.omega = Float64Array.from(all, (x) => x[0]);
77+ this.m = Int32Array.from(all, (x) => x[1]);
78+ this.n = Int32Array.from(all, (x) => x[2]);
79+ }
80+ get omegaMax(): number {
81+ return this.omega[this.count - 1];
82+ }
83+}
84+export const omegaOf = (g: Grid, m: number, n: number) => (lam1(g.Nx, m) / g.ax ** 2 + lam1(g.Ny, n) / g.ay ** 2) * Math.sqrt(D / RHO_H);
85+
86+// Unit-length shape of mode (m, n) on the nodes (sum of squares 1).
87+export function modeShape(g: Grid, m: number, n: number): Float64Array {
88+ const v = new Float64Array(g.n), cx = Math.sqrt(2 / (g.Nx + 1)), cy = Math.sqrt(2 / (g.Ny + 1));
89+ for (let j = 0; j < g.Ny; j++) {
90+ const sy = cy * Math.sin((n * (j + 1) * Math.PI) / (g.Ny + 1));
91+ for (let i = 0; i < g.Nx; i++) v[i + g.Nx * j] = sy * cx * Math.sin((m * (i + 1) * Math.PI) / (g.Nx + 1));
92+ }
93+ return v;
94+}
95+
96+// Forces and potential energy. One per grid; holds the scratch arrays.
97+export class Plate {
98+ g: Grid;
99+ private Sx: Float64Array;
100+ private Sy: Float64Array;
101+ private inv2: Float64Array; // 1 / Lam_mn^2, at m + Nx n (0-based)
102+ private c: Float64Array;
103+ private xx: Float64Array;
104+ private yy: Float64Array;
105+ private xy: Float64Array;
106+ private s: Float64Array;
107+ private phi: Float64Array;
108+ private t: Float64Array;
109+ private t2: Float64Array;
110+ private pa: Float64Array; // padded, (Nx + 2)(Ny + 2)
111+ private pb: Float64Array;
112+ private pc: Float64Array;
113+ constructor(g: Grid) {
114+ this.g = g;
115+ const { Nx, Ny, n } = g, f = () => new Float64Array(n), pad = () => new Float64Array((Nx + 2) * (Ny + 2));
116+ this.Sx = sineMatrix(Nx);
117+ this.Sy = sineMatrix(Ny);
118+ this.inv2 = f();
119+ for (let q = 0; q < Ny; q++) for (let p = 0; p < Nx; p++) this.inv2[p + Nx * q] = 1 / (lam1(Nx, p + 1) / g.ax ** 2 + lam1(Ny, q + 1) / g.ay ** 2) ** 2;
120+ this.c = f();
121+ this.xx = f();
122+ this.yy = f();
123+ this.xy = f();
124+ this.s = f();
125+ this.phi = f();
126+ this.t = f();
127+ this.t2 = f();
128+ this.pa = pad();
129+ this.pb = pad();
130+ this.pc = pad();
131+ }
132+
133+ // out = (Sy x Sx) X: the sine transform along x then y (its own inverse).
134+ sine2(X: Float64Array, out: Float64Array): void {
135+ const { Nx, Ny } = this.g, Sx = this.Sx, Sy = this.Sy, t = this.t;
136+ for (let j = 0; j < Ny; j++) {
137+ const r = Nx * j;
138+ for (let m = 0; m < Nx; m++) {
139+ let acc = 0;
140+ for (let i = 0; i < Nx; i++) acc += Sx[m * Nx + i] * X[r + i];
141+ t[r + m] = acc;
142+ }
143+ }
144+ for (let m = 0; m < Nx; m++) {
145+ for (let q = 0; q < Ny; q++) {
146+ let acc = 0;
147+ for (let j = 0; j < Ny; j++) acc += Sy[q * Ny + j] * t[m + Nx * j];
148+ out[m + Nx * q] = acc;
149+ }
150+ }
151+ }
152+
153+ // L x into out (positive: minus the Laplacian, x = 0 at the supports).
154+ private laplacian(x: Float64Array, out: Float64Array): void {
155+ const { Nx, Ny, ax, ay } = this.g, bx = 1 / (ax * ax), by = 1 / (ay * ay);
156+ for (let j = 0, k = 0; j < Ny; j++) {
157+ for (let i = 0; i < Nx; i++, k++) {
158+ let sx = 2 * x[k], sy = 2 * x[k];
159+ if (i > 0) sx -= x[k - 1];
160+ if (i < Nx - 1) sx -= x[k + 1];
161+ if (j > 0) sy -= x[k - Nx];
162+ if (j < Ny - 1) sy -= x[k + Nx];
163+ out[k] = bx * sx + by * sy;
164+ }
165+ }
166+ }
167+
168+ // Forces at u into F; returns the potential energy (bending + stretching).
169+ forces(u: Float64Array, F: Float64Array, linear = false): number {
170+ const g = this.g, n = g.n, A = g.ax * g.ay, c = this.c;
171+ this.laplacian(u, c);
172+ let Vb = 0;
173+ for (let k = 0; k < n; k++) Vb += c[k] * c[k];
174+ this.laplacian(c, F);
175+ for (let k = 0; k < n; k++) F[k] *= -D * A;
176+ if (linear) return (D * A * Vb) / 2;
177+ return (D * A * Vb) / 2 + this.stretching(u, F);
178+ }
179+
180+ // Stretching energy at u; with F given, its force is added.
181+ stretching(u: Float64Array, F?: Float64Array): number {
182+ const { Nx, Ny, n, ax, ay } = this.g, R = Nx + 2;
183+ const { pa: pad, xx, yy, xy, s, phi, t2 } = this;
184+ const bx = 1 / (ax * ax), by = 1 / (ay * ay), bxy = 1 / (ax * ay);
185+ for (let j = 0; j < Ny; j++) for (let i = 0; i < Nx; i++) pad[i + 1 + R * (j + 1)] = u[i + Nx * j];
186+ for (let j = 0, k = 0; j < Ny; j++) {
187+ for (let i = 0, K = 1 + R * (j + 1); i < Nx; i++, k++, K++) {
188+ const cc = 2 * pad[K];
189+ xx[k] = bx * (pad[K - 1] - cc + pad[K + 1]);
190+ yy[k] = by * (pad[K - R] - cc + pad[K + R]);
191+ xy[k] = 0.25 * bxy * (pad[K + R + 1] - pad[K - R + 1] - pad[K + R - 1] + pad[K - R - 1]);
192+ s[k] = 2 * (xx[k] * yy[k] - xy[k] * xy[k]);
193+ }
194+ }
195+ this.sine2(s, t2);
196+ for (let k = 0; k < n; k++) t2[k] *= this.inv2[k];
197+ this.sine2(t2, phi);
198+ let V = 0;
199+ for (let k = 0; k < n; k++) V += s[k] * phi[k];
200+ const A = ax * ay;
201+ if (F) {
202+ // F -= (Eh / 2) A (Dxx (phi Y) bx + Dyy (phi X) by - 2 Dxy (phi Z) bxy);
203+ // the difference operators are symmetric.
204+ // a reuses the padded copy of u (no longer needed); the rings of zeros
205+ // are never written.
206+ const a = this.pa, b = this.pb, c = this.pc;
207+ for (let j = 0, k = 0; j < Ny; j++) {
208+ for (let i = 0, K = 1 + R * (j + 1); i < Nx; i++, k++, K++) {
209+ a[K] = phi[k] * yy[k];
210+ b[K] = phi[k] * xx[k];
211+ c[K] = phi[k] * xy[k];
212+ }
213+ }
214+ const h = 0.5 * EH * A;
215+ for (let j = 0, k = 0; j < Ny; j++) {
216+ for (let i = 0, K = 1 + R * (j + 1); i < Nx; i++, k++, K++) {
217+ const gk =
218+ bx * (a[K - 1] - 2 * a[K] + a[K + 1]) +
219+ by * (b[K - R] - 2 * b[K] + b[K + R]) -
220+ 0.5 * bxy * (c[K + R + 1] - c[K - R + 1] - c[K + R - 1] + c[K - R - 1]);
221+ F[k] -= h * gk;
222+ }
223+ }
224+ }
225+ return (EH / 8) * A * V;
226+ }
227+
228+ private Fs: Float64Array | null = null;
229+ potential(u: Float64Array): number {
230+ return this.forces(u, (this.Fs ??= new Float64Array(this.g.n)));
231+ }
232+}
233+
234+// Nonlinear scale of mode (m, n): the energy at which its stretching energy
235+// equals its bending energy, for the mode shape alone (bending ~ q^2,
236+// stretching ~ q^4): E = Vb(1)^2 / Vs(1) at unit amplitude.
237+export function nonlinearScale(plate: Plate, m: number, n: number): number {
238+ const v = modeShape(plate.g, m, n), F = new Float64Array(plate.g.n);
239+ const Vb = plate.forces(v, F, true), Vs = plate.stretching(v);
240+ return (Vb * Vb) / Vs;
241+}
242+
243+// Integrator step: STEP / omega at the top of the band, as in ../plate. The
244+// in-plane force per length N (from the stretching) raises the top mode's
245+// omega^2 by about N (4 / ax^2 + 4 / ay^2) / rho h; N is at most about
246+// sqrt(2 Eh E / (ax ay)) (all of E as stretching in one cell), taken twice
247+// for margin.
248+const STEP = 0.2;
249+export function stepFor(g: Grid, modes: Modes, E: number): number {
250+ const N = 2 * Math.sqrt((2 * EH * Math.max(E, 0)) / (g.ax * g.ay));
251+ return STEP / Math.sqrt(modes.omegaMax ** 2 + (N * (4 / g.ax ** 2 + 4 / g.ay ** 2)) / RHO_H);
252+}
src/driven-plate/style.cssadded+14−0View file
@@ -0,0 +1,14 @@
1+/* Styles for this app only, after ../common/style.css. */
2+
3+.row { display: flex; align-items: center; gap: 8px; flex-wrap: wrap; }
4+.row + .row { margin-top: 10px; }
5+canvas { display: block; width: 100%; }
6+#plate { aspect-ratio: 1.45; }
7+#chart { height: 240px; }
8+#readout { margin-top: 8px; font-size: 12px; color: var(--muted); white-space: pre-wrap; line-height: 1.6; }
9+.column { display: flex; flex-direction: column; gap: 16px; min-width: 0; }
10+.card-head { font-size: 12px; margin-bottom: 12px; display: flex; align-items: center; }
11+input[type="number"] { width: 4.5em; font-family: var(--mono); }
12+.label-w { min-width: 7.5em; }
13+.value-w { min-width: 6em; text-align: right; }
14+#params { margin-top: 12px; font: 12px/1.6 var(--mono); color: var(--muted); white-space: pre; overflow-x: auto; }
src/index.htmlmodified+1−0View file
@@ -35,6 +35,7 @@
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>
3737 <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, clamped or free; modes (m, n) at m² + n², pitch glides when struck hard.</div></a></li>
38+ <li><a href="./driven-plate/"><strong>Driven plate</strong><div class="muted">A rectangular plate driven in one mode, with no damping, its energy raised slowly: watch the energy move from the driven mode into the others.</div></a></li>
3839 </ul>
3940 </main>
4041 </body>
vite.config.tsmodified+1−0View file
@@ -20,6 +20,7 @@ export default defineConfig({
2020 "fput-chain": page("fput-chain/index.html"),
2121 beam: page("beam/index.html"),
2222 plate: page("plate/index.html"),
23+ "driven-plate": page("driven-plate/index.html"),
2324 },
2425 // Name the chunk holding src/common/ "common" (otherwise it takes the
2526 // name of one of its modules).