System of equations - unsymmetric matrix need help!

Hello,

I have a system of linear equations with 45 equations and 39 unkowns. I found a Mathworks Help Center article saying that "To solve a system of linear equations involving ill-conditioned (large condition number) non-square matrices, you must use QR decomposition." The condition number of my A matrix is 2.45e19, which I am assuming is considered large, so I believe I need to use this method. I am trying to use the "qrmatrixsolve" function to solve my X=A^-1*b system. I have created 39 symbolic variables, and when I try to run my script I keep getting an error saying that "qrmatrixsolve function does not accept arguments of type sym". I have tried changing my input matrices to "single" and "double", but then I get an error saying qrmatrixsolve cannot handle either of those either.

Does anyone know what I can do to get this to work? Code is included below:

% System of Equations from T4 Oil Victaulic 

%Variables 
f = 0.2;
R1 = 2.013;
R2 = 1.913;
R3 = 2;
R4 = 2.013;
R51 = 3.15;
R52 = 1.929;
R5o = 3.346; 
theta = 90;

syms A B C D E F G H I J K L M N O P Q R S T U V W X Y Z AA BB CC DD EE FF HH GG II JJ KK LL MM 

format short

%Pipe 1 Equations
eqn1 = (A) + (Y) + (CC) == 0;
eqn2 = (B) + (Z) - (DD) == -40;
eqn3 = (C) + (EE) == 40;
eqn4 = (-D) + (Z)*3.697 - (FF) - (DD)*(15.675) == -34920.77;
eqn5 = (E) - (Y)*3.697 - (GG) - (CC)*(15.675) == 0;
eqn6 = (CC)*(4.875) + (F) == 0;
% eqn7 = (-B)*(3.697) - (DD)*(9.978) + (D) - (FF) == -22228.99;
% eqn8 = (E) + (A)*(3.697) - (CC)*(9.978) - (GG) == 0;
% eqn9 = (CC)*4.875 + (F) == 0;
eqn10 = (-B)*(15.675) - (Z)*(9.978) - (FF) + (C)*(4.875) == 12252.9;
% eqn11 = (A)*(15.675) + (Y)*(9.978) + (E) - (GG) == 0;
eqn12 = (-Y)*(4.875) + (F) - (A)*(4.875) == 0;
eqn13 = (-GG) == (DD)*f*R1;
eqn14 = (E) == (C)*f*R1;

%Pipe 2 Equations
eqn15 = (G) - (2011.95) - (A/cos(11)) == -10544.31;
eqn16 = (-H) + (B) == 0;
eqn17 = (-C/cos(11)) - (I) + (A)/(cos(11)) == -1974.98;
eqn18 = (J) - (H)*(7.6053) + (D)/cos(11) + (L)/sin(11) == 0;
eqn19 = (-C)*(4.9969) + (K) + (-E) + (A)/cos(11)*7.6053 - (A)/sin(11)*6.5983 == -11079.61;
eqn20 = (L) + ((-F)*cos(11)) - (-B)*(6.5983) + (D)/sin(11) == 0;
% eqn21 = ((D)/cos(11)) - (-H)*(7.6053) + (J) - (-F)/sin(11) == 0;
eqn22 = (G)*(7.6053) + (I)*(6.5983) + (-E) + (K) == -15301.48;
% eqn23 = (-F)/cos(11) + (-B)*(6.5983) + (L) + (D)/sin(11) == 0;
eqn24 = (E) == (C)*f*R2;
eqn25 = (K) == (G)*f*R2;

%Pipe 3 Equations
eqn26 = (M) - (G) == 0;
eqn27 = (N) - (H) == 0;
eqn28 = (-O) + (I) == 0;
eqn29 = (-P) + (J) == 0;
eqn30 = (I)*(30.4035) + (Q) + (K) == 0;
eqn31 = (R) + (H)*(30.4035) + (L) == 0;
% eqn32 = (-P) + (J) == 0;
% eqn33 = (O)*(30.4035) + (K) + Q == 0;
% eqn34 = (-L) + (N)*(30.4035) + (R) == 0;
eqn35 = (J) == (G)*f*R3;
eqn36 = (P) == (M)*f*R3;

%Pipe 4 Equations
eqn37 = (-M)*sin(theta) + (S)*sin(theta) == -2227.8;
eqn38 = (-N) - (T) == -2227.8;
eqn39 = (U) + (O) == 0;
eqn40 = (V)*sin(theta) + (O)*(4.5) + (-P)*sin(theta) == 0;
eqn41 = (W) + (-Q) + (O)*(4.875)*sin(theta) == 0;
eqn42 = (X) + (M)*(4.5)*sin(theta) + (N)*(4.875) - (2227.8)*(4.5)*sin(theta) + (-R) == 0;
% eqn43 = (V)*sin(theta) - (U)*(4.5) + (-P)*sin(theta) == 0;
% eqn44 = (-Q) - (U)*(4.875) + (W) == 0;
eqn45 = (-R) + (S)*(4.5)*sin(theta) + (X) - (T)*(4.875) == -10860.53;
eqn46 = (P) == (M)*f*R4;
eqn47 = (V) == (S)*f*R4;

%Pipe 5 Equations
eqn48 = (-U)*sin(11) + (-S)*cos(11) + (HH) == 0;
eqn49 = (T)/cos(30) + (II) == -3792.94;
eqn50 = [(-U)*cos(11)]/cos(30) - [(-S)*sin(11)]/cos(30) + (JJ) == 0;
eqn51 = (KK) + (-V)*cos(11) + (-U)*cos(11)*(9.9698) - (T)*(3.6774) == -7502.58;
eqn52 = (LL) + (W)/cos(30) - [(-X)*cos(11)]/sin(30) - (-U)*sin(11)*(6.7824) - (-S)*cos(11)*(6.7824) == 0;
eqn53 = (MM) + [(-X)*cos(11)]/cos(30) - (T)/sin(30) - (-U)*sin(11)*(8.4821) - (-S)*cos(11)*(8.4821) == 0;
eqn54 = (KK) + (-V)*cos(11) + (-X)*sin(11) - (JJ)*(8.4821) - (II)*(6.7824) == 41746.83;
% eqn55 = (LL) + (W)/cos(30) + (-X)*cos(11)/sin(30) + (HH)*(6.7824) == 0;
% eqn56 = (MM) + (-X)*cos(11)/cos(30) - (W)/sin(30) + (HH)*(8.4821) == 0;
eqn57 = (W) == (T)*f*R51;
eqn58 = (LL) == (II)*f*R52; 

%Taking all equations and converting to matrix form
[Am,Bm] = equationsToMatrix([eqn1, eqn2, eqn3, eqn4, eqn5, eqn6, eqn10, eqn12, eqn13, eqn14, eqn15, eqn16, eqn17, eqn18, eqn19, eqn20, eqn22, eqn24, eqn25, eqn26, eqn27, eqn28, eqn29, eqn30, eqn31, eqn35, eqn36, eqn37, eqn38, eqn39, eqn40, eqn41, eqn42, eqn45, eqn46, eqn47, eqn48, eqn49, eqn50, eqn51, eqn52, eqn53, eqn54, eqn57, eqn58],[A, B, C, D, E, F, G, H, I, J, K, L, M, N, O, P, Q, R, S, T, U, V, W, X, Y, Z, AA, BB, CC, DD, EE, FF, GG, HH, II, JJ, KK, LL, MM]);

% Ad = vpa(Am)
% Bd = vpa(Bb)
% size(Ad)
% rank(Ad)

%Creating single floating point matrices??
Amd = single(Am);
Bmd = single(Bm);


%Creating symmetric matrix system 
Aam = (Amd.'*Amd);
Bam = (Amd.'*Bmd);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Choosing an iterative solver (A is a large, sparse matrix?)
%is A square?
size(Am,1) == size(Am,2)

%no, says to use least squares
size(Aam,1) == size(Aam,2)
%Aam matrix is square since we did an orthogonal transform
%now is matrix symmetric?
% tf = issymmetric(Am)
isequal(Am,Am.')
isequal(Aam,Aam.')
%Am is not symmetric, Aam is symmetric
%now is matrix positive definite?
% chol(Aam)
%says no not positive definite
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Least Squares Method
%large residual indicates preconditioner matrix is required
% [Xlsqr,flag,relres,iter,resvec,lsvec] = lsqr(Amd, Bmd, 1e-10, 10000)
% N = length(resvec);
% semilogy(0:N-1,lsvec,'--o',0:N-1,resvec,'-o')
% legend('Least-squares residual','Relative residual')%plot another variable here, once with big jump
%a large residual means that more iterations are needed


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%%%%other methods%%%%%


% Using qr() method
% [Cm, Rm,] = qr(Am, Bm)
% Cmm = vpa(Cm)
% Rmm = vpa(Rm)
% Xmm = Rmm\Cmm
rng default
Xsrs = qrmatrixsolve(Am, Bm)
% cond(Am)

Answer 1

[Anshika Chaurasia]

Hi,

You need to replace:

Xsrs = qrmatrixsolve(Am, Bm)

With following code:

Xsrs = fixed.qrMatrixSolve(double(Am), double(Bm))

Refer to fixed.qrMatrixSolve documentation for more information.

Hope it helps!

[Leah Berry]

Hello Anshika, thank you for this answer! I tried this code and am getting the following error

"Undefined variable "fixed" or class "fixed.qrmatrixsolve"."....please let me know if you have any additonal help on this error

Thank you,

Leah

[Leah Berry]

@Anshika Chaurasia tagging you in case this doesn't make it to you again! :)

[Anshika Chaurasia]

Hi Leah,

Which MATLAB version you are using?

[Walter Roberson]

fixed.qrMatrixSolve is case sensitive, and requires the Fixed Point Designer Toolbox

[Leah Berry]

@Walter Roberson Hi Walter, I have the fixed point designer toolbox installed. Please see snapshot attached. I also updated my function as you mentioned it is case sensitive and am still getting an error that says:

"Undefined variable "fixed" or class "fixed.qrMatrixSolve".

Error in T4OilVicRev1 (line 144)

Xsrs = fixed.qrMatrixSolve(double(Am), double(Bm))"

@Anshika Chaurasia Hi Anshika, I am using Matlab R2014b, is this version maybe not compatible with fixed.qrMatrixSolve or the Fixed Point Designer toolbox?

[Walter Roberson]

needs R2020b

[Leah Berry]

@Walter Roberson I am now getting these errors after using 2020b, any thoughts on this?

Error using fixed.cordicConstants>getNumberOfIterations (line 43)

Inputs to cordicConstants must be numeric

Error in fixed.cordicConstants (line 9)

niter = coder.const(getNumberOfIterations(x));

Error in fixed.qr.qrUpdate (line 12)

[niter, Kn] = fixed.cordicConstants(y);

Error in fixed.qrAB (line 31)

[C,R] = fixed.qr.qrUpdate(C,R,A(i,:),B(i,:));

Error in fixed.qrMatrixSolve (line 29)

[C,R] = fixed.qrAB(A,B,regularizationParameter);

Error in T4OilVicRev1 (line 144)

Xsrs = fixed.qrMatrixSolve(Am, Bm)

[Walter Roberson]

You took out the double() around Am and Bm .

[Leah Berry]

Thanks @Walter Roberson. I am now getting a solution but most of the numbers are NaN. I don't think this is possible as these are supposed to be forces and moments. Do you think it would be best to use an iterative solver instead?

Xsrs =

1.0e+05 *

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

Inf

0.0081

0.2419

0.0005

-2.5807

-0.0974

0.0055

-0.1104

-0.5242

-0.0350

8.8675

-0.2443