Find all roots of non linear equation

Hi,

i want to find the all roots of an non linear equation without using the symbolic toolbox because of runtime issues. So i tried to use fzero and i found this helpful thread: https://de.mathworks.com/matlabcentral/answers/160117-how-to-use-fzero-and-choose-correctly-initial-guess

So i added my equation to it. The problem is that it doesn't fin initial guesses because of the step size. if i am using x = linspace(-1,1); it finds the root but not if i have a bigger intervall. And the problem is that normally i don't know the intervall that exactly.

1. Has sb an idea how i could improve the code to make it work for me ? Or has sb any idea how i could solve this problem with another method ? I cannot use roots() or sth else for polynoms because the equation is no polynom.

if true
  b1=0
b2=0
v1=0
v2=0
s1=0
s2=0.175
j2=0.28
tau =300
x = linspace(-10,10);               % Interval To Evaluate Over
f = @(x)-s2 + s1 + (v1 + (b1 - j2.*x).^2/(2.*j2) + (j2.*x.^2)/2 - (b1.*(b1 - j2.*x))/j2).*(tau - x - (2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(2.*j2) + (b1 - j2.*x)/j2 + (2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(2.*j2)) + (2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2)).^3/(48.*j2.^2) - (b1 - j2.*x).^3/(6.*j2.^2) + (j2.*x.^3)/3 + x.*(v1 + (b1 - j2.*x).^2/(2.*j2) - (b1.*(b1 - j2.*x))/j2) - (v1.*(b1 - j2.*x))/j2 + (2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(3/2))/(24.*j2.^2) + (b1.*(b1 - j2.*x).^2)/(2.*j2.^2) + ((2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2)).*(v1 - (j2.*v1 - j2.*v2 - b1.^2/2 + b2.^2/2 + j2.^2.*x.^2)/(2.*j2) + (b1 - j2.*x).^2/(2.*j2) + (j2.*x.^2)/2 - (b1.*(b1 - j2.*x))/j2))/(2.*j2) + (2.^(1/2).*(2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2)).^2.*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(16.*j2.^2) - (2.^(1/2).*(v1 + (b1 - j2.*x).^2/(2.*j2) + (j2.*x.^2)/2 - (b1.*(b1 - j2.*x))/j2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(2.*j2);
fplot(f)
fx = f(x);                          % Function Evaluated Over ‘x’
cs = fx.*circshift(fx,-1,2);        % Product Negative At Zero-Crossings
xc = x(cs <= 0);                    % Values Of ‘x’ Near Zero Crossings
for k1 = 1:length(xc)
  fz(k1) = fzero(f, xc(k1));      % Use ‘xc’ As Initial Zero Estimate
end
syms x
f = -s2 + s1 + (v1 + (b1 - j2.*x).^2/(2.*j2) + (j2.*x.^2)/2 - (b1.*(b1 - j2.*x))/j2).*(tau - x - (2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(2.*j2) + (b1 - j2.*x)/j2 + (2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(2.*j2)) + (2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2)).^3/(48.*j2.^2) - (b1 - j2.*x).^3/(6.*j2.^2) + (j2.*x.^3)/3 + x.*(v1 + (b1 - j2.*x).^2/(2.*j2) - (b1.*(b1 - j2.*x))/j2) - (v1.*(b1 - j2.*x))/j2 + (2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(3/2))/(24.*j2.^2) + (b1.*(b1 - j2.*x).^2)/(2.*j2.^2) + ((2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2)).*(v1 - (j2.*v1 - j2.*v2 - b1.^2/2 + b2.^2/2 + j2.^2.*x.^2)/(2.*j2) + (b1 - j2.*x).^2/(2.*j2) + (j2.*x.^2)/2 - (b1.*(b1 - j2.*x))/j2))/(2.*j2) + (2.^(1/2).*(2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2)).^2.*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(16.*j2.^2) - (2.^(1/2).*(v1 + (b1 - j2.*x).^2/(2.*j2) + (j2.*x.^2)/2 - (b1.*(b1 - j2.*x))/j2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(2.*j2);
vpa(solve(f,x))
end

[Dimitris Kalogiros]

Lets clarify something.

Is this the equation you want to solve ?

[IlPadrino]

Yes that's the equation... And the problem is that my x is in the square root. This means that i can't build the coefficients for the polynom. That's why i am searching an solution without polynoms

[Walter Roberson]

If you use syms tau instead of a numeric value, and you

simplify(solve(f, x))

then you get back two specific roots,

 -(10^(1/2)*(1/tau)^(1/2))/4
  (10^(1/2)*(1/tau)^(1/2))/4

You can work with that formula to establish the code, and then if calculation of the value for any given tau is proving too slow, or you need to be able to compile, then you can use the expression (or a variation of it.)

[IlPadrino]

Hi thanks for the answer :) but i cannot use solve if the parameters are inserted.

[Walter Roberson]

Which are the parameters?

[IlPadrino]

b1 b2 v1 v2 s1 s2 j2 tau. The parameters which are in the equation

Answer 1

[Star Strider]

That algorithm will work when the increment in your ‘x’ vector is small enough to detect the sign change.

The solution is to simply add more points to increase the resolution of the vector:

x = linspace(-10,10,1000);          % Interval To Evaluate Over

that then produces:

fz =
  -0.045636604654443   0.045643546458764

as desired, and acceptably close to the symbolic result.

[IlPadrino]

Thanks star strider for the help. But my problem is that it could be that the root is much smaller if for example tau is much bigger. Because normally the interval in which i want to search ist 0 to tau.this means that the resolution hast to be much higher... But i can't set the resolution that high. Also it can happen that i don't get the root if the resolution is to low. Isn't another way to solve this? Thanks in advance. I really appreciate your help.

[Star Strider]

After experimenting a bit, I came up with this solution (that I successfully tested with ‘tau’ up to 3E+017 and as low as 3E-17):

[xmin,fxmin] = fminsearch(f, 10)                            % Find Minimum Of Quadratic Equation
for k1 = 1:2
    [xv(k1),fxval(k1)] = fzero(f, (fxmin*(-1)^k1)/tau);     % Find Both Roots, Scale Initial Estimate By ‘1/tau’
end
xroots = xv                                                 % Display Roots (Optional)
function_values = fxval                                     % Display Roots Function Values (Optional)

The solution to having it work with the fzero call is to scale the initial estimate of the fzero result by the inverse of ‘tau’.

[IlPadrino]

Thank you very much :) That helped me a lot but didn't solved my problem :/ because the parameters can vary: For example: b1=0.3 b2=0.5 v1=0.2 v2=0.5 s1=0.001 s2=0.175 j2=0.28 tau =20 Then i have the problem that it doesn't find the zeros :/

[Star Strider]

My pleasure. With your new parameter set, replacing the fzero function with the Optimization Toolbox fsolve function produces the roots:

for k1 = 1:2
    [xv(k1,:),fxval(k1,:)] = fsolve(f, (-1*(-1)^k1)/tau);   % Find Both Roots, Scale Initial Estimate By ‘1/tau’
end
xroots = xv                                                 % Display Roots (Optional)
function_values = fxval                                     % Display Roots Function Values (Optional)

Note that they are complex with these parameters:

xroots =
 -0.010218893299935 + 0.395164734474393i
 -0.010218910315512 - 0.395164738332085i

With the value of the function at those roots also being complex. You need to decide how you want to deal with that in your code.

[IlPadrino]

But how do you decide how many roots he has to search? Because right now you set k1 to 2. But this cozld be changing.

[Star Strider]

It could change. To search for more roots, increase the ‘k1’ limit, and provide new initial estimates:

for k1 = 1:4
    [xv(k1,:),fxval(k1,:)] = fsolve(f, (-1/k1)*((-1)^k1)/tau);  % Find Both Roots, Scale Initial Estimate By ‘1/tau’
end

xroots =
 -0.010218893299935 + 0.395164734474393i
  0.002137898025294 - 0.395703760614642i
 -0.010218911868143 + 0.395164745783536i
  0.002137898022517 - 0.395703760622581i

These (using your most recent set of constants) are acceptably close to the symbolic solve results.

The initial estimates are calculated as ‘(-1/k1)*((-1)^k1)/tau’.

Answer 2

[Dimitris Kalogiros]

Dear IlPadrino

Your equation is not a "nonlinear" equation. If you try to search its roots fist at the interval [0, +inf) and then at (-inf, 0) , you can rid of this annoying square root.

Have a look at the following piece of code :

clear all;
clc;

% Initial problem
b1=0;
b2=0;
v1=0;
v2=0;
s1=0;
s2=0.175;
j2=0.28;
tau =300;

syms x

% initial form of f
f = -s2 + s1 + (v1 + (b1 - j2.*x).^2/(2.*j2) + (j2.*x.^2)/2 - (b1.*(b1 - j2.*x))/j2).*(tau - x - (2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(2.*j2) + (b1 - j2.*x)/j2 + (2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(2.*j2)) + (2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2)).^3/(48.*j2.^2) - (b1 - j2.*x).^3/(6.*j2.^2) + (j2.*x.^3)/3 + x.*(v1 + (b1 - j2.*x).^2/(2.*j2) - (b1.*(b1 - j2.*x))/j2) - (v1.*(b1 - j2.*x))/j2 + (2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(3/2))/(24.*j2.^2) + (b1.*(b1 - j2.*x).^2)/(2.*j2.^2) + ((2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2)).*(v1 - (j2.*v1 - j2.*v2 - b1.^2/2 + b2.^2/2 + j2.^2.*x.^2)/(2.*j2) + (b1 - j2.*x).^2/(2.*j2) + (j2.*x.^2)/2 - (b1.*(b1 - j2.*x))/j2))/(2.*j2) + (2.^(1/2).*(2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2)).^2.*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(16.*j2.^2) - (2.^(1/2).*(v1 + (b1 - j2.*x).^2/(2.*j2) + (j2.*x.^2)/2 - (b1.*(b1 - j2.*x))/j2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(2.*j2)

%simplify f

f=expand(f)

% searching for roots at interval [0 +inf)

assume(x>=0);

f1=simplify(f)

sol_Right=solve(f1==0)

%searching for roots at (-inf , 0)

assume(x<0)

f2=simplify(f)

sol_left=vpasolve(f2==0)

% Graphical demonstration
fplot(f,[-.1 .1]); grid on; hold on; fplot(0,[-.1, .1]);

I've run this , and you can find results of it, inside the file I've attached.

[James Tursa]

It certainly is nonlinear. I think what you meant to write is that with the proper interpretation it can be reduced to a polynomial.

[IlPadrino]

Hi Dimitris Kalogiros big thanks for your Help I really appreciate it. But the problem is that the parameter are inserted at the very end of my process. So that means i need general expression of the polynom where i just insert the parameters and then calculate my x variable with an numerical approximation (without symbolic Toolbox) Tha's why i can't do it that way. I tried it with the symbolic expression of my non linear equation to simplyfy and then expand it to generate a polynom but without success. In Conclusion if i want to create an polynom i have to create it with the parameters in it.

Answer 3

[Walter Roberson]

The solution involves the roots of a quartic.

P4 = -9*sqrt(2)*(-2*j2^2*tau^2+((-4*b1+4*b2)*tau-6*v1+6*v2)*j2+(b1+5*b2)*(b1-b2));
P3 = (72*tau*v2+72*s1-72*s2)*j2^2+((-72*v1+72*v2)*b1-36*b2^2*tau)*j2+24*b1^3-36*b1*b2^2+12*b2^3;
P2 = -(3*(-24*tau*(tau*v2+s1-s2)*j2^3+(((24*v1-48*v2)*tau-24*s1+24*s2)*b1+12*b2^2*tau^2+(24*tau*v2+24*s1-24*s2)*b2+36*(v1-v2)^2)*j2^2-(8*(tau*b1^2+(b2*tau+3*v1*(1/2)-3*v2*(1/2))*b1-(2*(b2*tau-9*v1*(1/4)+9*v2*(1/4)))*b2))*(b1-b2)*j2+(b1^2+10*b1*b2+13*b2^2)*(b1-b2)^2))*sqrt(2);
P1 = 0;
P0 = -sqrt(2)*(-72*(tau*v2+s1-s2)^2*j2^4+((144*(v1-v2))*(tau*v2+s1-s2)*b1+72*tau*(tau*v2+s1-s2)*b2^2-72*(v1-v2)^3)*j2^3+((-48*tau*v2-48*s1+48*s2)*b1^3+36*(v1-v2)^2*b1^2+72*b2^2*(-tau*v1+2*tau*v2+s1-s2)*b1-(18*(b2^2*tau^2+(4*v2*tau*(1/3)+4*s1*(1/3)-4*s2*(1/3))*b2+6*(v1-v2)^2))*b2^2)*j2^2-(6*(b1-b2))*((v1-v2)*b1^3-4*b2*(b2*tau-(1/4)*v1+(1/4)*v2)*b1^2-(4*(b2*tau+5*v1*(1/4)-5*v2*(1/4)))*b2^2*b1+(2*(b2*tau-9*v1*(1/2)+9*v2*(1/2)))*b2^3)*j2+(b1^4+2*b1^3*b2-10*b1*b2^3-11*b2^4)*(b1-b2)^2);
RRR = roots([P4, P3, P2, P1, P0]);

xroots = (3*RRR.^3*sqrt(2) + (12*tau*v2+12*s1-12*s2)*j2^2 + ((-12*v1+12*v2)*b1 - 6*b2^2*tau + 6*RRR.^2*tau)*j2 + 4*b1^3 + (6*RRR.^2-6*b2^2)*b1 + 2*b2^3 - 6*RRR.^2*b2) ./ (6*j2*((-2*v1+2*v2)*j2 + b1^2 - b2^2 + RRR.^2));

These are not guaranteed to be real-valued, and are not guaranteed to be unique.

[IlPadrino]

Wow, How did you do this ? Because i have some other Equation which i have also to transform.

[Walter Roberson]

I dropped it into Maple and used solve(). That gave me an expression in terms of the Maple RootOf() operator, which is the general form for roots of equations. The inside of the RootOf() operator was a quartic:

RootOf(24*sqrt(2)*b2^3*j2^2*s1+(36*j2*tau+36*b1-36*b2)*_Z^5+(-36*b2^2*j2*tau+72*j2^2*tau*v2+24*b1^3-36*b1*b2^2-72*b1*j2*v1+72*b1*j2*v2+12*b2^3+72*j2^2*s1-72*j2^2*s2)*_Z^3-24*sqrt(2)*b2^3*j2^2*s2-54*sqrt(2)*b2^4*j2*v2-36*sqrt(2)*b1^2*j2^2*v2^2+108*sqrt(2)*b2^2*j2^2*v1^2+54*sqrt(2)*b2^4*j2*v1+72*sqrt(2)*b1*b2^2*j2^2*tau*v1-144*sqrt(2)*b1*b2^2*j2^2*tau*v2-144*sqrt(2)*b1*j2^3*tau*v1*v2-72*sqrt(2)*j2^3*v2^3+72*sqrt(2)*j2^4*s1^2+3*sqrt(2)*b1^4*b2^2-9*sqrt(2)*b1^2*b2^4+8*sqrt(2)*b1^3*b2^3-12*sqrt(2)*b1*b2^5+72*sqrt(2)*j2^4*s2^2+72*sqrt(2)*j2^3*v1^3+144*sqrt(2)*j2^4*s1*tau*v2+72*sqrt(2)*b2^2*j2^3*s2*tau-72*sqrt(2)*b2^2*j2^3*s1*tau+144*sqrt(2)*b1*j2^3*s2*v1-36*sqrt(2)*b1^2*j2^2*v1^2-48*sqrt(2)*b1^3*j2^2*s2+11*sqrt(2)*b2^6-sqrt(2)*b1^6+(36*sqrt(2)*b1^2*j2*v1-36*sqrt(2)*b1^2*j2*v2+72*sqrt(2)*b1*j2^2*s1-72*sqrt(2)*b1*j2^2*s2-108*sqrt(2)*b2^2*j2*v1+108*sqrt(2)*b2^2*j2*v2-72*sqrt(2)*b2*j2^2*s1+72*sqrt(2)*b2*j2^2*s2+216*sqrt(2)*j2^2*v1*v2-36*sqrt(2)*b2^2*j2^2*tau^2+72*sqrt(2)*j2^3*tau^2*v2+24*sqrt(2)*b1^3*j2*tau+48*sqrt(2)*b2^3*j2*tau+72*sqrt(2)*j2^3*s1*tau-72*sqrt(2)*j2^3*s2*tau-3*sqrt(2)*b1^4-39*sqrt(2)*b2^4-24*sqrt(2)*b1^3*b2-108*sqrt(2)*j2^2*v2^2+48*sqrt(2)*b1*b2^3+18*sqrt(2)*b1^2*b2^2-108*sqrt(2)*j2^2*v1^2-72*sqrt(2)*b1*b2^2*j2*tau-72*sqrt(2)*b1*j2^2*tau*v1+144*sqrt(2)*b1*j2^2*tau*v2-72*sqrt(2)*b2*j2^2*tau*v2+72*sqrt(2)*b1*b2*j2*v1-72*sqrt(2)*b1*b2*j2*v2)*_Z^2+(18*sqrt(2)*j2^2*tau^2+36*sqrt(2)*b1*j2*tau-36*sqrt(2)*b2*j2*tau-9*sqrt(2)*b1^2-36*sqrt(2)*b1*b2+45*sqrt(2)*b2^2+54*sqrt(2)*j2*v1-54*sqrt(2)*j2*v2)*_Z^4+144*sqrt(2)*b1*j2^3*s1*v2-144*sqrt(2)*b1*j2^3*s1*v1+72*sqrt(2)*b1*b2^2*j2^2*s2-72*sqrt(2)*b1*b2^2*j2^2*s1+24*sqrt(2)*b1*b2^3*j2*v2-24*sqrt(2)*b1*b2^3*j2*v1+48*sqrt(2)*b1^3*j2^2*s1-6*sqrt(2)*b1^4*j2*v2+6*sqrt(2)*b1^4*j2*v1-12*sqrt(2)*b2^5*j2*tau+72*sqrt(2)*j2^4*tau^2*v2^2+18*sqrt(2)*b2^4*j2^2*tau^2+216*sqrt(2)*j2^3*v1*v2^2-216*sqrt(2)*j2^3*v1^2*v2-144*sqrt(2)*j2^4*s1*s2+108*sqrt(2)*b2^2*j2^2*v2^2+72*sqrt(2)*b1^2*j2^2*v1*v2+36*sqrt(2)*b1^2*b2^2*j2*v2-36*sqrt(2)*b1^2*b2^2*j2*v1-216*sqrt(2)*b2^2*j2^2*v1*v2-144*sqrt(2)*b1*j2^3*s2*v2+48*sqrt(2)*b1^3*j2^2*tau*v2-24*sqrt(2)*b1^3*b2^2*j2*tau-72*sqrt(2)*b2^2*j2^3*tau^2*v2+24*sqrt(2)*b2^3*j2^2*tau*v2+144*sqrt(2)*b1*j2^3*tau*v2^2+36*sqrt(2)*b1*b2^4*j2*tau-144*sqrt(2)*j2^4*s2*tau*v2)

... Ah, I missed the _Z^5 term! It is a quintic, not a quartic! So correction, add

P5 = 36*j2*tau+36*b1-36*b2;

and then

RRR = roots([P5, P4, P3, P2, P1, P0]);

Anyhow, used Maple's coeff() to extract the coefficients for each polynomial term and assigned them into variables for clarity.

[IlPadrino]

ok, cool thx ! Is it maybe possible to append your Maple worksheet ? i never worked with Maple so it would be cool to have one working sheet :) Because i tried it now with my second equation and i get an error from Maple :invalid arguments to coeffs

    s2 := s1+th3*(v1+(-j2*x+b1)^2/(2*j2)+(1/2)*j2*x^2-b1*(-j2*x+b1)/j2+j2*x*th1)+((1/2)*j2*th3+(1/2)*(4*j2^2*th1*x+j2^2*th3^2+4*j2^2*x^2-2*b1^2+2*b2^2+4*j2*v1-4*j2*v2)^(1/2))^3/(6*j2^2)+x*(v1+(-j2*x+b1)^2/(2*j2)-b1*(-j2*x+b1)/j2+j2*x*th1)-(j2*th3-2*b2+(4*j2^2*th1*x+j2^2*th3^2+4*j2^2*x^2-2*b1^2+2*b2^2+4*j2*v1-4*j2*v2)^(1/2))^3/(48*j2^2)-(-j2*x+b1)^3/(6*j2^2)+(1/3)*j2*x^3+th1*(v1+(-j2*x+b1)^2/(2*j2)-b1*(-j2*x+b1)/j2)+(1/2)*th3^2*((1/2)*j2*th3+(1/2)*(4*j2^2*th1*x+j2^2*th3^2+4*j2^2*x^2-2*b1^2+2*b2^2+4*j2*v1-4*j2*v2)^(1/2))+(v1+(-j2*x+b1)^2/(2*j2)+(1/2)*j2*x^2-b1*(-j2*x+b1)/j2+j2*x*th1)*(tau-x-th1-th3+(j2*th3-2*b2+(4*j2^2*th1*x+j2^2*th3^2+4*j2^2*x^2-2*b1^2+2*b2^2+4*j2*v1-4*j2*v2)^(1/2))/(2*j2)+(-j2*x+b1)/j2+((1/2)*j2*th3+(1/2)*(4*j2^2*th1*x+j2^2*th3^2+4*j2^2*x^2-2*b1^2+2*b2^2+4*j2*v1-4*j2*v2)^(1/2))/j2)+(1/2)*j2*x*th1^2-v1*(-j2*x+b1)/j2-(j2*th3-2*b2+(4*j2^2*th1*x+j2^2*th3^2+4*j2^2*x^2-2*b1^2+2*b2^2+4*j2*v1-4*j2*v2)^(1/2))*(v1-((1/2)*j2*th3+(1/2)*(4*j2^2*th1*x+j2^2*th3^2+4*j2^2*x^2-2*b1^2+2*b2^2+4*j2*v1-4*j2*v2)^(1/2))^2/(2*j2)+th3*((1/2)*j2*th3+(1/2)*(4*j2^2*th1*x+j2^2*th3^2+4*j2^2*x^2-2*b1^2+2*b2^2+4*j2*v1-4*j2*v2)^(1/2))+(-j2*x+b1)^2/(2*j2)+(1/2)*j2*x^2-b1*(-j2*x+b1)/j2+j2*x*th1)/(2*j2)-((1/2)*j2*th3+(1/2)*(4*j2^2*th1*x+j2^2*th3^2+4*j2^2*x^2-2*b1^2+2*b2^2+4*j2*v1-4*j2*v2)^(1/2))*(v1+(-j2*x+b1)^2/(2*j2)+(1/2)*j2*x^2-b1*(-j2*x+b1)/j2+j2*x*th1)/j2+b1*(-j2*x+b1)^2/(2*j2^2)+((1/2)*j2*th3+(1/2)*(4*j2^2*th1*x+j2^2*th3^2+4*j2^2*x^2-2*b1^2+2*b2^2+4*j2*v1-4*j2*v2)^(1/2))*(j2*th3-2*b2+(4*j2^2*th1*x+j2^2*th3^2+4*j2^2*x^2-2*b1^2+2*b2^2+4*j2*v1-4*j2*v2)^(1/2))^2/(8*j2^2)

[Walter Roberson]

EQN := s1-s2+(2*b2-2^(1/2)*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(1/2))^3/(48*j2^2)-(-j2*x+b1)^3/(6*j2^2)+(v1+(-j2*x+b1)^2/(2*j2)+(1/2)*j2*x^2-b1*(-j2*x+b1)/j2)*(tau-x-(2*b2-2^(1/2)*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(1/2))/(2*j2)+(-j2*x+b1)/j2+2^(1/2)*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(1/2)/(2*j2))+(1/3)*j2*x^3+x*(v1+(-j2*x+b1)^2/(2*j2)-b1*(-j2*x+b1)/j2)+(2*b2-2^(1/2)*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(1/2))*(v1-(j2^2*x^2+j2*v1-j2*v2-(1/2)*b1^2+(1/2)*b2^2)/(2*j2)+(-j2*x+b1)^2/(2*j2)+(1/2)*j2*x^2-b1*(-j2*x+b1)/j2)/(2*j2)-v1*(-j2*x+b1)/j2+2^(1/2)*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(3/2)/(24*j2^2)+b1*(-j2*x+b1)^2/(2*j2^2)+2^(1/2)*(2*b2-2^(1/2)*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(1/2))^2*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(1/2)/(16*j2^2)-2^(1/2)*(v1+(-j2*x+b1)^2/(2*j2)+(1/2)*j2*x^2-b1*(-j2*x+b1)/j2)*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(1/2)/(2*j2);
Rx := solve(EQN,x);
Rt := indets(rx, RootOf);
Rexp := op(1, Rt[1]);

P5 := simplify(coeff(Rexp, _Z, 5)) assuming real;
P4 := simplify(coeff(Rexp, _Z, 4)) assuming real;
P3 := simplify(coeff(Rexp, _Z, 3)) assuming real;
P2 := simplify(coeff(Rexp, _Z, 2)) assuming real;
P1 := simplify(coeff(Rexp, _Z, 1)) assuming real;
P0 := simplify(coeff(Rexp, _Z, 0)) assuming real;
xroots := simplify( subs(Rt[1]=RRR, Rx)) assuming real;

after which in MATLAB you copy P5 through P0, put them in a vector and roots() them and assign to RRR, and use the vectorized form of xroots

[IlPadrino]

I will try it ! Big Thanks!

[Walter Roberson]

Note that it takes a while to calculate Rx.

[IlPadrino]

Okay Thx :) So i calculated the coeffs like you did and get this :

if true
  EQN := s1-s2+(2*b2-2^(1/2)*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(1/2))^3/(48*j2^2)-(-j2*x+b1)^3/(6*j2^2)+(v1+(-j2*x+b1)^2/(2*j2)+(1/2)*j2*x^2-b1*(-j2*x+b1)/j2)*(tau-x-(2*b2-2^(1/2)*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(1/2))/(2*j2)+(-j2*x+b1)/j2+2^(1/2)*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(1/2)/(2*j2))+(1/3)*j2*x^3+x*(v1+(-j2*x+b1)^2/(2*j2)-b1*(-j2*x+b1)/j2)+(2*b2-2^(1/2)*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(1/2))*(v1-(j2^2*x^2+j2*v1-j2*v2-(1/2)*b1^2+(1/2)*b2^2)/(2*j2)+(-j2*x+b1)^2/(2*j2)+(1/2)*j2*x^2-b1*(-j2*x+b1)/j2)/(2*j2)-v1*(-j2*x+b1)/j2+2^(1/2)*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(3/2)/(24*j2^2)+b1*(-j2*x+b1)^2/(2*j2^2)+2^(1/2)*(2*b2-2^(1/2)*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(1/2))^2*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(1/2)/(16*j2^2)-2^(1/2)*(v1+(-j2*x+b1)^2/(2*j2)+(1/2)*j2*x^2-b1*(-j2*x+b1)/j2)*(2*j2^2*x^2-b1^2+b2^2+2*j2*v1-2*j2*v2)^(1/2)/(2*j2)
P5 = 36*j2*tau+36*b1-36*b2
P4 = 18*sqrt(2)*j2^2*tau^2+36*sqrt(2)*b1*j2*tau-36*sqrt(2)*b2*j2*tau-9*sqrt(2)*b1^2-36*sqrt(2)*b1*b2+45*sqrt(2)*b2^2+54*sqrt(2)*j2*v1-54*sqrt(2)*j2*v2
P3 = -36*b2^2*j2*tau+72*j2^2*tau*v2+24*b1^3-36*b1*b2^2-72*b1*j2*v1+72*b1*j2*v2+12*b2^3+72*j2^2*s1-72*j2^2*s2
P2 = -36*sqrt(2)*b1^2*j2*v2-72*sqrt(2)*b1*j2^2*tau*v1+144*sqrt(2)*b1*j2^2*tau*v2-72*sqrt(2)*b2*j2^2*tau*v2-108*sqrt(2)*j2^2*v2^2+18*sqrt(2)*b1^2*b2^2+48*sqrt(2)*b1*b2^3-108*sqrt(2)*j2^2*v1^2+216*sqrt(2)*j2^2*v1*v2-36*sqrt(2)*b2^2*j2^2*tau^2+72*sqrt(2)*j2^3*tau^2*v2-24*sqrt(2)*b1^3*b2-72*sqrt(2)*j2^3*s2*tau+36*sqrt(2)*b1^2*j2*v1-72*sqrt(2)*b1*b2^2*j2*tau+48*sqrt(2)*b2^3*j2*tau+72*sqrt(2)*j2^3*s1*tau+72*sqrt(2)*b1*b2*j2*v1-72*sqrt(2)*b1*b2*j2*v2-3*sqrt(2)*b1^4-39*sqrt(2)*b2^4+24*sqrt(2)*b1^3*j2*tau+108*sqrt(2)*b2^2*j2*v2-72*sqrt(2)*b2*j2^2*s1+72*sqrt(2)*b2*j2^2*s2+72*sqrt(2)*b1*j2^2*s1-72*sqrt(2)*b1*j2^2*s2-108*sqrt(2)*b2^2*j2*v1
P1 = 0
P0 = 18*sqrt(2)*b2^4*j2^2*tau^2+72*sqrt(2)*j2^4*tau^2*v2^2+8*sqrt(2)*b1^3*b2^3-144*sqrt(2)*b1*j2^3*tau*v1*v2+36*sqrt(2)*b1^2*b2^2*j2*v2-144*sqrt(2)*j2^4*s2*tau*v2-36*sqrt(2)*b1^2*b2^2*j2*v1+216*sqrt(2)*j2^3*v1*v2^2+48*sqrt(2)*b1^3*j2^2*s1-48*sqrt(2)*b1^3*j2^2*s2-36*sqrt(2)*b1^2*j2^2*v1^2+144*sqrt(2)*j2^4*s1*tau*v2-144*sqrt(2)*b1*b2^2*j2^2*tau*v2-sqrt(2)*b1^6+11*sqrt(2)*b2^6-144*sqrt(2)*j2^4*s1*s2-216*sqrt(2)*j2^3*v1^2*v2-36*sqrt(2)*b1^2*j2^2*v2^2-12*sqrt(2)*b2^5*j2*tau+6*sqrt(2)*b1^4*j2*v1-6*sqrt(2)*b1^4*j2*v2-72*sqrt(2)*b2^2*j2^3*s1*tau+72*sqrt(2)*b2^2*j2^3*s2*tau+144*sqrt(2)*b1*j2^3*tau*v2^2+24*sqrt(2)*b2^3*j2^2*tau*v2+72*sqrt(2)*b1^2*j2^2*v1*v2-24*sqrt(2)*b1*b2^3*j2*v1+54*sqrt(2)*b2^4*j2*v1-54*sqrt(2)*b2^4*j2*v2+24*sqrt(2)*b2^3*j2^2*s1-24*sqrt(2)*b2^3*j2^2*s2-9*sqrt(2)*b1^2*b2^4+72*sqrt(2)*b1*b2^2*j2^2*tau*v1+72*sqrt(2)*b1*b2^2*j2^2*s2-72*sqrt(2)*b1*b2^2*j2^2*s1-12*sqrt(2)*b1*b2^5+3*sqrt(2)*b1^4*b2^2+72*sqrt(2)*j2^4*s1^2+72*sqrt(2)*j2^3*v1^3-72*sqrt(2)*j2^3*v2^3+72*sqrt(2)*j2^4*s2^2+24*sqrt(2)*b1*b2^3*j2*v2+108*sqrt(2)*b2^2*j2^2*v1^2+108*sqrt(2)*b2^2*j2^2*v2^2-216*sqrt(2)*b2^2*j2^2*v1*v2-72*sqrt(2)*b2^2*j2^3*tau^2*v2-24*sqrt(2)*b1^3*b2^2*j2*tau+144*sqrt(2)*b1*j2^3*s1*v2+144*sqrt(2)*b1*j2^3*s2*v1-144*sqrt(2)*b1*j2^3*s2*v2+48*sqrt(2)*b1^3*j2^2*tau*v2+36*sqrt(2)*b1*b2^4*j2*tau-144*sqrt(2)*b1*j2^3*s1*v1
RRR = roots([P5,P4, P3, P2, P1, P0])
xroots = (1/6).*(3.*sqrt(2).*RRR.^3+6.*RRR.^2.*j2.*tau-6.*tau.*j2.*b2.^2+12.*tau.*j2.^2.*v2+6.*RRR.^2.*b1-6.*RRR.^2.*b2+4.*b1.^3-6.*b1.*b2.^2-12.*b1.*j2.*v1+12.*b1.*j2.*v2+2.*b2.^3+12.*s1.*j2.^2-12.*s2.*j2.^2)/(j2.*(RRR.^2+b1.^2-b2.^2-2.*j2.*v1+2.*j2.*v2))

end

As Output i Get this :

if true
  RRR =

-3.8253 + 0.0000i
 0.0007 + 0.1804i
 0.0007 - 0.1804i
 0.0035 + 0.1802i
 0.0035 - 0.1802i

xroots =

   9.6630 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
  -0.0000 - 0.0007i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
  -0.0000 + 0.0007i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0001 + 0.0007i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0001 - 0.0007i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
  end

The problem is that i compared it to the vpa solve and the roots are completely different:

if true
  syms x
f =  -s2 + s1 + (v1 + (b1 - j2.*x).^2/(2.*j2) + (j2.*x.^2)/2 - (b1.*(b1 - j2.*x))/j2).*(tau - x - (2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(2.*j2) + (b1 - j2.*x)/j2 + (2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(2.*j2)) + (2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2)).^3/(48.*j2.^2) - (b1 - j2.*x).^3/(6.*j2.^2) + (j2.*x.^3)/3 + x.*(v1 + (b1 - j2.*x).^2/(2.*j2) - (b1.*(b1 - j2.*x))/j2) - (v1.*(b1 - j2.*x))/j2 + (2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(3/2))/(24.*j2.^2) + (b1.*(b1 - j2.*x).^2)/(2.*j2.^2) + ((2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2)).*(v1 - (j2.*v1 - j2.*v2 - b1.^2/2 + b2.^2/2 + j2.^2.*x.^2)/(2.*j2) + (b1 - j2.*x).^2/(2.*j2) + (j2.*x.^2)/2 - (b1.*(b1 - j2.*x))/j2))/(2.*j2) + (2.^(1/2).*(2.*b2 - 2.^(1/2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2)).^2.*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(16.*j2.^2) - (2.^(1/2).*(v1 + (b1 - j2.*x).^2/(2.*j2) + (j2.*x.^2)/2 - (b1.*(b1 - j2.*x))/j2).*(2.*j2.*v1 - 2.*j2.*v2 - b1.^2 + b2.^2 + 2.*j2.^2.*x.^2).^(1/2))/(2.*j2);
vpa(solve(f,x))
end

Results:

if true
 ans =

0.0021378980246453118999483001395169 + 0.39570376058876405371605295973236i
 0.002137898024645311899948300139516 - 0.39570376058876405371605295973235i
- 0.010218911817920429505659801516899 + 0.39516474594473700606935881923862i
 - 0.0102189118179204295056598015169 - 0.39516474594473700606935881923862i

end

[IlPadrino]

no idea why it doesn't work ? :/

[Walter Roberson]

Sorry, I have been busy.