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 of plate_matrices,
15% which also supplies the diagonal mass matrix M: the identity on the interior
16% nodes for simply supported or clamped edges, the trapezoidal node weights on
17% every node for free edges. Multiplying by dt^2 and collecting terms gives
18% one linear system per step with a constant matrix,
19%
20% A w+ = C w - A_old w- + dt^2 F M s / mu
21% A = (1 + sigma0 dt) M - sigma1 dt L + theta kappa^2 dt^2 B
22% C = 2 M - (1 - 2 theta) kappa^2 dt^2 B
23% A_old = (1 - sigma0 dt) M + sigma1 dt L + theta kappa^2 dt^2 B.
24%
25% For free edges the three rigid-body motions are projected out of w+ after
26% every solve (plate_rigid_modes).
27%
28% A is symmetric positive definite, so it is factored once by sparse Cholesky
29% and each step costs two sparse triangular solves plus two sparse products.
30% (R' is formed once; forming it every step would cost more than the solves.)
31% For theta >= 1/4 the scheme is unconditionally stable (von Neumann analysis
32% of the amplification polynomial); theta = 1/4 is the trapezoidal rule. The
33% step is therefore chosen for accuracy, dt = 1 / (fs steps_per_sample).
34%
35% Approximations. Second order in space and time. Implicit averaging
36% lengthens the period of a mode by about (omega dt)^2 / 12, so the highest
37% modes are detuned rather than unstable. Probe, frames, stress envelope and
38% audio follow the conventions of plate_fd (see plate_result_schema).
40t_start = tic;
41P = plate_setup(params);
42plan = plate_implicit_plan(params);
43nx = plan.nx; ny = plan.ny; dx = plan.dx; dy = plan.dy; dt = plan.dt; nsteps = plan.nsteps;
44theta = P.implicit.theta;
45x = linspace(0, P.plate.Lx, nx);
46y = linspace(0, P.plate.Ly, ny);
47[X, Y] = meshgrid(x, y);
48fs = P.output.fs;
49notes = {};
50if theta < 0.25
51 notes{end + 1} = sprintf('theta = %g is below 1/4, so the scheme is only conditionally stable', theta);
52end
54% ---- Matrices on the unknowns, factored once
55free = strcmp(P.bc, 'free');
56[L, B, M, idx, mx, my] = plate_matrices(nx, ny, dx, dy, P.bc, P.material.nu);
57m = numel(idx);
58mvec = full(diag(M));
59if free
60 rr = 1:ny; cc = 1:nx; % every node is an unknown
61else
62 rr = 2:ny - 1; cc = 2:nx - 1; % the edge nodes carry w = 0
63end
64k2dt2 = P.kappa^2 * dt^2;
65A = (1 + P.damping.sigma0 * dt) * M - P.damping.sigma1 * dt * L + theta * k2dt2 * B;
66C = 2 * M - (1 - 2 * theta) * k2dt2 * B;
67A_old = (1 - P.damping.sigma0 * dt) * M + P.damping.sigma1 * dt * L + theta * k2dt2 * B;
68t_factor = tic;
69R = chol(A); % A = R' * R
70Rt = R'; % transposed once, not once per step
71factor_s = toc(t_factor);
73% ---- Load, probe, frames, listener (as in plate_fd)
74s = plate_force_footprint(P.excitation, X, Y, dx * dy);
75s_unk = s(rr, cc);
76load = (dt^2 / P.mu) * (mvec .* s_unk(:)); % M s on the unknowns
77if free
78 [Rb, Qb] = plate_rigid_modes(X(rr, cc), Y(rr, cc), mvec);
79 notes{end + 1} = 'free edges: the rigid-body translation and tilts are removed from the displacement after every step';
80end
81tt = (0:nsteps) * dt;
82F = plate_force_time(P.excitation, tt);
84[~, jp] = min(abs(x - P.probe.x)); jp = min(max(jp, 2), nx - 1);
85[~, ip] = min(abs(y - P.probe.y)); ip = min(max(ip, 2), ny - 1);
86if abs(x(jp) - P.probe.x) > 0.5 * dx || abs(y(ip) - P.probe.y) > 0.5 * dy
87 notes{end + 1} = sprintf('probe moved to the nearest interior node (%.4g, %.4g) m', x(jp), y(ip));
88end
90frame_steps = unique(max(1, min(nsteps, round(linspace(0, P.implicit.duration, P.output.nframes) / dt))));
91nf = numel(frame_steps);
92Wf = zeros(ny, nx, nf); VMf = zeros(ny, nx, nf); SPf = zeros(ny, nx, nf);
93tf = zeros(1, nf);
95Lis = plate_listener(P.listener, X, Y, dx * dy, fs);
96bins = Lis.delay(:) + 1;
97wgt = Lis.weight(:);
98na = ceil(nsteps * dt * fs) + Lis.nbins + 1;
99Pbins = zeros(Lis.nbins, na);
101% ---- Time stepping on the vector of unknowns
102w = zeros(m, 1); wm = zeros(m, 1);
103wfull = zeros(ny, nx);
104w_p = zeros(1, nsteps + 1); wxx_p = zeros(1, nsteps + 1); wyy_p = zeros(1, nsteps + 1); wxy_p = zeros(1, nsteps + 1);
105blowup = 1e3 * P.plate.h;
106diverged_step = 0;
107kf = 1;
109for nstep = 1:nsteps
110 rhs = C * w - A_old * wm + F(nstep) * load;
111 wn = R \ (Rt \ rhs); % forward then back substitution
112 if free
113 wn = wn - Rb * (Qb' * wn);
114 end
116 wfull(rr, cc) = reshape(wn, my, mx);
117 w_p(nstep + 1) = wfull(ip, jp);
118 wxx_p(nstep + 1) = (wfull(ip, jp + 1) - 2 * wfull(ip, jp) + wfull(ip, jp - 1)) / dx^2;
119 wyy_p(nstep + 1) = (wfull(ip + 1, jp) - 2 * wfull(ip, jp) + wfull(ip - 1, jp)) / dy^2;
120 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);
122 acc = zeros(ny, nx);
123 acc(rr, cc) = reshape((wn - 2 * w + wm) / dt^2, my, mx);
124 ia = floor(nstep * dt * fs) + 1;
125 Pbins(:, ia) = Pbins(:, ia) + accumarray(bins, wgt .* acc(:), [Lis.nbins 1]) * (dt * fs);
127 if kf <= nf && nstep == frame_steps(kf)
128 if free
129 [wxx, wyy, wxy] = plate_curvature_free(wfull, dx, dy, P.material.nu);
130 else
131 [wxx, wyy, wxy] = plate_curvature(plate_pad(wfull, P.bc), dx, dy);
132 end
133 S = plate_stress(P, wxx, wyy, wxy);
134 Wf(:, :, kf) = wfull; VMf(:, :, kf) = S.vm; SPf(:, :, kf) = S.sp; tf(kf) = nstep * dt;
135 kf = kf + 1;
136 end
138 peak = max(abs(wn));
139 if ~(peak < blowup)
140 diverged_step = nstep;
141 notes{end + 1} = sprintf('solution diverged at step %d (t = %.3g s); results are truncated there', nstep, nstep * dt);
142 break
143 end
144 wm = w; w = wn;
145end
147if diverged_step > 0
148 nkeep = diverged_step;
149 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);
150 keepf = frame_steps < diverged_step;
151 Wf = Wf(:, :, keepf); VMf = VMf(:, :, keepf); SPf = SPf(:, :, keepf); tf = tf(keepf);
152end
154p_audio = zeros(1, na);
155for d = 1:Lis.nbins
156 p_audio(d:end) = p_audio(d:end) + Pbins(d, 1:na - d + 1);
157end
158n_audio = min(na, ceil(numel(tt) * dt * fs));
159p_audio = p_audio(1:n_audio);
161Sp = plate_stress(P, wxx_p, wyy_p, wxy_p);
162[fw, Aw] = plate_spectrum(tt, w_p);
163[~, Avm] = plate_spectrum(tt, Sp.vm);
165result.method = 'implicit';
166result.grid = struct('x', x, 'y', y);
167result.t = tt;
168result.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);
169result.frames = struct('t', tf, 'w', Wf, 'vm', VMf, 'sp', SPf);
170if isempty(tf)
171 result.envelope = struct('vm', zeros(ny, nx), 'sp', zeros(ny, nx));
172else
173 result.envelope = struct('vm', max(VMf, [], 3), 'sp', max(abs(SPf), [], 3));
174end
175result.spectrum = struct('f', fw, 'w', Aw, 'vm', Avm);
176result.audio = struct('fs', fs, 'p', p_audio, 't0', Lis.t0);
177result.info = struct('nx', nx, 'ny', ny, 'dx', dx, 'dy', dy, 'dt', dt, 'theta', theta, ...
178 'steps_per_sample', P.implicit.steps_per_sample, 'nunknowns', m, 'nnz_factor', nnz(R), 'factor_s', factor_s, ...
179 'nsteps', nsteps, 'ncells', nx * ny, 'cell_steps', plan.cell_steps, 'fs', fs, ...
180 'stable', diverged_step == 0, 'diverged_step', diverged_step, 'wall_s', toc(t_start), 'notes', {notes});
181end