Introduction to Matrix inversion and finite element analysis. Is this code right?
Mostrar comentarios más antiguos
clc;
clear;
close all;
%% =========================
% MATERIAL PROPERTIES
%% =========================
E = 300e9; % Young modulus (Pa)
d = 0.25e-2; % Diameter (m)
A = pi*d^2/4; % Area
%% =========================
% NODE COORDINATES
%% =========================
nodes = [0 0;
1 0;
1 1];
%% =========================
% ELEMENTS
%% =========================
elements = [1 3;
2 3];
%% =========================
% GLOBAL MATRIX
%% =========================
K = zeros(6);
for e = 1:2
n1 = elements(e,1);
n2 = elements(e,2);
x1 = nodes(n1,1);
y1 = nodes(n1,2);
x2 = nodes(n2,1);
y2 = nodes(n2,2);
L = sqrt((x2-x1)^2 + (y2-y1)^2);
c = (x2-x1)/L;
s = (y2-y1)/L;
% Element stiffness matrix
k = (E*A/L)* ...
[ c^2 c*s -c^2 -c*s
c*s s^2 -c*s -s^2
-c^2 -c*s c^2 c*s
-c*s -s^2 c*s s^2 ];
% DOF locations
if e == 1
dof = [1 2 5 6];
else
dof = [3 4 5 6];
end
% Assembly
K(dof,dof) = K(dof,dof) + k;
fprintf('Element %d stiffness matrix:\n',e)
disp(k)
end
%% =========================
% CHECKS BEFORE BC
%% =========================
disp('Global stiffness matrix:')
disp(K)
disp('Determinant:')
disp(det(K))
disp('Eigenvalues:')
disp(eig(K))
disp('Condition number:')
disp(cond(K))
%% =========================
% FORCE VECTOR
%% =========================
F = [0;
0;
0;
0;
100;
0];
%% =========================
% BOUNDARY CONDITIONS
%% =========================
fixed = [1 2 3 4];
free = [5 6];
Kff = K(free,free);
Ff = F(free);
%% =========================
% AFTER BC CHECKS
%% =========================
disp('Reduced stiffness matrix:')
disp(Kff)
disp('Determinant after BC:')
disp(det(Kff))
disp('Eigenvalues after BC:')
disp(eig(Kff))
disp('Condition number after BC:')
disp(cond(Kff))
%% =========================
% DISPLACEMENTS
%% =========================
Uf = inv(Kff)*Ff;
U = zeros(6,1);
U(free) = Uf;
disp('Displacements:')
disp(U)
%% =========================
% REACTION FORCES
%% =========================
R = K*U - F;
disp('Reaction forces:')
disp(R)
%% =========================
% PLOT
%% =========================
scale = 1000;
new_nodes = nodes;
new_nodes(1,1) = nodes(1,1) + scale*U(1);
new_nodes(1,2) = nodes(1,2) + scale*U(2);
new_nodes(2,1) = nodes(2,1) + scale*U(3);
new_nodes(2,2) = nodes(2,2) + scale*U(4);
new_nodes(3,1) = nodes(3,1) + scale*U(5);
new_nodes(3,2) = nodes(3,2) + scale*U(6);
figure
hold on
grid on
axis equal
% Original truss
for i = 1:2
n1 = elements(i,1);
n2 = elements(i,2);
plot([nodes(n1,1) nodes(n2,1)], ...
[nodes(n1,2) nodes(n2,2)], ...
'b','LineWidth',2)
end
% Deformed truss
for i = 1:2
n1 = elements(i,1);
n2 = elements(i,2);
plot([new_nodes(n1,1) new_nodes(n2,1)], ...
[new_nodes(n1,2) new_nodes(n2,2)], ...
'r--','LineWidth',2)
end
legend('Original','Deformed')
title('Truss Deformation')
xlabel('X')
ylabel('Y')
Respuestas (1)
Yes, the code looks correct. Here is a quick sanity-check of each key piece:
DOF numbering is consistent throughout: node i owns DOFs 2i-1 (x) and 2i (y)
NodeDOFs
(0,0) -> 1,2
(1,0) -> 3,4
(1,1) -> 5,6
Element DOF mapping matches that convention:
- Element 1 (nodes 1 -> 3): [1,2,5,6]
- Element 2 (nodes 2 -> 3): [3,4,5,6]
Boundary conditions pin nodes 1 and 2 (DOFs 1–4), leaving node 3 free (DOFs 5–6), correct for a two-bar truss with a tip load.
Applied force F(5)=100 puts a horizontal load at node 3 in the x-direction, which is consistent with the problem setup.
new_nodes reshape is correct and maps cleanly onto the nodes matrix layout.
Physics check element 2 is vertical (c=0, s=1), so it can only resist vertical displacement at node 3, not horizontal. The diagonal element 1 (c=s=1/sqrt(2)) is the only one bracing the horizontal load, which is why node 3 will show a large horizontal displacement and a small vertical one. That is the expected behaviour for this geometry.
Everything looks good.
I made a few tweaks to your code (e.g. fixed the legend, replaced INV with MLDIVIDE):
%% =========================
% MATERIAL PROPERTIES
%% =========================
E = 300e9; % Young modulus (Pa)
d = 0.25e-2; % Diameter (m)
A = pi*d^2/4; % Area
%% =========================
% NODE COORDINATES
%% =========================
nodes = [0,0;...
1,0;...
1,1];
%% =========================
% ELEMENTS
%% =========================
elements = [1,3;...
2,3];
%% =========================
% GLOBAL MATRIX
%% =========================
K = zeros(6);
for e = 1:2
n1 = elements(e,1);
n2 = elements(e,2);
x1 = nodes(n1,1);
y1 = nodes(n1,2);
x2 = nodes(n2,1);
y2 = nodes(n2,2);
L = sqrt((x2-x1)^2 + (y2-y1)^2);
c = (x2-x1)/L;
s = (y2-y1)/L;
% Element stiffness matrix
k = (E*A/L)* ...
[ c^2, c*s, -c^2, -c*s;...
c*s, s^2, -c*s, -s^2;...
-c^2, -c*s, c^2, c*s;...
-c*s, -s^2, c*s, s^2 ];
% DOF locations
if e == 1
dof = [1,2,5,6];
else
dof = [3,4,5,6];
end
% Assembly
K(dof,dof) = K(dof,dof) + k;
fprintf('Element %d stiffness matrix:\n',e)
disp(k)
end
%% =========================
% CHECKS BEFORE BC
%% =========================
disp('Global stiffness matrix:')
disp(K)
disp('Determinant:')
disp(det(K))
disp('Eigenvalues:')
disp(eig(K))
disp('Condition number:')
disp(cond(K))
%% =========================
% FORCE VECTOR
%% =========================
F = [0; 0; 0; 0; 100; 0];
%% =========================
% BOUNDARY CONDITIONS
%% =========================
fixed = [1,2,3,4];
free = [5,6];
Kff = K(free,free);
Ff = F(free);
%% =========================
% AFTER BC CHECKS
%% =========================
disp('Reduced stiffness matrix:')
disp(Kff)
disp('Determinant after BC:')
disp(det(Kff))
disp('Eigenvalues after BC:')
disp(eig(Kff))
disp('Condition number after BC:')
disp(cond(Kff))
%% =========================
% DISPLACEMENTS
%% =========================
%Uf = inv(Kff)*Ff;
Uf = Kff \ Ff; % better than INV
U = zeros(6,1);
U(free) = Uf;
disp('Displacements:')
disp(U)
%% =========================
% REACTION FORCES
%% =========================
R = K*U - F;
disp('Reaction forces:')
disp(R)
%% =========================
% PLOT
%% =========================
scale = 1000;
% new_nodes = nodes;
% new_nodes(1,1) = nodes(1,1) + scale*U(1);
% new_nodes(1,2) = nodes(1,2) + scale*U(2);
% new_nodes(2,1) = nodes(2,1) + scale*U(3);
% new_nodes(2,2) = nodes(2,2) + scale*U(4);
% new_nodes(3,1) = nodes(3,1) + scale*U(5);
% new_nodes(3,2) = nodes(3,2) + scale*U(6);
new_nodes = nodes + scale * reshape(U,2,[]).';
figure
hold on
grid on
axis equal
% Original truss
for i = 1:2
n1 = elements(i,1);
n2 = elements(i,2);
ho = plot([nodes(n1,1) nodes(n2,1)], ...
[nodes(n1,2) nodes(n2,2)], ...
'b','LineWidth',2);
end
% Deformed truss
for i = 1:2
n1 = elements(i,1);
n2 = elements(i,2);
hd = plot([new_nodes(n1,1) new_nodes(n2,1)], ...
[new_nodes(n1,2) new_nodes(n2,2)], ...
'r--','LineWidth',2);
end
legend([ho,hd],'Original','Deformed', 'location','northwest')
title('Truss Deformation')
xlabel('X')
ylabel('Y')
Categorías
Más información sobre Simulation Setup en Centro de ayuda y File Exchange.
el 7 de Mayo de 2026 a las 11:49
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
