Issues with solving a quartic function analytically

6 次查看（过去 30 天）
Baas28
Baas28 2018-2-13
回答： Yassine Damerdji 2022-4-20
Dear all,
I'm solving for the variable "x" that is plugged into an optimizer, after it has been solved for. Because of the nature of the optimizer, I'm not able to use the "roots"-function. Because of this I opted to solve the 4th order polynomial analytically, with the following scheme, as can be found publically on Wikipedia:
I have implemented this as follows with the relevant part or the script below:
% func = ((x)^4) + (2*(UC)*((x)^3)) + (((UC)^2) + ((UT)^2))*((x)^2) - 1;
UC = [-8.848312431038208e-17];
UT = [0.550226217684170];
a = 1;
b = 2*(UC);
c = ( ((UC).^2) + ((UT).^2) );
d = 0;
e = -1;
D_0 = ((c).^2) - (3*b*d) + (12*a*e);
D_1 = 2*((c).^3) - (9*b.*c*d) + (27*((b).^2)*e) + (27*a*((d)^2)) - (72*a*c*e);
Q = ( (D_1 + sqrt(((D_1).^2) - 4*((D_0).^3) ) ) ./ 2).^(1/3);
p = ((8*a*c) - (3*((b).^2))) ./ (8*((a).^2));
q = (((b).^3) - (4*a*b.*c) + (8*((a).^2).*d)) / (8*((a).^3));
S = (1/2) * sqrt( ((-2/3).*p) + (1/(3*a)).*(Q + (D_0)./Q) );
x_1 = (-b/(4*a)) - S + (1/2)*sqrt((-4*((S).^2)) - (2*p) + (q./S));
x_2 = (-b/(4*a)) - S - (1/2)*sqrt((-4*((S).^2)) - (2*p) + (q./S));
x_3 = (-b/(4*a)) + S + (1/2)*sqrt((-4*((S).^2)) - (2*p) - (q./S));
x_4 = (-b/(4*a)) + S - (1/2)*sqrt((-4*((S).^2)) - (2*p) - (q./S));
x = [x_1 x_2 x_3 x_4]
I've cross-checked the script for a known function, however, for some values the script produces values that do not make any sense. I've noticed that this happens especially with either very large or small values, but as I figured still well within eps.
Example:
For UC = [-8.848312431038208e-17] and UT =[0.550226217684170],
x =
Columns 1 through 4
[-0.0000 + 0.3891i -0.0000 - 0.3891i 0.0000 + 0.3891i 0.0000 - 0.3891i]
whereas, checking the answer with "roots" reveals the correct answer, which is reflected by Wolframalpha.
x = roots([a b c d e])
x =
0.9274 + 0.0000i
-0.0000 + 1.0783i
-0.0000 - 1.0783i
-0.9274 + 0.0000i
Is anyone able to tell me whether Matlab rounds off, without me knowing it, or if there is something else that I'm not seeing. I have already checked the code multiple times, but cant seem to be able to find a mistake.
Thank you!
  3 个评论
显示 1更早的评论
John D'Errico
John D'Errico 2018-2-13
编辑：John D'Errico 2018-2-13
MATLAB does not "round off" without telling you. Of course, it does use double precision, so all computations are effectively "rounded off". But that is no different than what happens inside roots. With UC=1e-17 essentially, note that UC^2 is into floating point trash, so some computations will produce garbage.
For example, with the values that you are using:
((UC)^2) + ((UT)^2))
is identical to
((UT)^2))
when computed in double precision. Don't believe me? Try it.
UC = [-8.848312431038208e-17]; UT =[0.550226217684170];
UC^2 + UT^2 == UT^2
ans =
logical
1
More to the point is Matt's question, why do you think this will be usable in an optimizer when roots causes a problem?
Baas28
Baas28 2018-2-13
编辑：Baas28 2018-2-13
Thanks for the feedback! I'm using a third party optimal control software program, GPOPS, that is used for trajectory opimization. When the program starts doing its countless iterations, I couldn't get the roots-function to work (yet) with the gradients that are being plugged back in to it iteratively.
Matt, I agree on the point you make on the point you make, but even with UC = 0, another answer should popup, than it is doing now.
Anyway, as the analytical method was forwarded to me by a prof at uni, and I cross checked it with known values, I trust the method. Therefore, I was wondering if I'm doing something wrong within the Matlab environment. I'll give "roots" a final shot tomorrow though.

请先登录，再进行评论。

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

回答（1 个）

Yassine Damerdji
Yassine Damerdji 2022-4-20
Matlab "roots" function uses the eigenvalue decomposition to derives the roots. Eig function is an iterative procedure that derives the roots in descending order of magnitude. Differences with Ferrari analytical method could be a problem of roundoff errors.

类别

在 Help Center 和 File Exchange 中查找有关 Choose a Solver 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by