Roots of five degree equation with variable coefficients

syms x y P
a  = 3.70575*10^7;    %6.448*10^10/1740;
b = 1 +(x^2/y^2);     %(*we have two variables omega \ and Ohm*)
c = 3.67471*10^6;       %.6394*10^10/1740;
d = 3.67471*10^6;      %.6394*10^10/1740;
e = 6.54253*10^6 ;    %1.1384*10^10/1740;
f = 2* 1i* x/y;
g = 9.41954*10^6;   %1.639*10^10/1740;
h = 2.6028*10^12;      %0.6394*300*10^10/1740*2361;
I = 66.7593 *(0.05 +1i/y);    %49.2/(1740)*(2361)*(.05 +1i/x);
J = 20.3534;        % .05*300/1740*2361;
k = 146.05 ;        % 2*10^6*300/1740*2361;
l = 3.197*10^8/(0.05 +1i/x) ;   % .05*0.6394*10^10/(.05 + 1i/x);
m = 0.025/(.05 + 1i/x);
n = 0.025/(.05 + 1i/x);
q = 0.005/(.05 + 1i/x);
 r = 0.0000277624;            %1.1384*10^10/1740*1.753*10^(-15)*(49.2)^2;
s = 0.00634751   ;          %h/1740*1.753*10^(-15)*(49.2)^2;
 t = 2.43873*10^(-17);      %0.01/1740*1.753*10^(-15)*(49.2)^2;
u=1.9977*10^18;
A=[ -a+b*P -c -d e f*P; -f*P 0 0 0 -g+b*P; h*P I-P -J*P k*P 0; l -m n-P q 0; r -s -t -u+P 0];     %Matrix of order 5 
 D=det(A);%determinant of above matrix 
    %writing determinant as a polynomial in P(five degree polynomial, coefficients are as follow)

      b0= (- b^2-f^2);
     b1=(a*b + b*g + I*b^2 + I*f^2 + b^2*n + f^2*n + b^2*u + f^2*u + f^2*k*s + b*c*h + J*b^2*m + J*f^2*m + b^2*k*s);
     b2=(- a*g - b^2*n*u + b^2*q*t - f^2*n*u + f^2*q*t - I*a*b - I*b*g - a*b*n + b*d*l - c*g*h - b*g*n - a*b*u + b*e*r - b*g*u - I*b^2*n - I*f^2*n - I*b^2*u - I*f^2*u - a*b*k*s - b*c*h*u - b*c*k*r - b*e*h*s - b*g*k*s - J*b^2*m*u - J*b^2*q*s - J*f^2*m*u - J*f^2*q*s - b^2*k*m*t - b^2*k*n*s - f^2*k*m*t - f^2*k*n*s - J*a*b*m - J*b*c*l - J*b*g*m - b*c*h*n - b*d*h*m);
     b3= (I*a*g + a*g*n - d*g*l + a*g*u - e*g*r + a*b*n*u + a*g*k*s - b*d*l*u - b*e*l*t - b*e*n*r + c*g*h*u + c*g*k*r + e*g*h*s - a*b*q*t - b*d*q*r + b*g*n*u - b*g*q*t + I*b^2*n*u - I*b^2*q*t + I*f^2*n*u - I*f^2*q*t + I*a*b*n - I*b*d*l + J*a*g*m + J*c*g*l + I*b*g*n + I*a*b*u - I*b*e*r + I*b*g*u + c*g*h*n + d*g*h*m + J*a*b*m*u + J*b*c*l*u + J*b*e*l*s - J*b*e*m*r + J*a*b*q*s + J*b*c*q*r + J*b*g*m*u + J*b*g*q*s + a*b*k*m*t + a*b*k*n*s + b*c*h*n*u + b*c*k*l*t + b*c*k*n*r + b*d*h*m*u - b*d*k*l*s + b*d*k*m*r + b*e*h*m*t + b*e*h*n*s - b*c*h*q*t + b*d*h*q*s + b*g*k*m*t + b*g*k*n*s);
     b4= (- a*g*n*u + d*g*l*u + e*g*l*t + e*g*n*r + a*g*q*t + d*g*q*r - I*a*g*n + I*d*g*l - I*a*g*u + I*e*g*r - I*a*b*n*u + I*b*d*l*u + I*b*e*l*t + I*b*e*n*r + I*a*b*q*t + I*b*d*q*r - J*a*g*m*u - J*c*g*l*u - J*e*g*l*s + J*e*g*m*r - I*b*g*n*u - J*a*g*q*s - J*c*g*q*r + I*b*g*q*t - a*g*k*m*t - a*g*k*n*s - c*g*h*n*u - c*g*k*l*t - c*g*k*n*r - d*g*h*m*u + d*g*k*l*s - d*g*k*m*r - e*g*h*m*t - e*g*h*n*s + c*g*h*q*t - d*g*h*q*s);
     b5=(I*a*g*n*u - I*d*g*l*u - I*e*g*l*t - I*e*g*n*r - I*a*g*q*t - I*d*g*q*r);
  Y=[b0 b1 b2 b3 b4 b5];
      R  =roots(Y)

I am getting roots in terms of new variable 'z', I need roots in terms of x and y only. As i have to plot roots(P1,P2,P3,P4,P5) vs x vs y

Answer 1

[Torsten]

In general, there is no analytical formula for the roots of polynomials of degree > 4. You will have to assign values to x and y and solve for each combination of x and y separately using the "roots" command.

Furthermore, I doubt that "roots" works with symbolic expressions anyway.

Best wishes

Torsten.

Answer 2

[Walter Roberson]

Do not worry about it. Those z* variables you see are "bound" variables, dummy variables used to express the root placeholder. MATLAB will take care of their value.

However, computation of the roots will be slow.

Plotting will be a bit tricky. With the constants that you provide, when you start to evaluate to a high enough precision, all 5 of the roots are complex. If your precision is lower, such as the default 30, then you could be fooled into thinking you had a valid real-valued root or perhaps even think you had two.

[Himanshu SINGLA]

I need roots in terms of x and y, So that i can plot graph of ''real part of square root of 'root 1' vs x vs y'' in three dimension. But matlab says z is undefined variable.

[Walter Roberson]

R1 = R(1);
[X, Y] = ndgrid(linspace(-1, 1.01, 33));
Z = real(double(subs(sqrt(R1), {x, y}, {X, Y})));
surf(X, Y, Z)

Two notes about this:

1) You need to avoid some combinations of values. In particular, X = -1, Y = +1 results in the ^5 coefficient canceling, which leads to

Error using symengine
First argument must not be 0.

This is why I used 1.01 in the upper bound, so that the Y = 1 exactly is not generated so the class is side-stepped.

2) For reasons I am not clear on at the moment, if you put the real() on the symbolic expression instead of the double expression, the result can end up complex. That is a bug, which I will investigate further and report to MATLAB.

[Walter Roberson]

Yeah, using real() and imag() on a symbolic root() expression can lead to completely wrong results in R2018a, and I have just reported that. In different sessions with the same formula, real() returned a complex value and imag() returned 0 in one of the sessions, and in the other session, real() returned 0 and imag() returned 1i times the complex value.

[Himanshu SINGLA]

I tried above but my program did not run. So i tried as follow by taking bit simple values.

 x=linspace(1,100,50);
 for y=linspace(1,10,50);

x1=x./y;
x2=x.^2;
x3=y.^2;
x4=x2./x3; 

a=3.70575;
b = 1+x4;
c = 3.67471;
d = 3.67471;
e = 6.54253;
f = (2.i.*x1);  
h = 4.66927;
I=1.19762./(0.05 + 1i./y);
J=3.65129;
k=.0001460;
l = 3.197/(0.05 +1i./y);
m=  2.5/(0.05 + 1i./y);
n= 2.5/(0.05 + 1i./y);
q= 0.5/(0.05 + 1i./y);
r= 1.54182;
s=2.70875;
t= 1.31138;
u=  1.20929;

   b0= (- b^2-f^2);
   b1=(a*b + b*g + I*b^2 + I*f^2 + b^2*n + f^2*n + b^2*u + f^2*u + f^2*k*s + b*c*h + J*b^2*m + J*f^2*m + b^2*k*s);
   b2=(- a*g - b^2*n*u + b^2*q*t - f^2*n*u + f^2*q*t - I*a*b - I*b*g - a*b*n + b*d*l - c*g*h - b*g*n - a*b*u + b*e*r - b*g*u - I*b^2*n - I*f^2*n - I*b^2*u - I*f^2*u - a*b*k*s - b*c*h*u - b*c*k*r - b*e*h*s - b*g*k*s - J*b^2*m*u - J*b^2*q*s - J*f^2*m*u - J*f^2*q*s - b^2*k*m*t - b^2*k*n*s - f^2*k*m*t - f^2*k*n*s - J*a*b*m - J*b*c*l - J*b*g*m - b*c*h*n - b*d*h*m);
   b3= (I*a*g + a*g*n - d*g*l + a*g*u - e*g*r + a*b*n*u + a*g*k*s - b*d*l*u - b*e*l*t - b*e*n*r + c*g*h*u + c*g*k*r + e*g*h*s - a*b*q*t - b*d*q*r + b*g*n*u - b*g*q*t + I*b^2*n*u - I*b^2*q*t + I*f^2*n*u - I*f^2*q*t + I*a*b*n - I*b*d*l + J*a*g*m + J*c*g*l + I*b*g*n + I*a*b*u - I*b*e*r + I*b*g*u + c*g*h*n + d*g*h*m + J*a*b*m*u + J*b*c*l*u + J*b*e*l*s - J*b*e*m*r + J*a*b*q*s + J*b*c*q*r + J*b*g*m*u + J*b*g*q*s + a*b*k*m*t + a*b*k*n*s + b*c*h*n*u + b*c*k*l*t + b*c*k*n*r + b*d*h*m*u - b*d*k*l*s + b*d*k*m*r + b*e*h*m*t + b*e*h*n*s - b*c*h*q*t + b*d*h*q*s + b*g*k*m*t + b*g*k*n*s);
   b4= (- a*g*n*u + d*g*l*u + e*g*l*t + e*g*n*r + a*g*q*t + d*g*q*r - I*a*g*n + I*d*g*l - I*a*g*u + I*e*g*r - I*a*b*n*u + I*b*d*l*u + I*b*e*l*t + I*b*e*n*r + I*a*b*q*t + I*b*d*q*r - J*a*g*m*u - J*c*g*l*u - J*e*g*l*s + J*e*g*m*r - I*b*g*n*u - J*a*g*q*s - J*c*g*q*r + I*b*g*q*t - a*g*k*m*t - a*g*k*n*s - c*g*h*n*u - c*g*k*l*t - c*g*k*n*r - d*g*h*m*u + d*g*k*l*s - d*g*k*m*r - e*g*h*m*t - e*g*h*n*s + c*g*h*q*t - d*g*h*q*s);
   b5=(I*a*g*n*u - I*d*g*l*u - I*e*g*l*t - I*e*g*n*r - I*a*g*q*t - I*d*g*q*r);
   Y=[b0 b1 b2 b3 b4 b5];
    R=roots(Y);
    R1=R(1);
    R2=R(2);
    R3=R(3);
    R4=R(4);
    %R5=R(5);
    R6 = vpa(R2);
r2=sqrt(R6);
r3=real(r2);
plot3(r3,x1,y,'o')
    hold on
  end

If i remove comment before the R5, the program show error for R5. Why does it happen?

its graph is full of dots, I want that as curve or surface.

[Walter Roberson]

You were missing a definition for g; I copied it from what you had previously posted.

Because of your x=linspace(1,100,50) a number of your variables are numeric vectors. A number of places in your code had to be changed from ^ to .^ and probably a number had to be changed from * to .* .

Because of the vector x, your b0, b1, b2, etc are all numeric vectors of length 50, except that b5 happens to come out as a scalar. When you then do

      Y=[b0 b1 b2 b3 b4 b5];

you are constructing a numeric vector of length 251. roots() of that is then the roots of a polynomial of degree 250, giving you 250 answers.

You need to reconsider your use of x.

[Himanshu SINGLA]

i will retry it more carefully by taking care of all things that you mentioned, and then will tell you if the problem get solved.