System of equations - unsymmetric matrix need help!
3 views (last 30 days)
Show older comments
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)
0 Comments
Answers (1)
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))
Hope it helps!
9 Comments
See Also
Categories
Find more on Logical in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!