Matrix is singular to working precision

Attached is a pdf of the code that I have with a warning of matrix is singular to working precision that results in NaN for all values being solved for. Would love to know what is going on and how to fix it to produce numerical results

[Voss]

clc;
clear;
close all

% given data
f1= 79.7;
f2 = 14;
f3 = 80;
r4 = 51.26;
links= [f2 f3 r4 f1];

% Theta angles given
theta1= 126.582;
theta2 = 45;

%Initial guesses for Theta 3&4
theta3 = 30;
theta4 = 30;
Theta0 = [theta3 theta4];

%functions
myfunction = @(Theta0) linkeqnn(Theta0, links, theta1, theta2);

[angles, fval] = fsolve(myfunction, Theta0);

Unrecognized function or variable 'linkeqnn'.

Error in solution>@(Theta0)linkeqnn(Theta0,links,theta1,theta2) (line 22)
myfunction = @(Theta0) linkeqnn(Theta0, links, theta1, theta2);

Error in fsolve (line 167)
fuser = feval(funfcn{3},x,varargin{:});

fprintf('The numerical solution yields Theta 3 = %4.2f deg. and Theta 4 = %4.2f deg.\n', angles(1), angles(2))

%Knowns
omega2 = 10;
theta = [theta1, theta2, angles(1), angles(2)];
links = [79.7, 14, 80, 51.26];

%Initial Guesses
omega3 = 0;
omega4 = 0;
omega0 = [omega3 omega4];

%Call Solver
myv = @(omega0) veloeqn(omega0, theta, links, omega2);
[velos, vval] = fsolve(myv, omega0);

%Angular Acceleration Work

%Knowns
alpha2 = 0;
omega = [0, omega2, velos(1), velos(2)];

%Initial Guesses
alpha3 = 70;
alpha4 = 70;
alpha0 = [alpha3 alpha4];

%Call Solver
mya = @(alpha0) acceleqn(alpha0, omega, theta, links);
mya(alpha0)
[accels, aval] = fsolve(mya, alpha0);

% Matrix
f1 = 13.2;
f2 = links(2);
f3 = (1/2)*links(3);
f41 = 79.22;
f42 = links(4);
fp = 124.44;
f34 = 32;
R12X = f2*cosd(theta(2));
R12Y = f2*sind(theta(2));
R32X = (links(2)-f1)*cosd(theta(2));
R32Y = (links(2)-f1)*sind(theta(2));
R23X = f3*cosd(theta(3));
R23Y = f3*sind(theta(3));
R43X = f3*cosd(theta(3));
R43Y = f3*sind(theta(3));
R14X = f41*cosd(theta(4));
R14Y = f41*sind(theta(4));
R34X = f34*cosd(theta(4));
R34Y = f34*sind(theta(4));
RPX = fp*cosd(theta(4));
RPY = fp*sind(theta(4));
G = [1,0,1,0,0,0,0,0,0,0,0,0,0;
    0,1,0,1,0,0,0,0,0,0,0,0,0;
    R12Y,-R12X,-R32Y,R32X,1,0,0,0,0,0,0,0,0;
    0,0,0,0,0,1,0,1,0,0,0,0,0;
    0,0,0,0,0,0,1,0,1,0,0,0,0;
    0,0,0,0,0,R23Y,R23X,-R43Y,-R43X,0,0,0,0;
    0,0,0,0,0,0,0,0,0,1,0,1,0;
    0,0,0,0,0,0,0,0,0,0,1,0,1;
    0,0,0,0,0,0,0,0,0,R14Y,-R14X,-R34Y,-R34X;
    0,0,1,0,0,1,0,0,0,0,0,0,0;
    0,0,0,1,0,1,0,0,0,0,0,0,0;
    0,0,0,0,0,0,0,0,1,0,1,0,0;
    0,0,0,0,0,0,0,0,0,1,0,1,0;];
m2 = (598)/(32.2);
m3 = (193)/(32.2);
m4 = (2706)/(32.2);
ICG2 = 50;
ICG3 = 5;
ICG4 = 40;
ACG2X = -f2*(alpha2)*sind(theta2) - f2*(omega2^2)*cosd(theta2);
ACG2Y = f2*(alpha2)*cosd(theta2) - f2*(omega2^2)*sind(theta2);
ACG3X = -f2*(0)*sind(theta2) - f2*(10^2)*cosd(theta2) - f3*(accels(1)*sind(angle(1)) - f3*(omega(3)^2)*cosd(theta(3)));
ACG3Y = -f2*(0)*cosd(theta2) - f2*(10^2)*sind(theta2) - f3*(accels(1)*cosd(angle(1)) - f3*(omega(3)^2)*sind(theta(3)));
ACG4X = -f42*accels(2)*sind(angle(2)) - f42*(omega0(2)^2)*cosd(angle(2));
ACG4Y = f42*accels(2)*cosd(angle(2)) - f42*(omega0(2)^2)*sind(angle(2));
W2 = 598;
W3 = 193;
W4 = 2706;
Fcable = 2970;
A = [m2*ACG2X;
    m2*ACG2Y + W2;
    ICG2*alpha2;
    m3*ACG3X;
    m3*ACG3Y + W3;
    ICG3*alpha0(1);
    m4*ACG4X
    m4*ACG4Y + W4 + Fcable;
    ICG4*alpha0(2) - (124.44 + Fcable);
    0;
    0;
    0;
    0];

T = G\A;

% Display Numerical Solution for Theta 3 and Theta 4
fprintf('The numerical solution yields Theta 3 = %.2f degrees and Theta 4 = %.2f degrees.\n', angles(1), angles(2));

% Display Angular Velocities
fprintf('\nAngular Velocities:\n');
fprintf('Omega 3: %.2f rad/s\n', velos(1));
fprintf('Omega 4: %.2f rad/s\n', velos(2));

% Display Angular Accelerations
fprintf('\nAngular Accelerations:\n');
fprintf('Alpha 3: %.2f rad/s^2\n', accels(1));
fprintf('Alpha 4: %.2f rad/s^2\n', accels(2));

% Display Force Results from the Matrix Solution
fprintf('\nForce Results from the Matrix Solution:\n');
fprintf('Force 1: %.2f N\n', T(1));
fprintf('Force 2: %.2f N\n', T(2));
fprintf('Force 3: %.2f N\n', T(3));
fprintf('Force 4: %.2f N\n', T(4));
fprintf('Force 5: %.2f N\n', T(5));
fprintf('Force 6: %.2f N\n', T(6));
fprintf('Force 7: %.2f N\n', T(7));
fprintf('Force 8: %.2f N\n', T(8));
fprintf('Force 9: %.2f N\n', T(9));
fprintf('Force 10: %.2f N\n', T(10));
fprintf('Force 11: %.2f N\n', T(11));
fprintf('Force 12: %.2f N\n', T(12));
fprintf('Force 13: %.2f N\n', T(13));

Answer 1

[John D'Errico]

I don't even need to run your code, to see an obvious mistake. In fact, I knew what I expect to be a problem, even before I looked at your code.

When you set your parameter to have the SAME starting values, this creates a situation where the solver cannot essentially decide where where to look, which direction to go. That is, should I increase theta3 or theta4? Decrease one of the other? Which one? You might imagine fsolve as going into a funk, unable to make a decision.

Yes, you might just scream now, for gods sake, fsolve should just pick one! But it cannot. The mathematics generates a singular matrix. Now look at the errors you are getting. SINGULAR MATRIX.

Then I looked at your code.

%Initial guesses for Theta 3&4
theta3 = 30;
theta4 = 30;
Theta0 = [theta3 theta4];

%Initial Guesses
omega3 = 0;
omega4 = 0;
omega0 = [omega3 omega4];

Etc. You keep on doing exactly the same thing. Now, if I look at your code, at where those variables appear, everything is quite symmetrical. Again, this causes a singular matrix, when you start it out exactly there.

The answer is simple. NEVER do that. Never duplicate starting values in a problem like this.