Linear Equation Solving, I am getting the wrong output

I have two pretty easy equations to solve for two matricies within these equations. They are easy to solve by hand but I am having trouble finding the best approach to solving it right now I tried to limit the outputs but that just gave me some wierd values. Any ideas are appreciated.

clear;
clc;
close all
spring_Const = [500,600,400,400];
nodal_conn = [1,2;2,3;3,4;2,4];
disp = [0;NaN;NaN;0];
force = [NaN;NaN;1000;NaN];
Num_node = unique(nodal_conn);
Num_Elements = size(nodal_conn,1);
Global_K = zeros(length(Num_node),length(Num_node));
for count = 1:Num_Elements
    node_1_pos = nodal_conn(count,1);
    node_2_pos = nodal_conn(count,2);
    Global_K(node_1_pos,node_1_pos) = Global_K(node_1_pos,node_1_pos) + spring_Const(count);
    Global_K(node_2_pos,node_2_pos) = Global_K(node_2_pos,node_2_pos) + spring_Const(count);
    Global_K(node_2_pos,node_1_pos) = Global_K(node_2_pos,node_1_pos) - spring_Const(count);
    Global_K(node_1_pos,node_2_pos) = Global_K(node_1_pos,node_2_pos) - spring_Const(count);
end 

% Causing Issue %
syms 'F' 'U' [1 4] real
U(1)=0;
U(4)=0;
F(3)=1000;
F = transpose(F);
U = transpose(U);
EQN_1 = F == Global_K*U %% equations set up properly but not solving how I want 

EQN_1 = 

EQN_2 = sum(F) == 0

EQN_2 = 

solve([EQN_1;EQN_2]) % Answers should be as follows F1 = -263, F2 = around 0, F4 = -736, U2 = .526, U3 = 1.316

ans = struct with fields:
    F1: 2500/3
    F2: -2500
    F4: 2000/3
    U2: -5/3
    U3: 0

Answer 1

[Matt J]

Your equations system is singular, so no surprise that there is more than one solution.

[A,b]=equationsToMatrix([EQN_1;EQN_2]) 

A = 

b = 

det(A)

ans = 

0

Answer 2

[Animesh]

Hi @Jonathan Leutz

The issue with the code is that you are trying to solve a system of linear equations using symbolic variables, but you are not specifying which variables you want to solve for. In this case, you have more unknowns than equations, which leads to an underdetermined system. Additionally, you are not using the known force at node 3 (1000 N) effectively in the symbolic system.

The correct approach is to reduce the system to only the unknowns and then solve for them. Since the displacements at nodes 1 and 4 are fixed (0), you only need to solve for the displacements at nodes 2 and 3 (U(2) and U(3)). Similarly, the force at node 3 is known (1000 N), so you only need to solve for the forces at nodes 1, 2, and 4 (F(1), F(2), and F(4)).

Here's how you can modify the code to solve the system correctly:

clear;
clc;
close all;
spring_Const = [500,600,400,400];
nodal_conn = [1,2;2,3;3,4;2,4];
Num_node = unique(nodal_conn);
Num_Elements = size(nodal_conn,1);
Global_K = zeros(length(Num_node),length(Num_node));
for count = 1:Num_Elements
    node_1_pos = nodal_conn(count,1);
    node_2_pos = nodal_conn(count,2);
    Global_K(node_1_pos,node_1_pos) = Global_K(node_1_pos,node_1_pos) + spring_Const(count);
    Global_K(node_2_pos,node_2_pos) = Global_K(node_2_pos,node_2_pos) + spring_Const(count);
    Global_K(node_2_pos,node_1_pos) = Global_K(node_2_pos,node_1_pos) - spring_Const(count);
    Global_K(node_1_pos,node_2_pos) = Global_K(node_1_pos,node_2_pos) - spring_Const(count);
end

% Known displacements and forces
U_known = [0; NaN; NaN; 0];
F_known = [NaN; NaN; 1000; NaN];

% Separate the knowns and unknowns
known_disp_indices = ~isnan(U_known);
unknown_force_indices = isnan(F_known);

% Reduced stiffness matrix for unknown displacements
K_uu = Global_K(~known_disp_indices, ~known_disp_indices);

% Solve for unknown displacements
F_ext = zeros(length(Num_node),1);
F_ext(3) = 1000; % External force applied at node 3
U_unknown = K_uu \ F_ext(~known_disp_indices);

% Display results for displacements
disp('Displacements at nodes 2 and 3:');

Displacements at nodes 2 and 3:

disp(U_unknown);

    0.5263
    1.3158

% Update the displacement vector with the calculated unknowns
U_known(~known_disp_indices) = U_unknown;

% Calculate the reaction forces at the supports
F_reaction = Global_K * U_known;

% Display results for forces
disp('Forces at nodes 1, 2, and 4:');

Forces at nodes 1, 2, and 4:

disp(F_reaction(unknown_force_indices));

 -263.1579
    0.0000
 -736.8421