Solution of irradical equation in matlab

I have a 4 times 4 matrix. The matrix contains constant and a variable x. Some terms of matrix contains squra root terms of x. Determinant of this matrix is not polynomial but an equation containg radicals in it. i want to solve determinant for x. I applied vpasolve(det), this gives me a only single value of x.

As the matrix is very complex so it is not possible to manually square it and then reduce to polynomial.

How can i solve this problem??

[Walter Roberson]

Does this happen to be effectively a process of finding eigenvalues ?

Could you post your matrix for us to experiment with?

Do you need to find the indefinitely precise solutions (formulas), or is it good enough to find all of the numeric approximations? Do you need to find all roots (complex) or only the real roots?

How many roots do you need? With the kind of setup you describe, you could quite easily end up with something equivalent to a polynomial of degree 2^(number of terms that include sqrt)

[Himanshu SINGLA]

clc;
clear all;
format long;
RR=[];
X1=[];
E=[];
S1=[]; 
OMGG=[];
syms K Y1 Y2 Y3 Y4 Y5
 OMG2=2:2:10;
%  OMG=8.2:.5:20;
  for  OMG1=1:size(OMG2,2); 
       OMG=OMG2(1,OMG1);
       OMGG=[OMGG;OMG];
OMGS=OMG*OMG;


AT0=0.005;
OMEGA=1;
AT1=0;


ET=1i;
T0=293; %%in K
ALM1 =7.59 ;
ALM2=2.14;
MU1=1.89 ;
MU2=0.45 ;
RHO1=2192;
RHO2=1010;
C0=96.3;   %specific heat
KS =2.51;  %classical fourier constant
S=0.021;   %heat generation due to velocity diff
XI=750;    %momentum generatio coeff due to vel diff
BETA0 =0.00005 ;  % elastic constant of isotropic solid
DELTA=0.0001;  %helmholtz free energy function
N=0.15;   %porosity of mixture
RHO= (1-N)*RHO1 + N*RHO2;


C11=(ALM1+2*MU1)./RHO1;
C12= MU1./RHO1;
C21=(ALM2+2*MU2)./RHO2;
C22= MU2./RHO2;
C12D=ET*OMG*C12;
C22D=ET*OMG*C22;
XI1B= ET*XI./OMG*RHO1;  
XI2B= ET*XI./OMG*RHO2;
T0S=AT0+(ET./OMG);    % (tau0 star)
T1S=AT1+(ET./OMG);    %(tau1 star)
BETA1D= S./T0 + BETA0;
BETA1= BETA1D*ET*T1S./OMG*RHO1;
BETA2= BETA1D*ET*T1S./OMG*RHO1;



SB= S./RHO*C0; %S BAR 
KB=KS./RHO*C0*T0;
BETA0B=BETA0./RHO*C0;
DELTAB=RHO*DELTA/C0;

KD= KB./(1-ET*OMG*T0*OMEGA)*(ET*OMG);
TAUM = T0S./(1-ET*OMG*T0*OMG)*(ET*OMG);
SD= SB./(1-ET*OMG*T0*OMEGA); %S DASH
A= SD - BETA0;
B= SD + DELTAB; 




A0=(-C11*C12D*KD);
A1=(C11*KD*OMGS - C12D*KD*OMGS - A*BETA1*C12D*OMGS + B*BETA2*C11*OMGS - C11*C12D*OMGS*TAUM + C11*KD*OMGS*XI2B - C12D*KD*OMGS*XI1B);
A2=(KD*OMGS^2 + A*BETA1*OMGS^2 + B*BETA2*OMGS^2 + C11*OMGS^2*TAUM - C12D*OMGS^2*TAUM + KD*OMGS^2*XI1B + KD*OMGS^2*XI2B + A*BETA1*OMGS^2*XI2B - A*BETA2*OMGS^2*XI1B - B*BETA1*OMGS^2*XI2B + B*BETA2*OMGS^2*XI1B + C11*OMGS^2*TAUM*XI2B - C12D*OMGS^2*TAUM*XI1B); 
A3=  OMGS^3*TAUM + OMGS^3*TAUM*XI1B + OMGS^3*TAUM*XI2B;
Y=[A0 A1 A2 A3];

r=roots(Y);
K1S=r(1);
K2S=r(2);
K3S=r(3);

B0= C21*C22D;
B1= C21*XI2B+C21-C22D*XI1B-C22D;
B2=(1+XI1B+XI2B);
Z=[B0 B1 B2];
R=roots(Z);
K4S=R(1);
K5S=R(2);
KSQ=K*K;
%  M1 = [C11*P +(XI1B + 1)*OMGS -XI1B*OMGS BETA1*OMGS; XI2B*OMGS -C21D*P + (XI2B + 1)*OMGS -BETA2*OMGS; -A*P B*P KD*P +OMGS*TAUM]
%  M2= [C21*P +(XI1B + 1)*OMGS -XI1B*OMGS; XI2B*OMGS -C21D*P + (XI2B + 1)*OMGS];
%  
% K1S= ;
% K2S= ;
% K3S= ;
% K4S= ;
% K5S= ;

 M1= sqrt(KSQ-K1S);
  M2= sqrt(KSQ-K2S);
  M3= sqrt(KSQ-K3S);
  M4= sqrt(KSQ-K4S);
  M5= sqrt(KSQ-K5S);
  M1S=M1*M1;
  M2S=M2*M2;
  M3S=M3*M3;
  M4S=M4*M4;
  M5S=M5*M5;
 
  
  ZETA1=(C11*K1S+(XI1B+1)*OMGS)*BETA2*OMGS+XI2B*BETA1*OMGS*OMGS./XI1B*BETA2*OMGS*OMGS+(-C12D*K1S+ (XI2B+1)*OMGS)*BETA1*OMGS;
  ZETA2=(C11*K2S+(XI1B+1)*OMGS)*BETA2*OMGS+XI2B*BETA1*OMGS*OMGS./XI1B*BETA2*OMGS*OMGS+(-C12D*K2S+ (XI2B+1)*OMGS)*BETA1*OMGS;
  ZETA3=(C11*K3S+(XI1B+1)*OMGS)*BETA2*OMGS+XI2B*BETA1*OMGS*OMGS./XI1B*BETA2*OMGS*OMGS+(-C12D*K3S+ (XI2B+1)*OMGS)*BETA1*OMGS;
  
 ETA1=XI1B*OMGS*ZETA1-(C11*K1S+(XI1B+1)*OMGS);
 ETA2=XI1B*OMGS*ZETA2-(C11*K2S+(XI1B+1)*OMGS);
 ETA3=XI1B*OMGS*ZETA3-(C11*K3S+(XI1B+1)*OMGS);
 
 XI4= C21*K4S+ (XI1B+1)*OMGS;
 XI5= C21*K5S+ (XI1B+1)*OMGS;
 
  
  
  
  
  
  
  
  A11= ALM1*K1S+2*MU1*M1S+BETA0*T1S*ETA1;
  A12= ALM1*K2S+2*MU1*M2S+BETA0*T1S*ETA2;
  A13= ALM1*K3S+2*MU1*M3S+BETA0*T1S*ETA3;
  A14=-2*ET*MU1*K*M4S;
  A15=-2*ET*MU1*K*M5S;
  
  A21= 2*ET*K*M1S;
  A22= 2*ET*K*M2S;
  A23= 2*ET*K*M3S;
  A24=KSQ-M4S;
  A25=KSQ-M5S;
  
  A31=ETA1*M1S;
  A32=ETA2*M2S;
  A33=ETA3*M3S;
  A34=0;
  A35=0;
  
  A41= (ALM2*K1S+2*MU2*M1S)*ZETA1;
  A42= (ALM2*K2S+2*MU2*M2S)*ZETA2;
  A43= (ALM2*K3S+2*MU2*M3S)*ZETA3;
  A44=-2*ET*MU2*K*M4S*XI4;  
  A45=-2*ET*MU2*K*M5S*XI5;    
    
  A51=2*ET*ZETA1*K*M1S;
  A52=2*ET*ZETA2*K*M2S;
  A53=2*ET*ZETA3*K*M3S;
  A54=(KSQ-M4S)*XI4;
  A55=(KSQ-M5S)*XI5;  
  
  M=[A11 A12 A13 A14 A15; A21 A22 A23 A24 A25; A31 A32 A33 A34 A35; A41 A42 A43 A44 A45; A51 A52 A53 A54 A55];
  
  D=det(M);
   Y= collect(D,K);
  
%   RR=root(D);%
   RR=vpasolve(D==0,K);%
   
   
By running this, I only get single value of RR rather It should give more than 1 values.

Answer 1

[Walter Roberson]

No problem.

format long;

Q = @(v) sym(v);

RR=[];
X1=[];
E=[];
S1=[]; 
OMGG=[];
syms K Y1 Y2 Y3 Y4 Y5
 OMG2=Q(2:2:10);
%  OMG=8.2:.5:20;
  for  OMG1=1:size(OMG2,2) 
       OMG=OMG2(1,OMG1);
       OMGG=[OMGG;OMG];
  end
OMGS=OMG*OMG;
AT0=Q(0.005);
OMEGA=Q(1);
AT1=Q(0);
ET=Q(1i);
T0=Q(293); %%in K
ALM1 =Q(7.59) ;
ALM2=Q(2.14);
MU1=Q(1.89) ;
MU2=Q(0.45) ;
RHO1=Q(2192);
RHO2=Q(1010);
C0=Q(96.3);   %specific heat
KS =Q(2.51);  %classical fourier constant
S=Q(0.021);   %heat generation due to velocity diff
XI=Q(750);    %momentum generatio coeff due to vel diff
BETA0 =Q(0.00005) ;  % elastic constant of isotropic solid
DELTA=Q(0.0001);  %helmholtz free energy function
N=Q(0.15);   %porosity of mixture
RHO= (1-N)*RHO1 + N*RHO2;
C11=(ALM1+2*MU1)./RHO1;
C12= MU1./RHO1;
C21=(ALM2+2*MU2)./RHO2;
C22= MU2./RHO2;
C12D=ET*OMG*C12;
C22D=ET*OMG*C22;
XI1B= ET*XI./OMG*RHO1;  
XI2B= ET*XI./OMG*RHO2;
T0S=AT0+(ET./OMG);    % (tau0 star)
T1S=AT1+(ET./OMG);    %(tau1 star)
BETA1D= S./T0 + BETA0;
BETA1= BETA1D*ET*T1S./OMG*RHO1;
BETA2= BETA1D*ET*T1S./OMG*RHO1;
SB= S./RHO*C0; %S BAR 
KB=KS./RHO*C0*T0;
BETA0B=BETA0./RHO*C0;
DELTAB=RHO*DELTA/C0;
KD= KB./(1-ET*OMG*T0*OMEGA)*(ET*OMG);
TAUM = T0S./(1-ET*OMG*T0*OMG)*(ET*OMG);
SD= SB./(1-ET*OMG*T0*OMEGA); %S DASH
A= SD - BETA0;
B= SD + DELTAB; 
A0=(-C11*C12D*KD);
A1=(C11*KD*OMGS - C12D*KD*OMGS - A*BETA1*C12D*OMGS + B*BETA2*C11*OMGS - C11*C12D*OMGS*TAUM + C11*KD*OMGS*XI2B - C12D*KD*OMGS*XI1B);
A2=(KD*OMGS^2 + A*BETA1*OMGS^2 + B*BETA2*OMGS^2 + C11*OMGS^2*TAUM - C12D*OMGS^2*TAUM + KD*OMGS^2*XI1B + KD*OMGS^2*XI2B + A*BETA1*OMGS^2*XI2B - A*BETA2*OMGS^2*XI1B - B*BETA1*OMGS^2*XI2B + B*BETA2*OMGS^2*XI1B + C11*OMGS^2*TAUM*XI2B - C12D*OMGS^2*TAUM*XI1B); 
A3=  OMGS^3*TAUM + OMGS^3*TAUM*XI1B + OMGS^3*TAUM*XI2B;
Y=[A0 A1 A2 A3];
%r=roots(Y);
r = solve(poly2sym(Y), 'MaxDegree', 3);
K1S=r(1);
K2S=r(2);
K3S=r(3);
B0= C21*C22D;
B1= C21*XI2B+C21-C22D*XI1B-C22D;
B2=(1+XI1B+XI2B);
Z=[B0 B1 B2];
R = solve(poly2sym(Z), 'MaxDegree', 3);
K4S=R(1);
K5S=R(2);
KSQ=K*K;
%  M1 = [C11*P +(XI1B + 1)*OMGS -XI1B*OMGS BETA1*OMGS; XI2B*OMGS -C21D*P + (XI2B + 1)*OMGS -BETA2*OMGS; -A*P B*P KD*P +OMGS*TAUM]
%  M2= [C21*P +(XI1B + 1)*OMGS -XI1B*OMGS; XI2B*OMGS -C21D*P + (XI2B + 1)*OMGS];
%  
% K1S= ;
% K2S= ;
% K3S= ;
% K4S= ;
% K5S= ;
 M1= sqrt(KSQ-K1S);
  M2= sqrt(KSQ-K2S);
  M3= sqrt(KSQ-K3S);
  M4= sqrt(KSQ-K4S);
  M5= sqrt(KSQ-K5S);
  M1S=M1*M1;
  M2S=M2*M2;
  M3S=M3*M3;
  M4S=M4*M4;
  M5S=M5*M5;
 
  
  ZETA1=(C11*K1S+(XI1B+1)*OMGS)*BETA2*OMGS+XI2B*BETA1*OMGS*OMGS./XI1B*BETA2*OMGS*OMGS+(-C12D*K1S+ (XI2B+1)*OMGS)*BETA1*OMGS;
  ZETA2=(C11*K2S+(XI1B+1)*OMGS)*BETA2*OMGS+XI2B*BETA1*OMGS*OMGS./XI1B*BETA2*OMGS*OMGS+(-C12D*K2S+ (XI2B+1)*OMGS)*BETA1*OMGS;
  ZETA3=(C11*K3S+(XI1B+1)*OMGS)*BETA2*OMGS+XI2B*BETA1*OMGS*OMGS./XI1B*BETA2*OMGS*OMGS+(-C12D*K3S+ (XI2B+1)*OMGS)*BETA1*OMGS;
  
 ETA1=XI1B*OMGS*ZETA1-(C11*K1S+(XI1B+1)*OMGS);
 ETA2=XI1B*OMGS*ZETA2-(C11*K2S+(XI1B+1)*OMGS);
 ETA3=XI1B*OMGS*ZETA3-(C11*K3S+(XI1B+1)*OMGS);
 
 XI4= C21*K4S+ (XI1B+1)*OMGS;
 XI5= C21*K5S+ (XI1B+1)*OMGS;
 
  
  
  
  
  
  
  
  A11= ALM1*K1S+2*MU1*M1S+BETA0*T1S*ETA1;
  A12= ALM1*K2S+2*MU1*M2S+BETA0*T1S*ETA2;
  A13= ALM1*K3S+2*MU1*M3S+BETA0*T1S*ETA3;
  A14=-2*ET*MU1*K*M4S;
  A15=-2*ET*MU1*K*M5S;
  
  A21= 2*ET*K*M1S;
  A22= 2*ET*K*M2S;
  A23= 2*ET*K*M3S;
  A24=KSQ-M4S;
  A25=KSQ-M5S;
  
  A31=ETA1*M1S;
  A32=ETA2*M2S;
  A33=ETA3*M3S;
  A34=0;
  A35=0;
  
  A41= (ALM2*K1S+2*MU2*M1S)*ZETA1;
  A42= (ALM2*K2S+2*MU2*M2S)*ZETA2;
  A43= (ALM2*K3S+2*MU2*M3S)*ZETA3;
  A44=-2*ET*MU2*K*M4S*XI4;  
  A45=-2*ET*MU2*K*M5S*XI5;    
    
  A51=2*ET*ZETA1*K*M1S;
  A52=2*ET*ZETA2*K*M2S;
  A53=2*ET*ZETA3*K*M3S;
  A54=(KSQ-M4S)*XI4;
  A55=(KSQ-M5S)*XI5;  
  
  M=[A11 A12 A13 A14 A15; A21 A22 A23 A24 A25; A31 A32 A33 A34 A35; A41 A42 A43 A44 A45; A51 A52 A53 A54 A55];
  
  D=det(M);
   Y= collect(D,K);
  
%   RR=root(D);%
   RRR = solve(D==0, K, 'MaxDegree', 4);
   RR = vpa(RRR);%
   RRd = double(RR);
   RRd

RRd = 8×1

-0.488289840687524 + 0.482618182235875i -0.000000377515847 - 0.000001826506960i -0.000001826150043 + 0.000000377239853i -8.393573674339343 + 2.382886985058038i 0.488289840687524 - 0.482618182235875i 0.000000377515847 + 0.000001826506960i 0.000001826150043 - 0.000000377239853i 8.393573674339343 - 2.382886985058038i

●

[Walter Roberson]

The RRR variable will contain exact solutions.

The purpose of Q() is to convert floating point values into symbolic rationals.

[Himanshu SINGLA]

Okay, thankyou so much,

But why did you consider maximum degree 4 in RRR, degree is going to be more than 4.

[Himanshu SINGLA]

And how does VPA in next line working, instead of vpasolve in one single line?

[Himanshu SINGLA]

Will all the eight values satisfy determinant D?

[Walter Roberson]

If you look at

 solve(D==0, K)

you will see that there are 8 solutions, 4 of which are exact negative of the other 4. The groups of 4 are square roots of the four roots of polynomials of degree 4, such as

root(z^4 - z^3*(48145969404149762426206552992398560878657588144436438822799170711243321317125/7432292235941762209320847890919342320851111829324966121434199809438 - 20053905075213392701881239662596138077802525426623354113383688320619073082125i/4954861490627841472880565260612894880567407886216644080956133206292) - z^2*(20069142256998730275265101549258697051935917553714964060756313010805902439534235875/108500616583091419114172961911231274756950537654379067466192697948 + 2780256691066021693141105012180147846819732328291329018364913472994851228074256750i/9041718048590951592847746825935939563079211471198255622182724829) - z*(163097489106001911469345509927015924193085777780099370400967081046971103875/2739146000695155040750509557366287441269438538259219094518188728 + 1147329233840373042063728693689625035193499503163556429241550342559527625i/114131083362298126697937898223595310052893272427467462271591197) + (14521641018413932313181557563939213803237784799771426193074919280177000000/342393250086894380093813694670785930158679817282402386814773591 + 3158847732808475138233191505945938604400903363157835056636521134515200000i/342393250086894380093813694670785930158679817282402386814773591), z, 1)^(1/2)

[Walter Roberson]

RRR is all 8 exact solutions (all of which will satisfy the determinent being 0)

vpa(RRR) is numeric approximations of the values, in the same form you would get from vpasolve()

double(vpa(RRR)) gives the double precision values.

vpasolve() has two completely different modes of operations. If it detects that the input is a polynomial of a single variable, then vpasolve() in the background does a solve() and then vpa() so that it gives all of the outputs instead of just one output. But if the system does not look like a polynomial in one variable, then vpasolve() uses a modified Newton Raphson numeric method to find a single root.

If you were to vpasolve(D==0, K) at that point in the code, you would see all 8 values -- the only difference between you will see compared to vpa(solve(D==0,K)) is that the order is different.

[Himanshu SINGLA]

Thank you so much sir for your precious time.

I want to further plot graph with omega and V.

box on
figure(1)
 V=OMG./real(RRd);
S1=[S1 abs(V)];
plot(OMG,S1(1,:),'r','linewidth',2)

when I use this command graph is showing blank.

[Walter Roberson]

You have

 OMG2=Q(2:2:10);
%  OMG=8.2:.5:20;
  for  OMG1=1:size(OMG2,2) 
       OMG=OMG2(1,OMG1);
       OMGG=[OMGG;OMG];
  end

This leaves you with OMGG being a column vector holding the same thing as the values in OMG2, which is code that would be better written as

  OMGG = OMG2(:);

It also leaves you with OMG being a scalar that is the last value in OMG2 .

Then later you ask to plot() using that scalar as your x coordinate.

Meanwhile, RRd will be a column vector, so V will be a column vector. S1 is empty, so [S1 abs(V)] is the same as abs(V) and that will be a column vector because V is a column vector. Then you take the first row of S1, and as it is a column vector the first row is just the first element. You then ask to plot that as the y coordinate in plot.

So you are asking to plot() with a scalar x coordinate and a scalar y coordinate. You do not specify any linestyle or marker information, so those will default to linestyle '-' and no marker (unless you have drawn a lot of things already in the same axes, in which case a different linestyle might be used.) You have a scalar point, and requested no markers and requested a line be drawn. But lines can only be drawn when there are two more more adjacent finite points in the same input array.

If you were to change the 'r' to 'r*' then you would get a single point plotted somewhere.