Issues with solving a quartic function analytically

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!

[Matt J]

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.

Please clarify this part. Are you really sure all this is necessary?

[John D'Errico]

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]

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.

Answer 1

[Yassine Damerdji]

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.