Find all roots of non linear equation

Question

IlPadrino 2018-7-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/411013-find-all-roots-of-non-linear-equation

评论： Walter Roberson 2018-7-21

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

Walter Roberson 2018-7-19

Which are the parameters?

IlPadrino 2018-7-19

编辑：IlPadrino 2018-7-19

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2018-7-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/411013-find-all-roots-of-non-linear-equation#answer_329373

在 MATLAB Online 中打开

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.

6 个评论
显示 4更早的评论隐藏 4更早的评论

IlPadrino 2018-7-19

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 2018-7-19

在 MATLAB Online 中打开

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 2018-7-18

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/411013-find-all-roots-of-non-linear-equation#answer_329375

在 MATLAB Online 中打开

dkalο_equation.pdf

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.

2 个评论
显示无隐藏无

James Tursa 2018-7-18

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 2018-7-19

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 2018-7-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/411013-find-all-roots-of-non-linear-equation#answer_329490

在 MATLAB Online 中打开

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.

9 个评论
显示 7更早的评论隐藏 7更早的评论

Walter Roberson 2018-7-19

在 MATLAB Online 中打开

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 2018-7-19

编辑：Walter Roberson 2018-7-19

在 MATLAB Online 中打开

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 2018-7-19

在 MATLAB Online 中打开

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 2018-7-19

在 MATLAB Online 中打开

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 2018-7-21

no idea why it doesn't work ? :/

Walter Roberson 2018-7-21

Sorry, I have been busy.

请先登录，再进行评论。

Find all roots of non linear equation

6 个评论
显示 4更早的评论隐藏 4更早的评论

回答（3 个）

6 个评论
显示 4更早的评论隐藏 4更早的评论

2 个评论
显示无隐藏无

9 个评论
显示 7更早的评论隐藏 7更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Find all roots of non linear equation

6 个评论 显示 4更早的评论隐藏 4更早的评论

回答（3 个）

6 个评论 显示 4更早的评论隐藏 4更早的评论

2 个评论 显示 无隐藏 无

9 个评论 显示 7更早的评论隐藏 7更早的评论

另请参阅

类别

标签

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论

6 个评论
显示 4更早的评论隐藏 4更早的评论

2 个评论
显示无隐藏无

9 个评论
显示 7更早的评论隐藏 7更早的评论