17959acInteractive tutorial on the theory and the solverJeremy Magland 1// Why the stiff terms are integrated implicitly. The test problem is one
2// undamped mode, u'' = -w^2 u, stepped by the paper's scheme
3// (1 + g z) u+ = (2 - (1 - 2g) z) u - (1 + g z) u-, z = (w ht)^2,
4// which is the explicit central difference at g = 0. The slider pushes the
5// mode frequency past the explicit ceiling w ht = 2 while the implicit
6// scheme stays on the unit circle.
8import { useState } from 'react';
9import { useCanvas, drawFrame, PLOT_COLORS } from '../../components/plots/canvasUtils';
10import { Figure, FigControls, FigSlider, usePlayer } from '../figparts';
12const HEIGHT = 250;
13const G = 0.2498;
14const NSTEPS = 90;
16function stepScheme(z: number, g: number): { u: number[]; rho: number } {
17 const A = 1 + g * z;
18 const B = 2 - (1 - 2 * g) * z;
19 // exact start: u(-dt) = u(0) cos(w dt) with w dt = sqrt(z)
20 const u: number[] = [1, Math.cos(Math.sqrt(z))];
21 // recurrence runs forward from (u[1], u[0])
22 let um = u[0];
23 let uc = u[1];
24 for (let k = 2; k < NSTEPS; k++) {
25 const up = (B * uc - A * um) / A;
26 u.push(up);
27 um = uc;
28 uc = up;
29 }
30 // spectral radius of the two-step recurrence
31 const disc = B * B - 4 * A * A;
32 let rho: number;
33 if (disc <= 0) rho = 1;
34 else rho = (Math.abs(B) + Math.sqrt(disc)) / (2 * A);
35 return { u, rho };
36}
38// For the audio comparison: the frequency the implicit scheme actually
39// produces for this mode, from the phase angle of its recurrence roots.
40function schemeFreq(z: number, g: number, stepRate: number): number {
41 const A = 1 + g * z;
42 const B = 2 - (1 - 2 * g) * z;
43 const theta = Math.atan2(Math.sqrt(Math.max(0, 4 * A * A - B * B)), B);
44 return (theta / (2 * Math.PI)) * stepRate;
45}
47const AUDIO_SR = 24000;
48// A pretend step rate low enough that the slider's whole range is audible.
49const DEMO_STEP_RATE = 8000;
51function renderDetuneAudio(fTrue: number, fNum: number): Float32Array<ArrayBuffer> {
52 const seg = 0.4;
53 const gap = 0.1;
54 const n = Math.round(AUDIO_SR * (2 * seg + gap));
55 const y = new Float32Array(n);
56 const tone = (f: number, start: number) => {
57 const n0 = Math.round(start * AUDIO_SR);
58 const n1 = n0 + Math.round(seg * AUDIO_SR);
59 const fade = Math.round(0.012 * AUDIO_SR);
60 for (let k = n0; k < n1; k++) {
61 let g2 = 0.35;
62 if (k - n0 < fade) g2 *= (k - n0) / fade;
63 if (n1 - k < fade) g2 *= (n1 - k) / fade;
64 y[k] = g2 * Math.sin((2 * Math.PI * f * (k - n0)) / AUDIO_SR);
65 }
66 };
67 tone(fTrue, 0);
68 tone(fNum, seg + gap);
69 return y;
70}
72export default function StabilityFig() {
73 const [wht, setWht] = useState(2.2);
74 const play = usePlayer();
76 const z = wht * wht;
77 const ex = stepScheme(z, 0);
78 const im = stepScheme(z, G);
79 const fTrue = (wht / (2 * Math.PI)) * DEMO_STEP_RATE;
80 const fNum = schemeFreq(z, G, DEMO_STEP_RATE);
82 const ref = useCanvas(
83 (ctx, w, h) => {
84 drawFrame(ctx, w, h);
85 const panels: { u: number[]; rho: number; name: string; color: string }[] = [
86 { ...ex, name: 'explicit (g = 0)', color: PLOT_COLORS.warn },
87 { ...im, name: `implicit (g = ${G})`, color: PLOT_COLORS.trace },
88 ];
89 const pw = (w - 40) / 2;
90 panels.forEach((p, pi) => {
91 const x0 = 14 + pi * (pw + 12);
92 const y0 = 12;
93 const ph = h - 50;
94 ctx.strokeStyle = PLOT_COLORS.grid;
95 ctx.strokeRect(x0, y0, pw, ph);
96 ctx.beginPath();
97 ctx.moveTo(x0, y0 + ph / 2);
98 ctx.lineTo(x0 + pw, y0 + ph / 2);
99 ctx.stroke();
101 // clip to the panel so divergence does not spray over the figure
102 ctx.save();
103 ctx.beginPath();
104 ctx.rect(x0, y0, pw, ph);
105 ctx.clip();
106 ctx.strokeStyle = p.color;
107 ctx.lineWidth = 1.6;
108 ctx.beginPath();
109 for (let k = 0; k < NSTEPS; k++) {
110 const x = x0 + (k / (NSTEPS - 1)) * pw;
111 const y = y0 + ph / 2 - Math.max(-1e6, Math.min(1e6, p.u[k])) * (ph / 2.4);
112 if (k === 0) ctx.moveTo(x, y);
113 else ctx.lineTo(x, y);
114 }
115 ctx.stroke();
116 ctx.lineWidth = 1;
117 ctx.restore();
119 ctx.fillStyle = p.color;
120 ctx.font = '11px system-ui';
121 ctx.textAlign = 'left';
122 const growth =
123 p.rho <= 1 + 1e-12 ? 'bounded (|amplification| = 1)' : `grows ×${p.rho.toFixed(3)} per step`;
124 ctx.fillText(`${p.name}: ${growth}`, x0 + 6, y0 + ph + 16);
125 });
127 ctx.fillStyle = PLOT_COLORS.axis;
128 ctx.textAlign = 'center';
129 ctx.fillText(
130 `one mode at ω·ht = ${wht.toFixed(2)} (explicit ceiling: 2)`,
131 w / 2,
132 h - 8,
133 );
134 },
135 [wht, z],
136 HEIGHT,
137 );
139 return (
140 <Figure
141 caption={
142 <>
143 Ninety steps of one undamped mode under both schemes, at the frequency set by the
144 slider. The explicit central difference is exact in cost and cheap in memory, but
145 the moment ω·h<sub>t</sub> passes 2 its recurrence has a root outside
146 the unit circle and the mode explodes; a plate's biharmonic stiffness puts its
147 highest grid modes far past that ceiling at any audio time step. With the implicit
148 weighting at g ≥ 1/4 the roots stay on the unit circle at every frequency (the
149 paper's g = 0.2498 sits just below 1/4, and is stable to ω·h
150 <sub>t</sub> ≈ 71, far beyond anything the grid carries). The price of the
151 implicit form is a linear solve per step, which Section 7 makes affordable, and a
152 detuning that grows with frequency: the play button compares the mode's true pitch
153 against the pitch the scheme actually produces, scaled into the audible range by
154 pretending the step rate were 8 kHz. Near the left end of the slider the two tones
155 are indistinguishable, which is the regime the audible modes live in; the real
156 simulation steps at 96 kHz, so a mode this badly flattened would sit far above the
157 spectrum that matters.
158 </>
159 }
160 >
161 <canvas ref={ref} />
162 <FigControls>
163 <FigSlider
164 label="ω·ht"
165 value={wht}
166 min={0.3}
167 max={4.5}
168 step={0.01}
169 format={(v) => v.toFixed(2)}
170 onChange={setWht}
171 />
172 <button onClick={() => play(renderDetuneAudio(fTrue, fNum), AUDIO_SR)}>
173 ▶ Hear the detuning: {fTrue.toFixed(0)} Hz true, then {fNum.toFixed(0)} Hz as the
174 scheme plays it
175 </button>
176 </FigControls>
177 </Figure>
178 );
179}