matlab code to calculate analytical and fem modes for the rectangular plate simlpy supported on 4 edges in matlab

main.m
L = 1;  % length of the plate in x-direction
W = 0.5; % length of the plate in y-direction
h = 0.01; % thickness of the plate
E = 2e11; % Young's modulus of the material
nu = 0.3; % Poisson's ratio of the material
rho = 7850; % density of the material
252 chars
7 lines

main.m
nx = 20; % number of elements in x-direction
ny = 10; % number of elements in y-direction

x = linspace(0, L, nx + 1);
y = linspace(0, W, ny + 1);
[X, Y] = meshgrid(x, y);

mesh = zeros(3, 2, nx*ny);
m = 0;

for i = 1:nx
    for j = 1:ny
        m = m + 1;
        n1 = (ny + 1)*(i - 1) + j;
        n2 = (ny + 1)*i + j;
        mesh(:, :, m) = [n1, n2; n1 + 1, n2 + 1; n1, n2 + 1];
    end
end
395 chars
19 lines

main.m
K = zeros(2*(nx+1)*(ny+1), 2*(nx+1)*(ny+1));
M = zeros(2*(nx+1)*(ny+1), 2*(nx+1)*(ny+1));

for i = 1:size(mesh, 3)
    nodes = mesh(:, :, i);
    x1 = X(nodes(1, :), nodes(2, :));
    y1 = Y(nodes(1, :), nodes(2, :));
    [Ke, Me] = plate_element_stiffness(x1, y1, h, E, nu, rho);
    idx = [2*nodes(1, :) - 1, 2*nodes(1, :), 2*nodes(2, :) - 1, 2*nodes(2, :)];
    K(idx, idx) = K(idx, idx) + Ke;
    M(idx, idx) = M(idx, idx) + Me;
end
437 chars
13 lines

main.m
[V, D] = eig(K, M);

frequencies = sqrt(diag(D))/2/pi; % convert from angular frequency to Hz

modes = zeros(size(X(:), 1), size(V, 2));

for i = 1:size(V, 2)
    modes(:, i) = V(:, i);
end
190 chars
10 lines

main.m
num_modes = 3;

figure;
for i = 1:num_modes
    mode = modes(:, i);
    mode_x = [mode(1:2:end); 0];
    mode_y = [mode(2:2:end); 0];
    subplot(1, num_modes, i);
    surf(X, Y, reshape(mode, ny+1, nx+1));
    xlabel('x [m]');
    ylabel('y [m]');
    zlabel('Amplitude');
    title(['Mode shape ', num2str(i), ', f = ', num2str(frequencies(i)), ' Hz']);
    xlim([0 L]);
    ylim([0 W]);
    view(-120, 30);
    grid on;
end
427 chars
19 lines