Introduction to Matrix inversion and finite element analysis. Is this code right?
Show older comments
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')
Answers (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')
Categories
Find more on Simulation Setup in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
