2 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)

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!

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

Start Hunting!