1function [L, B, mx, my] = plate_matrices(nx, ny, dx, dy, bc)
2% PLATE_MATRICES Sparse Laplacian and biharmonic matrices on the interior nodes.
3%
4% [L, B, mx, my] = plate_matrices(nx, ny, dx, dy, bc) assembles, for the
5% mx = nx-2 by my = ny-2 interior nodes of the grid (edge nodes carry w = 0),
6% L the five-point Laplacian (5 nonzeros per row)
7% B the thirteen-point biharmonic (13 nonzeros per row)
8% as sparse matrices acting on w(:) with the interior field stored as an
9% my x mx array, so the y index runs fastest, matching w(2:end-1, 2:end-1)(:).
10%
11% The stencils are the same as the ghost-node scheme in plate_fd: the second
12% difference with w = 0 at the edge is D2; the fourth difference is D2*D2 for
13% simply supported edges (odd reflection, w_nn = 0), and for clamped edges the
14% even reflection w_ghost = w_1 changes the end coefficient of the fourth
15% difference from 5/dx^4 to 7/dx^4. The mixed term 2 w_xxyy never reaches a
16% ghost node. A test in the lab checks B against the explicit stencil.
17%
18% Everything is built from Kronecker products so the matrices stay sparse.
20mx = nx - 2;
21my = ny - 2;
22ex = ones(mx, 1);
23ey = ones(my, 1);
24D2x = spdiags([ex -2 * ex ex], -1:1, mx, mx) / dx^2;
25D2y = spdiags([ey -2 * ey ey], -1:1, my, my) / dy^2;
26Ix = speye(mx);
27Iy = speye(my);
29D4x = D2x * D2x;
30D4y = D2y * D2y;
31switch bc
32 case 'simply-supported'
33 % D2*D2 already encodes the odd reflection
34 case 'clamped'
35 D4x(1, 1) = 7 / dx^4; D4x(mx, mx) = 7 / dx^4;
36 D4y(1, 1) = 7 / dy^4; D4y(my, my) = 7 / dy^4;
37 otherwise
38 error('plate_matrices: unknown boundary condition "%s"', bc);
39end
41L = kron(Ix, D2y) + kron(D2x, Iy);
42B = kron(Ix, D4y) + 2 * kron(D2x, D2y) + kron(D4x, Iy);
43end