1function result = plate_implicit(params)
2% PLATE_IMPLICIT Implicit finite differences (theta scheme) for the damped Kirchhoff-Love plate.
3%
4% result = plate_implicit(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 the same grid and with the same spatial stencils as plate_fd, but the
9% stiffness term is averaged over three time levels,
10%
11% (w+ - 2w + w-)/dt^2 + sigma0 (w+ - w-)/dt - sigma1 L (w+ - w-)/dt
12% + kappa^2 B (theta w+ + (1 - 2 theta) w + theta w-) = F s / mu,
13%
14% with L and B the sparse Laplacian and biharmonic matrices on the interior
15% nodes (plate_matrices). Multiplying by dt^2 and collecting terms gives one
16% linear system per step with a constant matrix,
17%
18% A w+ = C w - A_old w- + dt^2 F s / mu
19% A = (1 + sigma0 dt) I - sigma1 dt L + theta kappa^2 dt^2 B
20% C = 2 I - (1 - 2 theta) kappa^2 dt^2 B
21% A_old = (1 - sigma0 dt) I + sigma1 dt L + theta kappa^2 dt^2 B.
22%
23% A is symmetric positive definite, so it is factored once by sparse Cholesky
24% and each step costs two sparse triangular solves plus two sparse products.
25% (R' is formed once; forming it every step would cost more than the solves.)
26% For theta >= 1/4 the scheme is unconditionally stable (von Neumann analysis
27% of the amplification polynomial); theta = 1/4 is the trapezoidal rule. The
28% step is therefore chosen for accuracy, dt = 1 / (fs steps_per_sample).
29%
30% Approximations. Second order in space and time. Implicit averaging
31% lengthens the period of a mode by about (omega dt)^2 / 12, so the highest
32% modes are detuned rather than unstable. Probe, frames, stress envelope and
33% audio follow the conventions of plate_fd (see plate_result_schema).
35t_start = tic;
36P = plate_setup(params);
37plan = plate_implicit_plan(params);
38nx = plan.nx; ny = plan.ny; dx = plan.dx; dy = plan.dy; dt = plan.dt; nsteps = plan.nsteps;
39theta = P.implicit.theta;
40x = linspace(0, P.plate.Lx, nx);
41y = linspace(0, P.plate.Ly, ny);
42[X, Y] = meshgrid(x, y);
43fs = P.output.fs;
44notes = {};
45if theta < 0.25
46 notes{end + 1} = sprintf('theta = %g is below 1/4, so the scheme is only conditionally stable', theta);
47end
49% ---- Matrices on the interior nodes, factored once
50[L, B, mx, my] = plate_matrices(nx, ny, dx, dy, P.bc);
51m = mx * my;
52I = speye(m);
53k2dt2 = P.kappa^2 * dt^2;
54A = (1 + P.damping.sigma0 * dt) * I - P.damping.sigma1 * dt * L + theta * k2dt2 * B;
55C = 2 * I - (1 - 2 * theta) * k2dt2 * B;
56A_old = (1 - P.damping.sigma0 * dt) * I + P.damping.sigma1 * dt * L + theta * k2dt2 * B;
57t_factor = tic;
58R = chol(A); % A = R' * R
59Rt = R'; % transposed once, not once per step
60factor_s = toc(t_factor);
62% ---- Load, probe, frames, listener (as in plate_fd)
63s = plate_force_footprint(P.excitation, X, Y, dx * dy);
64s_int = s(2:end - 1, 2:end - 1);
65load = (dt^2 / P.mu) * s_int(:);
66tt = (0:nsteps) * dt;
67F = plate_force_time(P.excitation, tt);
69[~, jp] = min(abs(x - P.probe.x)); jp = min(max(jp, 2), nx - 1);
70[~, ip] = min(abs(y - P.probe.y)); ip = min(max(ip, 2), ny - 1);
71if abs(x(jp) - P.probe.x) > 0.5 * dx || abs(y(ip) - P.probe.y) > 0.5 * dy
72 notes{end + 1} = sprintf('probe moved to the nearest interior node (%.4g, %.4g) m', x(jp), y(ip));
73end
75frame_steps = unique(max(1, min(nsteps, round(linspace(0, P.implicit.duration, P.output.nframes) / dt))));
76nf = numel(frame_steps);
77Wf = zeros(ny, nx, nf); VMf = zeros(ny, nx, nf); SPf = zeros(ny, nx, nf);
78tf = zeros(1, nf);
80Lis = plate_listener(P.listener, X, Y, dx * dy, fs);
81bins = Lis.delay(:) + 1;
82wgt = Lis.weight(:);
83na = ceil(nsteps * dt * fs) + Lis.nbins + 1;
84Pbins = zeros(Lis.nbins, na);
86% ---- Time stepping on the interior vector
87w = zeros(m, 1); wm = zeros(m, 1);
88wfull = zeros(ny, nx);
89w_p = zeros(1, nsteps + 1); wxx_p = zeros(1, nsteps + 1); wyy_p = zeros(1, nsteps + 1); wxy_p = zeros(1, nsteps + 1);
90blowup = 1e3 * P.plate.h;
91diverged_step = 0;
92kf = 1;
94for nstep = 1:nsteps
95 rhs = C * w - A_old * wm + F(nstep) * load;
96 wn = R \ (Rt \ rhs); % forward then back substitution
98 wfull(2:end - 1, 2:end - 1) = reshape(wn, my, mx);
99 w_p(nstep + 1) = wfull(ip, jp);
100 wxx_p(nstep + 1) = (wfull(ip, jp + 1) - 2 * wfull(ip, jp) + wfull(ip, jp - 1)) / dx^2;
101 wyy_p(nstep + 1) = (wfull(ip + 1, jp) - 2 * wfull(ip, jp) + wfull(ip - 1, jp)) / dy^2;
102 wxy_p(nstep + 1) = (wfull(ip + 1, jp + 1) - wfull(ip + 1, jp - 1) - wfull(ip - 1, jp + 1) + wfull(ip - 1, jp - 1)) / (4 * dx * dy);
104 acc = zeros(ny, nx);
105 acc(2:end - 1, 2:end - 1) = reshape((wn - 2 * w + wm) / dt^2, my, mx);
106 ia = floor(nstep * dt * fs) + 1;
107 Pbins(:, ia) = Pbins(:, ia) + accumarray(bins, wgt .* acc(:), [Lis.nbins 1]) * (dt * fs);
109 if kf <= nf && nstep == frame_steps(kf)
110 [wxx, wyy, wxy] = plate_curvature(plate_pad(wfull, P.bc), dx, dy);
111 S = plate_stress(P, wxx, wyy, wxy);
112 Wf(:, :, kf) = wfull; VMf(:, :, kf) = S.vm; SPf(:, :, kf) = S.sp; tf(kf) = nstep * dt;
113 kf = kf + 1;
114 end
116 peak = max(abs(wn));
117 if ~(peak < blowup)
118 diverged_step = nstep;
119 notes{end + 1} = sprintf('solution diverged at step %d (t = %.3g s); results are truncated there', nstep, nstep * dt);
120 break
121 end
122 wm = w; w = wn;
123end
125if diverged_step > 0
126 nkeep = diverged_step;
127 tt = tt(1:nkeep); w_p = w_p(1:nkeep); wxx_p = wxx_p(1:nkeep); wyy_p = wyy_p(1:nkeep); wxy_p = wxy_p(1:nkeep);
128 keepf = frame_steps < diverged_step;
129 Wf = Wf(:, :, keepf); VMf = VMf(:, :, keepf); SPf = SPf(:, :, keepf); tf = tf(keepf);
130end
132p_audio = zeros(1, na);
133for d = 1:Lis.nbins
134 p_audio(d:end) = p_audio(d:end) + Pbins(d, 1:na - d + 1);
135end
136n_audio = min(na, ceil(numel(tt) * dt * fs));
137p_audio = p_audio(1:n_audio);
139Sp = plate_stress(P, wxx_p, wyy_p, wxy_p);
140[fw, Aw] = plate_spectrum(tt, w_p);
141[~, Avm] = plate_spectrum(tt, Sp.vm);
143result.method = 'implicit';
144result.grid = struct('x', x, 'y', y);
145result.t = tt;
146result.probe = struct('x', x(jp), 'y', y(ip), 'w', w_p, 'sxx', Sp.sxx, 'syy', Sp.syy, 'txy', Sp.txy, 'vm', Sp.vm, 'sp', Sp.sp);
147result.frames = struct('t', tf, 'w', Wf, 'vm', VMf, 'sp', SPf);
148if isempty(tf)
149 result.envelope = struct('vm', zeros(ny, nx), 'sp', zeros(ny, nx));
150else
151 result.envelope = struct('vm', max(VMf, [], 3), 'sp', max(abs(SPf), [], 3));
152end
153result.spectrum = struct('f', fw, 'w', Aw, 'vm', Avm);
154result.audio = struct('fs', fs, 'p', p_audio, 't0', Lis.t0);
155result.info = struct('nx', nx, 'ny', ny, 'dx', dx, 'dy', dy, 'dt', dt, 'theta', theta, ...
156 'steps_per_sample', P.implicit.steps_per_sample, 'nunknowns', m, 'nnz_factor', nnz(R), 'factor_s', factor_s, ...
157 'nsteps', nsteps, 'ncells', nx * ny, 'cell_steps', plan.cell_steps, 'fs', fs, ...
158 'stable', diverged_step == 0, 'diverged_step', diverged_step, 'wall_s', toc(t_start), 'notes', {notes});
159end