When I try to constrain node 3 to a line my displacements become ridiculously large, what am I doing wrong?

clear;
clc;
%convert inputs into points
n1 = [12, 0];
n2 = [12, 6];
n3 = [0, 0];
n4 = [0, 10];
Angr3 = 240;
Angf = 135;
% calculating lengths of members
L1 = sqrt((n2(1) - n1(1))^2 + (n2(2) - n1(2))^2);
L2 = sqrt((n3(1) - n1(1))^2 + (n3(2) - n1(2))^2);
L3 = sqrt((n3(1) - n2(1))^2 + (n3(2) - n2(2))^2);
L4 = sqrt((n4(1) - n2(1))^2 + (n4(2) - n2(2))^2);
L5 = sqrt((n4(1) - n3(1))^2 + (n4(2) - n3(2))^2);
% calculating angle of each member
Ang1 = rad2deg(atan2(n2(2) - n1(2), n2(1) - n1(1)));
Ang2 = rad2deg(atan2(n3(2) - n1(2), n3(1) - n1(1)));
Ang3 = rad2deg(atan2(n3(2) - n2(2), n3(1) - n2(1)));
Ang4 = rad2deg(atan2(n4(2) - n2(2), n4(1) - n2(1)));
Ang5 = rad2deg(atan2(n4(2) - n3(2), n4(1) - n3(1)));
% Material properties
Est = 30*10^6;
Eal = 11*10^6;
Ast = pi*(0.25)^2;
Aal = pi*(0.2)^2;
% Calculating K values
k1 = (Ast*Est)/L1;   % steel
k2 = (Aal*Eal)/L2;   % aluminum
k3 = (Ast*Est)/L3;   % steel
k4 = (Aal*Eal)/L4;   % aluminum
k5 = (Ast*Est)/L5;   % steel
% calculating normalized k matrices
l1 = cosd(Ang1);
m1 = sind(Ang1);
K1 = k1 * [l1^2, l1*m1, -l1^2, -l1*m1;
    l1*m1, m1^2, -l1*m1, -m1^2;
    -l1^2, -l1*m1, l1^2, l1*m1;
    -l1*m1, -m1^2, l1*m1, m1^2];
l2 = cosd(Ang2);
m2 = sind(Ang2);
K2 = k2 * [l2^2, l2*m2, -l2^2, -l2*m2;
    l2*m2, m2^2, -l2*m2, -m2^2;
    -l2^2, -l2*m2, l2^2, l2*m2;
    -l2*m2, -m2^2, l2*m2, m2^2];
l3 = cosd(Ang3);
m3 = sind(Ang3);
K3 = k3 * [l3^2, l3*m3, -l3^2, -l3*m3;
    l3*m3, m3^2, -l3*m3, -m3^2;
    -l3^2, -l3*m3, l3^2, l3*m3;
    -l3*m3, -m3^2, l3*m3, m3^2];
l4 = cosd(Ang4);
m4 = sind(Ang4);
K4 = k4 * [l4^2, l4*m4, -l4^2, -l4*m4;
    l4*m4, m4^2, -l4*m4, -m4^2;
    -l4^2, -l4*m4, l4^2, l4*m4;
    -l4*m4, -m4^2, l4*m4, m4^2];
l5 = cosd(Ang5);
m5 = sind(Ang5);
K5 = k5 * [l5^2, l5*m5, -l5^2, -l5*m5;
    l5*m5, m5^2, -l5*m5, -m5^2;
    -l5^2, -l5*m5, l5^2, l5*m5;
    -l5*m5, -m5^2, l5*m5, m5^2];
% Initialize global stiffness matrix
GlobalK = zeros(8, 8);
% L1: between nodes 1 and 2
GlobalK([1,2,3,4], [1,2,3,4]) = GlobalK([1,2,3,4], [1,2,3,4]) + K1;
% L2: between nodes 1 and 3
GlobalK([1,2,5,6], [1,2,5,6]) = GlobalK([1,2,5,6], [1,2,5,6]) + K2;
% L3: between nodes 2 and 3
GlobalK([3,4,5,6], [3,4,5,6]) = GlobalK([3,4,5,6], [3,4,5,6]) + K3;
% L4: between nodes 2 and 4
GlobalK([3,4,7,8], [3,4,7,8]) = GlobalK([3,4,7,8], [3,4,7,8]) + K4;
% L5: between nodes 3 and 4
GlobalK([5,6,7,8], [5,6,7,8]) = GlobalK([5,6,7,8], [5,6,7,8]) + K5;
% Retain a copy of the original stiffness matrix for reactions
GlobalK_original = GlobalK;
% Force vector
P = 2000;
F = zeros(8, 1);
F(1) = P*cosd(Angf);   
F(2) = P*sind(Angf);    

% Enforce constraint: d3y
tan_angr3 = tand(Angr3);  % Slope of the constraint line
% Modify GlobalK to enforce d3y = tan_angr3 * d3x
GlobalK(6, :) = GlobalK(6, :) - tan_angr3 * GlobalK(5, :);  
GlobalK(:, 6) = GlobalK(:, 6) - tan_angr3 * GlobalK(:, 5);  
GlobalK(6, 6) = GlobalK(6, 6) + tan_angr3^2 * GlobalK(5, 5); 
% Adjust force vector 
F(6) = F(6) - tan_angr3 * F(5);
% Remove (d3y) as it's now dependent
GlobalK(6, :) = 0;  % Zero out row 6
GlobalK(:, 6) = 0;  % Zero out column 6
GlobalK(6, 6) = 1;  % Stability term 
F(6) = 0; 
% Solve for displacements
D = GlobalK \ F;

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.250473e-19.

F = GlobalK_original * D;
disp(D);

   1.0e+15 *

   -1.1014
    0.1166
   -1.1597
    0.1166
   -1.1014
         0
   -1.1986
    0.0000

disp(F);

   1.0e+04 *

         0
         0
   -4.8922
   -1.6461
    6.5536
   -0.4125
    1.6154
   -0.3991

Answer 1

[Walter Roberson]

The rank() of GlobalK is only 6 but the matrix is 8 x 8. Using GlobalK \ F is likely to generate garbage.

If you must use a rank 6 matrix GlobalK then you should calculate

D = pinv(GlobalK) * F;

[Emerson Hale]

It works!!!!! there is still definatley some rounding error but it works thank you Walter