AYUDA EN CODIGO DE MATLAB PARA UN ESTRUCTURA DE PUNETE EN 3D

Quiero que me funcione para porder hallar los desplazamiento y reacciones, pero me sale el siguiente error. Warning: Matrix is singular to working

precision.

> In TrabajoF (line 80)

ADJUNTO EL CÓDIGO MATLAB:

% Armadura espacial
% Módulo de elasticidad:
E = 210 * 1e6; % en N/m² (Pa)

% Área:
A = 400 / 1e6; % en m²

% Coordenadas de los nodos:
XYZ = [16 5 0; 12 5 0;8 5 0; 4 5 0; 0 5 0; 12 5 4; 4 5 4; 0 0 0; 4 0 0;...
       8 0 0; 12 0 0; 16 0 0; 12 0 4; 4 0 4];

x = XYZ(:,1);
y = XYZ(:,2);
z = XYZ(:,3);

% Topología:
IJ = [1	2;2	3;3	4;4	5;1	6;2	6;3	6;3	7;4	7;5	7;6	7;12 11;11 10;10 9;...
      9	8;12 13;11 13;10 13;10 14;9	14;8 14;13 14;6	13;14 7;5 8;1 12];

% Dibujo de la estructura:
figure
for e = 1:26
    Q = [XYZ(IJ(e,1),1) XYZ(IJ(e,1),2) XYZ(IJ(e,1),3); ...
         XYZ(IJ(e,2),1) XYZ(IJ(e,2),2) XYZ(IJ(e,2),3)];
    plot3(Q(:,1), Q(:,2), Q(:,3), "-b")
    hold on
end
xlabel("x")
ylabel("y")
zlabel("z")
axis equal

% Longitudes de las barras
l = zeros(26,1);
for e = 1:26
    l(e) = sqrt((x(IJ(e,2))-x(IJ(e,1)))^2 + (y(IJ(e,2))-y(IJ(e,1)))^2 + (z(IJ(e,2))-z(IJ(e,1)))^2);
end

% Matriz de rigidez global
K = zeros(42,42);
for e = 1:26
    eta = (x(IJ(e,2)) - x(IJ(e,1))) / l(e);
    mu = (y(IJ(e,2)) - y(IJ(e,1))) / l(e);
    nu = (z(IJ(e,2)) - z(IJ(e,1))) / l(e);
    
    % Matriz de rigidez local del elemento
    K_e = (E * A / l(e)) * [eta^2 eta*mu eta*nu -eta^2 -eta*mu -eta*nu; ...
                             eta*mu mu^2 mu*nu -eta*mu -mu^2 -mu*nu; ...
                             eta*nu mu*nu nu^2 -eta*nu -mu*nu -nu^2; ...
                             -eta^2 -eta*mu -eta*nu eta^2 eta*mu eta*nu; ...
                             -eta*mu -mu^2 -mu*nu eta*mu mu^2 mu*nu; ...
                             -eta*nu -mu*nu -nu^2 eta*nu mu*nu nu^2];
    
    % Grados de libertad globales del elemento
    g_e = [3*IJ(e,1)-2 3*IJ(e,1)-1 3*IJ(e,1) 3*IJ(e,2)-2 3*IJ(e,2)-1 3*IJ(e,2)];
    
    % Expande K_e en la matriz de rigidez global K
    DeltaK_e = zeros(42,42);
    DeltaK_e(g_e, g_e) = K_e;
    K = K + DeltaK_e;
end

% Definición de desplazamientos y reacciones
% Grados de libertad restringidos y libres
a = (5:42)'; % Grados de libertad libres
b = (1:4)';  % Grados de libertad restringidos

% Dividir la matriz de rigidez
K_aa = K(a, a);
K_ab = K(a, b);
K_ba = K(b, a);
K_bb = K(b, b);

% Vector de fuerzas
P = zeros(42,1);
P(3) = -2; 
P(10) =-2% Fuerza en el nodo 1, dirección z

% Calcular desplazamientos en los grados de libertad libres
U_a = K_aa \ (P(a) - K_ab * zeros(length(b), 1));

% Reacciones en los grados de libertad restringidos
R_b = K_ba * U_a + K_bb * zeros(length(b), 1);

% Mostrar resultados
disp('Desplazamientos en los grados de libertad libres:')
disp(U_a)

disp('Reacciones en los grados de libertad restringidos:')
disp(R_b)

Answer 1

[Torsten]

If you insert the lines

% Calcular desplazamientos en los grados de libertad libres
rank(K_aa)
size(K_aa)

before the command

U_a = K_aa \ (P(a) - K_ab * zeros(length(b), 1));

you will see that K_aa is 38x38, but that its rank is only 25 - thus deeply rank-deficient.