1function result = plate_fd(params)
2% PLATE_FD Explicit finite differences for the damped Kirchhoff-Love plate.
3%
4% result = plate_fd(params) advances
5%
6% w_tt + 2 sigma0 w_t - 2 sigma1 lap(w_t) + kappa^2 lap(lap(w)) = F(t) s(x,y) / mu
7%
8% on a uniform nx x ny grid (nodes include the edges) with the scheme
9%
10% (1 + sigma0 dt) w+ = 2 w - (1 - sigma0 dt) w- - dt^2 kappa^2 L(L w)
11% + 2 sigma1 dt (L w - L w-) + dt^2 F s / mu
12%
13% where w-, w, w+ are consecutive time levels and L is the five-point
14% Laplacian. Applying L twice gives the 13-point biharmonic stencil. The
15% scheme is second order in space and time (Bilbao, Numerical Sound
16% Synthesis, 2009, ch. 12). On a rectangle, simply supported and clamped
17% edges enter through one layer of ghost nodes (plate_pad). Free edges have
18% no ghost-layer form, so for them L and L(L .) are the sparse energy-method
19% operators of plate_matrices applied to every node, and the rigid-body
20% motions are projected out of the displacement after each step
21% (plate_rigid_modes). The same path serves the centre-pinned plate, whose
22% pinned node is zeroed after each step, and the ellipse (plate_domain),
23% whose edge is a staircase on the grid with every condition assembled as a
24% sparse operator.
25%
26% The step is dt = cfl * dt_max with dt_max from plate_fd_plan; the scheme
27% is unstable for cfl > 1 and the run stops as soon as the field diverges,
28% reporting the step at which it did (info.diverged_step).
29%
30% The work is in three stages, each its own file, so that a host can run the
31% time loop elsewhere (the lab has a WebGPU port of plate_fd_loop) while the
32% setup and the output stay this code:
33% S = plate_fd_setup(params) grid, step, load, probe, listener, operators
34% R = plate_fd_loop(S) the time stepping; records probe, frames, audio bins
35% result = plate_fd_output(S, R) stresses, audio, spectra, the result struct
36%
37% Approximations. Second-order truncation error in dx, dy, dt; the probe is
38% snapped to the nearest grid node; stress fields are evaluated only at the
39% saved frames, so the envelope is a sample of the true peak; the audio is
40% the Rayleigh integral box-averaged onto audio samples.
41%
42% See plate_result_schema for the fields of result.
44S = plate_fd_setup(params);
45R = plate_fd_loop(S);
46result = plate_fd_output(S, R);
47end