Plate: optional curvature (domed rest shape, Marguerre shallow shell): stretching measured from the rest shape, computed modes with the curvature stiffness, generalized nonlinear scale
4 changed files+161−55
src/plate/README.mdmodified+3−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; step 3 (2026-10-09), free edges; step 4 (2026-10-09), stretching with free edges (and the free layout moved to the grid vertices); step 5 (2026-10-09), clamped edges; step 6 (2026-10-09), Lyapunov exponent.
3+Struck square plate, step 4 of "Listening for chaos" (after double-pendulum, fput-chain and beam), the first two-dimensional system on the way to cymbal and gong. Work in progress: step 1 (2026-10-09), a linear, simply supported square plate with mallet strike, animation and sped-up note with spectrogram; step 2 (2026-10-09), von Kármán stretching; step 3 (2026-10-09), free edges; step 4 (2026-10-09), stretching with free edges (and the free layout moved to the grid vertices); step 5 (2026-10-09), clamped edges; step 6 (2026-10-09), Lyapunov exponent; step 7 (2026-10-09), curvature.
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: ghost, mode bars and histories, cross-section.
66
@@ -27,3 +27,5 @@ A membrane tension T for held edges (drum-like: springs between neighbours and t
2727 Cost issue, fixed (2026-10-09, `../common/note-render.ts`): the 4x oversampled note used to force at least 4 steps per audio sample, more than stability needs whenever the highest mode is low in sound (below about 6 kHz at 48 kHz: clamped above all). The renderer now steps at the stability step and reads the signal between steps by cubic interpolation through the four nearest steps. Checks: a linear undamped note (simply supported, 1 J, 10 s) has nothing above −140 dB away from the modes either way (old −144 dB, new −141 dB), with 1.9x fewer steps; damped nonlinear notes (2.51 J, h 5 cm) keep level and spectral centroid over time within the spread a 10% change of step gives anyway (chaotic notes: a different step gives a different, equivalent realization; for clamped, up to 2.7 dB and 100 Hz from a 10% step change, more than old against new). Browser, same conditions, default note: simply supported 7.2 → 5.1 s, clamped 11.9 → 5.2 s, free 9.3 → 6.1 s. Beam and FPUT chain unchanged in practice (their top modes sit above 6 kHz at the default pitch, so their steps were already shorter than the reading interval).
2828
2929 Lyapunov exponent (step 6, 2026-10-09): as the beam's (two-trajectory method, `../common/lyapunov.ts`), with the plate's shadow in `lyapunov.ts`: a second plate D0 = 1e-8 away, distances with displacements in units of A/√(N²) (A the lowest mode's amplitude at which stretching equals bending, `nonlinearAmplitude`; 1 m for a linear plate) and momenta in m ω₁ times that; the initial offset is a fixed pseudo-random pattern (free: its rigid part removed). Layout now as the beam's: Parameters folded under the plate, a Real time card on the right with Lyapunov (λ 10 s and since the strike, first 5 s left out of the y range) and Energy error charts. Sound panel: Lyapunov chart and λ vs E scatter, both off at first (the shadow doubles the time to compute a note: default note 4.3 s without, about 12 s with, in the browser). Exponents (2026-10-09, N = 8, kick at (1.3, 1.1), undamped, 300 s, λ since the strike, 1/s): linear plate 0.007 at every energy (finite-time; the last 100 s, −0.001: zero); simply supported, h 5 cm: 0.008, 0.065, 0.26, 0.53, 1.01 at 0.001, 0.1, 1, 10, 100 J; free, h 5 cm: 0.006, 0.015, 0.16, 1.54, 8.4. The free plate's exponent grows much faster with energy.
30+
31+Curvature (step 7, 2026-10-09; Parameters, rise, default 0): the plate rests, free of stress, in a shallow dome u0 of height H (held edges: H sin(πx/L) sin(πy/L), zero at the supports; free: the paraboloid H(1 − r²/c²) about the centre, a spherical cap), as the beam's curved rest shape. Marguerre shallow shell: bending unchanged (linear, from the rest shape), stretching from the rest shape, s = l(u0 + u, u0 + u) − l(u0, u0) = l(u, u) + 2 l(u0, u), the force its exact gradient (the second differences of the total shape in the scatter). The 2 l(u0, u) term adds linear stiffness, which lifts the low modes most (the cymbal effect), and a quadratic coupling. Acts only through the stretching (h > 0). With a rise the modes are computed (symmetry blocks; the rest shapes are mirror-symmetric), also simply supported, from the linearized forces (`forces(..., linear)` now includes the curvature stiffness; `linearEnergy` = bending + that quadratic part). Nonlinear scale and amplitude generalized: the lowest mode's amplitude at which the nonlinear stretching (cubic and quartic) equals the linear energy, by bisection (flat: unchanged, 0.131 J). Checks: force against central differences 2e-10 to 5e-10 (simply supported, clamped, free, H 0.1 m, N 6); linearized force against the force at 1e-7 amplitude 6e-8; free plate still invariant under rigid motion (2e-15); flat results unchanged. Continuum check (simply supported, H/h = 0.5, 1, 2, 4): the (1, 1) shape's Rayleigh quotient ω²/ω₀² is 1 + 0.836 (H/h)² in the continuum (single mode, exact stress solution, the same 0.209 coefficient as the flat check above, times 4); grid 1.2102, 1.8410, 4.3639, 14.456 at N = 32 against 1.2090, 1.8360, 4.3440, 14.376 (second order). Lowest modes, simply supported, N = 16, relative to the flat (1, 1): H = 0.05 m (= h): 1.35, 2.64, 2.64, 4.04; H = 0.1: 1.94, 3.03, 3.03, 4.24; H = 0.2: 2.90, 4.05, 4.05, 4.89 (flat: 1, 2.48, 2.48, 3.97). Energy error with a rise of 0.1 m at 100 J: 1e-3 (simply supported), 6e-4 (free).
src/plate/index.htmlmodified+3−2View file
@@ -56,6 +56,7 @@
5656 <label>m <input id="p-m" type="number" step="any" /> <span class="unit">kg</span></label>
5757 <label>κ <input id="p-kappa" type="number" step="any" /> <span class="unit">N/m (bending)</span></label>
5858 <label>h <input id="p-thick" type="number" step="any" /> <span class="unit">m (thickness, sets the stretching; 0 = linear)</span></label>
59+ <label>rise <input id="p-rise" type="number" step="any" /> <span class="unit">m (height of a domed rest shape; 0 = flat)</span></label>
5960 <label>a <input id="p-a" type="number" step="any" /> <span class="unit">m (spacing)</span></label>
6061 <label>τ <input id="p-tau" type="number" step="any" /> <span class="unit">s (decay, 0 = off)</span></label>
6162 <label>τ top <input id="p-tautop" type="number" step="any" /> <span class="unit">s (decay of the highest modes from the loss, 0 = off)</span></label>
@@ -117,10 +118,10 @@
117118 </template>
118119 <template id="help-params">
119120 <h2>Parameters</h2>
120- <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>
121+ <p>Mass, bending κ, the thickness h (it sets the stretching, as the amplitude at which it matters; 0 for a linear plate), the rise of a domed rest shape (0 for flat; it acts through the stretching, so it needs h > 0), 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>
121122 <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>
122123 <p>The stretching follows the von Kármán plate equations: the in-plane stress comes from a stress function, solved on the grid at every step. A simply supported edge may move outward but not along itself; a free edge is free in its own plane too. The nonlinear scale below is the energy of the lowest mode at which its stretching energy equals its bending energy.</p>
123- <p>With free edges the bending energy also includes the plate's twisting term (Poisson's ratio 0.3), which sets the edge conditions; the modes are computed, and their frequencies converge to plate theory (Leissa's tables) as the spacing squared (the lowest is 2.6% low at 9 × 9, 0.2% at 33 × 33). On a free plate the outer masses sit on the edge, carrying half a cell's mass (a quarter at the corners), so the free plate's side is one spacing less than the number of masses. Clamped edges hold the slope at zero too; the modes are computed (no (m, n)), and the plate is stiffer than a simply supported one, with overtones about as widely spread (plate theory: 1, 2.04, 3.01, 3.66, ...). 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>
124+ <p>With free edges the bending energy also includes the plate's twisting term (Poisson's ratio 0.3), which sets the edge conditions; the modes are computed, and their frequencies converge to plate theory (Leissa's tables) as the spacing squared (the lowest is 2.6% low at 9 × 9, 0.2% at 33 × 33). On a free plate the outer masses sit on the edge, carrying half a cell's mass (a quarter at the corners), so the free plate's side is one spacing less than the number of masses. Clamped edges hold the slope at zero too; the modes are computed (no (m, n)), and the plate is stiffer than a simply supported one, with overtones about as widely spread (plate theory: 1, 2.04, 3.01, 3.66, ...). A rise curves the plate at rest into a shallow dome (a shell, as a cymbal is): the plate rests in it without stress, and is drawn as its displacement from it. Bending is unchanged, but bending the dome stretches it in its own plane, which stiffens the plate: the lowest modes rise most (by about √(1 + 0.84 (rise / h)²) for the simply supported plate's lowest), so the low end of the spectrum crowds together, and the stretching gains a quadratic coupling besides the cubic one. With a rise the modes are computed. For a flat 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>
124125 </template>
125126 <template id="help-realtime">
126127 <h2>Real time</h2>
src/plate/main.tsmodified+7−6View file
@@ -42,7 +42,7 @@ const pauseButton = document.getElementById("pause") as HTMLButtonElement;
4242 // half an octave, 100 J is a broadband wash). base: the page's
4343 // parameters, for the plate at REF_N masses per side; P: the same plate
4444 // refined to the chosen N (see refine).
45-let base: Params = { N: REF_N, edges: "simply", m: 1, kappa: 20, thick: 0.05, a: 0.5, tau: 540, tauTop: 200 };
45+let base: Params = { N: REF_N, edges: "simply", m: 1, kappa: 20, thick: 0.05, rise: 0, a: 0.5, tau: 540, tauTop: 200 };
4646 let P: Params = base;
4747 // The mallet: a physical object, not refined with N.
4848 let mallet: MalletParams = { M: 1, K: 1000 };
@@ -185,6 +185,7 @@ const paramInputs = {
185185 m: document.getElementById("p-m") as HTMLInputElement,
186186 kappa: document.getElementById("p-kappa") as HTMLInputElement,
187187 thick: document.getElementById("p-thick") as HTMLInputElement,
188+ rise: document.getElementById("p-rise") as HTMLInputElement,
188189 a: document.getElementById("p-a") as HTMLInputElement,
189190 tau: document.getElementById("p-tau") as HTMLInputElement,
190191 tauTop: document.getElementById("p-tautop") as HTMLInputElement,
@@ -216,7 +217,7 @@ function distinctModes(): { w: number; count: number; k: number }[] {
216217 // and the decay times.
217218 function showDerived() {
218219 const w1 = modes.omega[0], hz = (w: number) => w / (2 * Math.PI);
219- const free = P.edges !== "simply"; // no (m, n)
220+ const free = modes.S !== null; // computed modes: no (m, n)
220221 const rows = distinctModes().slice(0, 8).map(({ w, count, k }) => {
221222 const m = modes.m[k], n = modes.n[k];
222223 return (
@@ -235,8 +236,8 @@ function showDerived() {
235236 (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`) +
236237 `side ${fmtLength(plateSide(P))}, ${modes.count} modes, ${hz(w1).toPrecision(3)} … ${hz(modes.omegaMax).toPrecision(3)} Hz ` +
237238 `(top ${(modes.omegaMax / w1).toFixed(1)}× the lowest)\n` +
238- `lowest modes${P.edges === "free" ? " (free edges, rigid motions removed)" : P.edges === "clamped" ? " (clamped, computed)" : " (m, n)"}:\n${rows.join("\n")}\n` +
239- `nonlinear scale ${Number.isFinite(nonlinearScale(P, modes)) ? fmtEnergy(nonlinearScale(P, modes)) : "none (linear)"} (${free ? "lowest mode" : "mode 1, 1"}: stretching = bending energy)\n` +
239+ `lowest modes${P.edges === "free" ? " (free edges, rigid motions removed)" : free ? " (computed)" : " (m, n)"}:\n${rows.join("\n")}\n` +
240+ `nonlinear scale ${Number.isFinite(nonlinearScale(P, modes)) ? fmtEnergy(nonlinearScale(P, modes)) : "none (linear)"} (${free ? "lowest mode" : "mode 1, 1"}: nonlinear stretching = linear energy)\n` +
240241 `decay time ${decay}\n` +
241242 `integrator step ${(1000 * H).toPrecision(3)} ms (from the strike energy; ` +
242243 `${(1000 * stepFor(P, modes, 0)).toPrecision(3)} ms when weak, ` +
@@ -441,10 +442,10 @@ edgesSelect.addEventListener("change", () => setParams({ ...base, edges: edgesSe
441442 document.getElementById("param-apply")!.addEventListener("click", () => {
442443 const v = Object.fromEntries(Object.entries(paramInputs).map(([key, el]) => [key, Number(el.value)]));
443444 const mv = { M: Number(malletInputs.M.value), K: Number(malletInputs.K.value) };
444- const ok = Object.entries(v).every(([key, x]) => Number.isFinite(x) && (x > 0 || (["tau", "tauTop", "thick"].includes(key) && x === 0))) &&
445+ const ok = Object.entries(v).every(([key, x]) => Number.isFinite(x) && (x > 0 || (["tau", "tauTop", "thick", "rise"].includes(key) && x === 0))) &&
445446 Object.values(mv).every((x) => Number.isFinite(x) && x > 0);
446447 if (!ok) {
447- paramStatus.textContent = "parameters must be positive (h, τ, τ top may be 0)";
448+ paramStatus.textContent = "parameters must be positive (h, rise, τ, τ top may be 0)";
448449 return;
449450 }
450451 paramStatus.textContent = "";
src/plate/plate.tsmodified+148−46View file
@@ -38,6 +38,19 @@
3838 // would need another solver. The stretching is set by the thickness h:
3939 // ks / kappa = 12 (1 - nu^2) / h^2 (from D = Eh^3 / 12 (1 - nu^2), nu = 0.3),
4040 // so it matters once the plate moves by about h.
41+// Curvature (rise H > 0): the plate rests, free of stress, in a dome u0 (a
42+// shallow shell, Marguerre's equations): held edges H sin(pi x / L) sin(pi
43+// y / L), zero at the supports; free, the paraboloid H (1 - r^2 / c^2)
44+// about the centre, c the half side (a spherical cap). Displacements u are
45+// measured from it. Bending is linear, so it is unchanged; the stretching is
46+// measured from the rest shape: s = l(u0 + u, u0 + u) - l(u0, u0) = l(u, u) +
47+// 2 l(u0, u), as the beam's e = d^2 + 2 D0 d. The term 2 l(u0, u) is linear
48+// in u: it adds linear stiffness, raising the low modes (in a cymbal, the
49+// dome lifts the lowest modes to a few hundred Hz), and a quadratic
50+// coupling besides the cubic one. Curvature acts only through the
51+// stretching: with h = 0 (no stretching) the rise has no effect. With a rise
52+// the modes are computed (L^2 is no longer diagonal in the coupled
53+// stiffness), also when simply supported.
4154 // Clamped: the supports also hold the slope at zero, an even ghost beyond
4255 // them (u(-2, j) = u(0, j)). The support nodes then have a bending term,
4356 // (L u) = -2 u(0, j) there, weighted 1/2 (the support's half cell): the
@@ -98,6 +111,7 @@ export interface Params {
98111 m: number; // kg, each mass
99112 kappa: number; // N/m, bending
100113 thick: number; // m, thickness, sets the stretching (von Karman); 0 = linear
114+ rise: number; // m, height of the curved rest shape (a dome); 0 = flat
101115 a: number; // m, rest spacing (for positions and drawing; it does not enter the equations)
102116 tau: number; // s, amplitude decay time from drag (all modes); 0 = no damping
103117 tauTop: number; // s, amplitude decay time of the top of the band from the loss alone; 0 = off
@@ -413,19 +427,24 @@ function inverseClamped(M: number, s: Float64Array, phi: Float64Array): void {
413427 for (let k = 0; k < M * M; k++) phi[k] = x[k] - w[k];
414428 }
415429
416-// Stretching energy of a free plate at u; with F given, its force is added.
417-function freeStretching(P: Params, ks: number, u: Float64Array, F?: Float64Array): number {
418- const N = P.N, M = N - 2, z = stretchScratch(M);
419- const { xx, yy, xy, s, phi } = z;
430+// Second differences of a free plate at its interior nodes (index q), from u.
431+function diffsFree(N: number, u: Float64Array, xx: Float64Array, yy: Float64Array, xy: Float64Array): void {
420432 for (let j = 1, q = 0; j < N - 1; j++) {
421433 for (let i = 1, k = 1 + N * j; i < N - 1; i++, k++, q++) {
422434 const cc = 2 * u[k];
423435 xx[q] = u[k - 1] - cc + u[k + 1];
424436 yy[q] = u[k - N] - cc + u[k + N];
425437 xy[q] = 0.25 * (u[k + N + 1] - u[k - N + 1] - u[k + N - 1] + u[k - N - 1]);
426- s[q] = 2 * (xx[q] * yy[q] - xy[q] * xy[q]);
427438 }
428439 }
440+}
441+
442+// Stretching energy of a free plate at u; with F given, its force is added.
443+function freeStretching(P: Params, ks: number, u: Float64Array, F: Float64Array | undefined, linear: boolean): number {
444+ const N = P.N, M = N - 2, z = stretchScratch(M);
445+ const { xx, yy, xy, s, phi } = z;
446+ diffsFree(N, u, xx, yy, xy);
447+ if (!sourceTerm(P, M * M, xx, yy, xy, s, linear)) return 0;
429448 inverseClamped(M, s, phi);
430449 let V = 0;
431450 for (let q = 0; q < M * M; q++) V += s[q] * phi[q];
@@ -451,6 +470,75 @@ function freeStretching(P: Params, ks: number, u: Float64Array, F?: Float64Array
451470 return (ks / 8) * V;
452471 }
453472
473+// The rest shape u0 (m) of each node; null when flat. Cached per Params.
474+const restCache = new WeakMap<Params, Float64Array | null>();
475+export function restShape(P: Params): Float64Array | null {
476+ if (restCache.has(P)) return restCache.get(P)!;
477+ let u0: Float64Array | null = null;
478+ if (P.rise) {
479+ const N = P.N, L = plateSide(P), c = L / 2;
480+ u0 = new Float64Array(N * N);
481+ for (let j = 0; j < N; j++) {
482+ for (let i = 0; i < N; i++) {
483+ const x = nodePosition(P, i), y = nodePosition(P, j);
484+ u0[i + N * j] = isFree(P)
485+ ? P.rise * (1 - ((x - c) ** 2 + (y - c) ** 2) / (c * c))
486+ : P.rise * Math.sin((Math.PI * x) / L) * Math.sin((Math.PI * y) / L);
487+ }
488+ }
489+ }
490+ restCache.set(P, u0);
491+ return u0;
492+}
493+// Second differences of the rest shape at the stretching's nodes, and its
494+// own s0 = l(u0, u0); null when flat.
495+interface RestDiffs { xx: Float64Array; yy: Float64Array; xy: Float64Array; s0: Float64Array }
496+const restDiffCache = new WeakMap<Params, RestDiffs | null>();
497+function restDiffs(P: Params): RestDiffs | null {
498+ if (restDiffCache.has(P)) return restDiffCache.get(P)!;
499+ const u0 = restShape(P);
500+ let r: RestDiffs | null = null;
501+ if (u0) {
502+ const N = P.N, n = isFree(P) ? (N - 2) * (N - 2) : N * N, f = () => new Float64Array(n);
503+ r = { xx: f(), yy: f(), xy: f(), s0: f() };
504+ if (isFree(P)) diffsFree(N, u0, r.xx, r.yy, r.xy);
505+ else diffsHeld(N, u0, new Float64Array((N + 2) * (N + 2)), r.xx, r.yy, r.xy);
506+ for (let k = 0; k < n; k++) r.s0[k] = 2 * (r.xx[k] * r.yy[k] - r.xy[k] * r.xy[k]);
507+ }
508+ restDiffCache.set(P, r);
509+ return r;
510+}
511+// The stretching's source s from the second differences of u (n values, in
512+// place): l(u, u) when flat; with a rest shape, l(u0 + u, u0 + u) - l(u0,
513+// u0), leaving the total shape's differences in xx, yy, xy for the force;
514+// linear: only the part linear in u, 2 l(u0, u), leaving the rest shape's
515+// differences (the force of the linearized energy). Returns false when s is
516+// zero (linear and flat).
517+function sourceTerm(P: Params, n: number, xx: Float64Array, yy: Float64Array, xy: Float64Array, s: Float64Array, linear: boolean): boolean {
518+ const r = restDiffs(P);
519+ if (!r) {
520+ if (linear) return false;
521+ for (let k = 0; k < n; k++) s[k] = 2 * (xx[k] * yy[k] - xy[k] * xy[k]);
522+ return true;
523+ }
524+ if (linear) {
525+ for (let k = 0; k < n; k++) {
526+ s[k] = 2 * (r.xx[k] * yy[k] + r.yy[k] * xx[k] - 2 * r.xy[k] * xy[k]);
527+ xx[k] = r.xx[k];
528+ yy[k] = r.yy[k];
529+ xy[k] = r.xy[k];
530+ }
531+ return true;
532+ }
533+ for (let k = 0; k < n; k++) {
534+ xx[k] += r.xx[k];
535+ yy[k] += r.yy[k];
536+ xy[k] += r.xy[k];
537+ s[k] = 2 * (xx[k] * yy[k] - xy[k] * xy[k]) - r.s0[k];
538+ }
539+ return true;
540+}
541+
454542 // Scratch for the stretching terms, by N. Padded arrays have a ring of
455543 // zeros beyond the edge, (N + 2)^2 values, node (i, j) at (i + 1) + (N + 2)(j + 1).
456544 interface Stretch { pad: Float64Array; xx: Float64Array; yy: Float64Array; xy: Float64Array; s: Float64Array; phi: Float64Array; a: Float64Array; b: Float64Array; c: Float64Array }
@@ -464,25 +552,32 @@ function stretchScratch(N: number): Stretch {
464552 }
465553 return z;
466554 }
467-// Stretching energy (J) at u; with F given, its force is added to F.
468-function stretching(P: Params, u: Float64Array, F?: Float64Array): number {
469- const ks = stretchCoefficient(P);
470- if (ks === 0) return 0;
471- if (isFree(P)) return freeStretching(P, ks, u, F);
472- const N = P.N, n = N * N, R = N + 2, z = stretchScratch(N);
473- const { pad, xx, yy, xy, s, phi, a, b, c } = z;
555+// Second differences of a held plate at every node (zero beyond the edge):
556+// xx, yy and the centred cross difference xy (pad: scratch, (N + 2)^2).
557+function diffsHeld(N: number, u: Float64Array, pad: Float64Array, xx: Float64Array, yy: Float64Array, xy: Float64Array): void {
558+ const R = N + 2;
474559 for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) pad[i + 1 + R * (j + 1)] = u[i + N * j];
475- // Second differences at every node (zero beyond the edge): xx, yy and the
476- // centred cross difference xy; s = l(u, u).
477560 for (let j = 0, k = 0; j < N; j++) {
478561 for (let i = 0, K = 1 + R * (j + 1); i < N; i++, k++, K++) {
479562 const cc = 2 * pad[K];
480563 xx[k] = pad[K - 1] - cc + pad[K + 1];
481564 yy[k] = pad[K - R] - cc + pad[K + R];
482565 xy[k] = 0.25 * (pad[K + R + 1] - pad[K - R + 1] - pad[K + R - 1] + pad[K - R - 1]);
483- s[k] = 2 * (xx[k] * yy[k] - xy[k] * xy[k]);
484566 }
485567 }
568+}
569+
570+// Stretching energy (J) at u; with F given, its force is added to F.
571+// linear: only the part quadratic in u (from 2 l(u0, u); zero when flat),
572+// and its (linear) force.
573+function stretching(P: Params, u: Float64Array, F?: Float64Array, linear = false): number {
574+ const ks = stretchCoefficient(P);
575+ if (ks === 0) return 0;
576+ if (isFree(P)) return freeStretching(P, ks, u, F, linear);
577+ const N = P.N, n = N * N, R = N + 2, z = stretchScratch(N);
578+ const { pad, xx, yy, xy, s, phi, a, b, c } = z;
579+ diffsHeld(N, u, pad, xx, yy, xy);
580+ if (!sourceTerm(P, n, xx, yy, xy, s, linear)) return 0;
486581 inverseL2(N, s, phi);
487582 let V = 0;
488583 for (let k = 0; k < n; k++) V += s[k] * phi[k];
@@ -511,11 +606,12 @@ function stretching(P: Params, u: Float64Array, F?: Float64Array): number {
511606 }
512607
513608 // Forces: bending, F = -kappa L (L u) (free: freeBending), plus stretching;
514-// c (if given) is scratch. linear: bending only.
609+// c (if given) is scratch. linear: the linearized forces (bending, and the
610+// curvature's linear stretching stiffness).
515611 export function forces(P: Params, u: Float64Array, F: Float64Array, c: Float64Array = new Float64Array(u.length), linear = false): void {
516612 if (isFree(P)) {
517613 freeBending(P, u, F);
518- if (!linear) stretching(P, u, F);
614+ stretching(P, u, F, linear);
519615 return;
520616 }
521617 laplacian(P.N, u, c);
@@ -525,7 +621,7 @@ export function forces(P: Params, u: Float64Array, F: Float64Array, c: Float64Ar
525621 const N = P.N;
526622 for (let j = 0, k = 0; j < N; j++) for (let i = 0; i < N; i++, k++) F[k] -= 2 * P.kappa * touched(N, i, j) * u[k];
527623 }
528- if (!linear) stretching(P, u, F);
624+ stretching(P, u, F, linear);
529625 }
530626
531627 export function kinetic(P: Params, s: State): number {
@@ -534,8 +630,8 @@ export function kinetic(P: Params, s: State): number {
534630 for (let k = 0; k < s.p.length; k++) T += (s.p[k] * s.p[k]) / M[k];
535631 return T / 2;
536632 }
537-// Linear energy at u: bending.
538-export function linearEnergy(P: Params, u: Float64Array): number {
633+// Bending energy at u.
634+export function bendingEnergy(P: Params, u: Float64Array): number {
539635 if (isFree(P)) return freeBending(P, u);
540636 const N = P.N, c = new Float64Array(u.length);
541637 laplacian(N, u, c);
@@ -544,8 +640,11 @@ export function linearEnergy(P: Params, u: Float64Array): number {
544640 if (isClamped(P)) for (let j = 0, k = 0; j < N; j++) for (let i = 0; i < N; i++, k++) V += 2 * touched(N, i, j) * u[k] * u[k];
545641 return (P.kappa * V) / 2;
546642 }
643+// Linear energy at u (the quadratic part of the potential): bending, plus
644+// the curvature's linear stretching stiffness.
645+export const linearEnergy = (P: Params, u: Float64Array) => bendingEnergy(P, u) + stretching(P, u, undefined, true);
547646 export const stretchEnergy = (P: Params, u: Float64Array) => stretching(P, u);
548-export const potential = (P: Params, s: State) => linearEnergy(P, s.u) + stretching(P, s.u);
647+export const potential = (P: Params, s: State) => bendingEnergy(P, s.u) + stretching(P, s.u);
549648 export const energy = (P: Params, s: State) => kinetic(P, s) + potential(P, s);
550649
551650 // Velocity Verlet, with drag and loss in half steps around it (Strang), as
@@ -642,7 +741,7 @@ export class Modes {
642741 constructor(P: Params) {
643742 const N = P.N;
644743 this.N = N;
645- if (isFree(P) || isClamped(P)) {
744+ if (isFree(P) || isClamped(P) || restDiffs(P) && stretchCoefficient(P) > 0) {
646745 const { lambda, S } = computedModes(P);
647746 this.count = lambda.length;
648747 this.omega = Float64Array.from(lambda, (l) => Math.sqrt(Math.max(l, 0)));
@@ -739,33 +838,36 @@ export function decayTimes(P: Params, modes: Modes): Float64Array {
739838 });
740839 }
741840
742-// Nonlinear scale (J): the energy of the lowest mode ((1, 1) when simply
743-// supported) at which its stretching energy equals its linear energy
744-// (Infinity when linear). Bending goes as the amplitude squared, stretching
745-// as its fourth power.
746-export function nonlinearScale(P: Params, modes: Modes): number {
747- const { V2, V4 } = lowestModeEnergies(P, modes);
748- return V4 > 0 ? (V2 * V2) / V4 : Infinity;
749-}
750-// Bending (V2) and stretching (V4) energy of the lowest mode's shape scaled
751-// to unit length (sum of u^2 = 1 m^2); 0 stretching when linear.
752-function lowestModeEnergies(P: Params, modes: Modes): { V2: number; V4: number } {
841+// Nonlinear amplitude (m) and scale (J): the lowest mode's shape u1 (unit
842+// length, sum of u^2 = 1 m^2) at amplitude q has linear energy V2(q) (bending
843+// and the curvature's stiffness, ~ q^2) and nonlinear stretching energy
844+// Vn(q) (the rest, ~ q^3 and q^4); the amplitude is the q at which |Vn| =
845+// V2, the scale V2 there. Flat: q^2 = V2(1) / Vn(1). Infinity when linear.
846+function nonlinearPoint(P: Params, modes: Modes): { q: number; E: number } {
847+ if (stretchCoefficient(P) === 0) return { q: Infinity, E: Infinity };
753848 const N = P.N;
754- let u = new Float64Array(N * N);
755- if (modes.S) u = Float64Array.from(modes.S[0]);
756- else for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) u[i + N * j] = Math.sin((Math.PI * (i + 1)) / (N + 1)) * Math.sin((Math.PI * (j + 1)) / (N + 1));
849+ let u1 = new Float64Array(N * N);
850+ if (modes.S) u1 = Float64Array.from(modes.S[0]);
851+ else for (let j = 0; j < N; j++) for (let i = 0; i < N; i++) u1[i + N * j] = Math.sin((Math.PI * (i + 1)) / (N + 1)) * Math.sin((Math.PI * (j + 1)) / (N + 1));
757852 let n2 = 0;
758- for (const x of u) n2 += x * x;
759- for (let k = 0; k < u.length; k++) u[k] /= Math.sqrt(n2);
760- return { V2: linearEnergy(P, u), V4: stretching(P, u) };
761-}
762-// Nonlinear amplitude (m): the length of the lowest mode's displacement
763-// vector at which its stretching energy equals its bending energy (Infinity
764-// when linear).
765-export function nonlinearAmplitude(P: Params, modes: Modes): number {
766- const { V2, V4 } = lowestModeEnergies(P, modes);
767- return V4 > 0 ? Math.sqrt(V2 / V4) : Infinity;
853+ for (const x of u1) n2 += x * x;
854+ for (let k = 0; k < u1.length; k++) u1[k] /= Math.sqrt(n2);
855+ const at = (q: number) => {
856+ const u = u1.map((x) => q * x), V2 = linearEnergy(P, u);
857+ return { V2, ratio: Math.abs(stretching(P, u) - stretching(P, u, undefined, true)) / V2 };
858+ };
859+ // Geometric bisection on q between 1e-9 and 1e3 m.
860+ let lo = 1e-9, hi = 1e3;
861+ if (at(hi).ratio < 1) return { q: Infinity, E: Infinity };
862+ for (let it = 0; it < 80; it++) {
863+ const mid = Math.sqrt(lo * hi);
864+ if (at(mid).ratio < 1) lo = mid;
865+ else hi = mid;
866+ }
867+ return { q: hi, E: at(hi).V2 };
768868 }
869+export const nonlinearScale = (P: Params, modes: Modes) => nonlinearPoint(P, modes).E;
870+export const nonlinearAmplitude = (P: Params, modes: Modes) => nonlinearPoint(P, modes).q;
769871
770872 // Integrator step (s) for motion with energy E: STEP / omega, omega at the
771873 // top of the band, as in the beam. Stretching adds an in-plane tension T