FEA for cantilever beam

what is the problem in this code?

clear; clc;
% Define beam parameters
L = 2;              % Length of the beam
num_elements = 100; % Number of elements
num_modes = 6;      % Number of modes to calculate
% Prompt user for beam properties
rho = 2711.518;
A = 0.00555;
I = 0.00003186625;
E = 71*10^9;
% Compute element size
dx = L / num_elements;
% Assemble mass and stiffness matrices
M = zeros(num_elements+1, num_elements+1);
K = zeros(num_elements+1, num_elements+1);
for i = 1:num_elements
    % Element mass matrix
    m = [156    22*dx   54      -13*dx;
        22*dx  4*dx^2  13*dx   -3*dx^2;
        54     13*dx   156     -22*dx;
        -13*dx  -3*dx^2 -22*dx   4*dx^2] * (rho * A * dx / 420);
    
    % Element stiffness matrix
    k = [12     6*dx    -12     6*dx;
        6*dx    4*dx^2  -6*dx   2*dx^2;
        -12     -6*dx    12     -6*dx;
        6*dx     2*dx^2 -6*dx   4*dx^2] * (E * I / dx^3);
    % Assemble into global mass and stiffness matrices
    idx = (i-1)*2 + 1;
    M(idx:idx+3, idx:idx+3) = M(idx:idx+3, idx:idx+3) + m;
    K(idx:idx+3, idx:idx+3) = K(idx:idx+3, idx:idx+3) + k;
end
% Apply boundary conditions (cantilever beam)
M = M(2:end, 2:end);
K = K(2:end, 2:end);
% Solve eigenvalue problem
[V, D] = eig(K, M);
% Sort eigenvalues and eigenvectors
[D, sort_idx] = sort(diag(D));
V = V(:, sort_idx);
% Extract first num_modes eigenvalues and eigenvectors
omega = sqrt(D(1:num_modes));
modes = V(:, 1:num_modes);
% Plot mode shapes
x = linspace(0, L, num_elements+1);
figure;
hold on;
title('Cantilever Beam Mode Shapes');
xlabel('Beam Length');
ylabel('Displacement');
for i = 1:num_modes
    plot(x, [0; modes(:,i)], 'DisplayName', ['Mode ' num2str(i)]);
end
legend;
grid on;
% Display natural frequencies
disp('Natural Frequencies:');
disp(omega);

Answer 1

[Diwakar Diwakar]

Try this code:

clear; clc;
% Define beam parameters
L = 2;              % Length of the beam
num_elements = 100; % Number of elements
num_modes = 6;      % Number of modes to calculate

% Prompt user for beam properties
rho = 2711.518;
A = 0.00555;
I = 0.00003186625;
E = 71 * 10^9;

% Compute element size
dx = L / num_elements;

% Assemble mass and stiffness matrices
M = zeros((num_elements+1) * 2, (num_elements+1) * 2);
K = zeros((num_elements+1) * 2, (num_elements+1) * 2);

for i = 1:num_elements
    % Element mass matrix
    m = [156    22*dx   54      -13*dx;
        22*dx  4*dx^2  13*dx   -3*dx^2;
        54     13*dx   156     -22*dx;
        -13*dx  -3*dx^2 -22*dx   4*dx^2] * (rho * A * dx / 420);
    
    % Element stiffness matrix
    k = [12     6*dx    -12     6*dx;
        6*dx    4*dx^2  -6*dx   2*dx^2;
        -12     -6*dx    12     -6*dx;
        6*dx     2*dx^2 -6*dx   4*dx^2] * (E * I / dx^3);
    
    % Assemble into global mass and stiffness matrices
    idx = (i-1)*2 + 1;
    M(idx:idx+3, idx:idx+3) = M(idx:idx+3, idx:idx+3) + m;
    K(idx:idx+3, idx:idx+3) = K(idx:idx+3, idx:idx+3) + k;
end

% Apply boundary conditions (cantilever beam)
M = M(3:end, 3:end);
K = K(3:end, 3:end);

% Solve eigenvalue problem
[V, D] = eig(K, M);

% Sort eigenvalues and eigenvectors
[D, sort_idx] = sort(diag(D));
V = V(:, sort_idx);

% Extract first num_modes eigenvalues and eigenvectors
omega = sqrt(D(1:num_modes));
modes = V(:, 1:num_modes);

% Plot x-y axis
x = linspace(0, L, (num_elements+1) * 2);
y = zeros(size(x));

figure;
hold on;
title('Cantilever Beam Mode Shapes');
xlabel('Beam Length');
ylabel('Displacement');
plot(x, y, 'k-');
axis equal;
grid on;

% Display natural frequencies
disp('Natural Frequencies:');

Natural Frequencies:

disp(omega);

   1.0e+04 *

    0.0341
    0.2136
    0.5981
    1.1720
    1.9373
    2.8941

[Muhammad Saad]

Can you explain how you applied the boundary conditions?

'% Apply boundary conditions (cantilever beam)

M = M(3:end, 3:end);

K = K(3:end, 3:end);'