System of equations - unsymmetric matrix need help!

4 views (last 30 days)
Leah Berry
Leah Berry on 12 Apr 2021
Commented: Leah Berry on 6 May 2021
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)

Answers (1)

Anshika Chaurasia
Anshika Chaurasia on 16 Apr 2021
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!
  9 Comments
Leah Berry
Leah Berry on 6 May 2021
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

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!