finite element to calculate deflection of fixed beam and distributed force in matlab

main.m
% Define properties
L = 5; % Length of beam
E = 200e9; % Young's modulus
I = 1e-4; % Moment of inertia
w = 1000; % Distributed force

% Create finite element model
num_elements = 10; % Number of elements
num_nodes = num_elements + 1;
node_coords = linspace(0, L, num_nodes)';
elements = [1:num_elements; 2:num_nodes];

% Assemble global stiffness matrix
K = zeros(num_nodes, num_nodes);
for e = 1:num_elements
    element_nodes = elements(:, e);
    Le = node_coords(element_nodes(2)) - node_coords(element_nodes(1));
    ke = E*I/Le^3 * [12, 6*Le, -12, 6*Le; 
                     6*Le, 4*Le^2, -6*Le, 2*Le^2; 
                     -12, -6*Le, 12, -6*Le; 
                     6*Le, 2*Le^2, -6*Le, 4*Le^2];
    K(element_nodes, element_nodes) = K(element_nodes, element_nodes) + ke;
end

% Apply boundary conditions (fixed end)
fixed_node = 1;
K(fixed_node, :) = 0;
K(:, fixed_node) = 0;
K(fixed_node, fixed_node) = 1;

% Apply distributed force
f = zeros(num_nodes, 1);
for e = 1:num_elements
    element_nodes = elements(:, e);
    Le = node_coords(element_nodes(2)) - node_coords(element_nodes(1));
    fe = w*Le/2 * [1; 2; 1; 2];
    f(element_nodes) = f(element_nodes) + fe;
end

% Solve for displacements
u = K\f;

% Compute deflection at mid-span
mid_node = ceil(num_nodes/2);
deflection = u(mid_node);

disp(['Deflection at mid-span: ', num2str(deflection)]);
1369 chars
48 lines