MATLAB: Solution of irradical equation in matlab

determinantequationMATLABmatrix

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??

Best Answer

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

●
Related Question