2

I have a question about solving quartic equation in MATLAB, please see here

The MATLAB function shown as below:

MATLAB

function [res] = solveQuartic(a, b, c, d, e)
 p = (8*a*c - 3*b^2)/(8*a^2);
 q = (b^3 - 4*a*b*c + 8*a^2*d)/(8*a^3);

 delta0 = c^2-3*b*d + 12*a*e;
 delta1 = 2*c^3 - 9*b*c*d + 27*b^2*e + 27*a*d^2 - 72*a*c*e;
 Q = ((delta1 + sqrt(delta1^2 - 4*delta0^3))*0.5)^(1/3);
 S = 0.5*sqrt(-2/3*p + (Q + delta0/Q)/(3*a));

 x1 = -b/(4*a) - S - 0.5*sqrt(-4*S^2-2*p + q/S);
 x2 = -b/(4*a) - S + 0.5*sqrt(-4*S^2-2*p + q/S);
 x3 = -b/(4*a) + S - 0.5*sqrt(-4*S^2-2*p - q/S);
 x4 = -b/(4*a) + S + 0.5*sqrt(-4*S^2-2*p - q/S);
 res = [x1; x2; x3; x4];
end

TEST

solveQuartic(-65, -312, -582, -488, -153)

(*
-1.400000000000000 - 0.627571632442189i
-1.000000000000000 + 0.000000000000000i
-1.400000000000000 + 0.627571632442189i
-1.000000000000000 + 0.000000000000000i
*)

enter image description here

Obviously, the quartic equation $$-65x^4-312x^3-582x^2-488x-153=0$$ owns two real roots

Now I would like to exact the real roots of the quartic equation, my trial as below

function [res] = solveQuarticReals(a, b, c, d, e)
 p = (8*a*c - 3*b^2)/(8*a^2);
 q = (b^3 - 4*a*b*c + 8*a^2*d)/(8*a^3);

 delta0 = c^2-3*b*d + 12*a*e;
 delta1 = 2*c^3 - 9*b*c*d + 27*b^2*e + 27*a*d^2 - 72*a*c*e;
 Q = ((delta1 + sqrt(delta1^2 - 4*delta0^3))*0.5)^(1/3);
 S = 0.5*sqrt(-2/3*p + (Q + delta0/Q)/(3*a));

 x1 = -b/(4*a) - S - 0.5*sqrt(-4*S^2-2*p + q/S);
 x2 = -b/(4*a) - S + 0.5*sqrt(-4*S^2-2*p + q/S);
 x3 = -b/(4*a) + S - 0.5*sqrt(-4*S^2-2*p - q/S);
 x4 = -b/(4*a) + S + 0.5*sqrt(-4*S^2-2*p - q/S);
 res = [x1; x2; x3; x4];
 %exact the real roots
 j = 0;
  for i = 1 : length(res)
   if imag(res(i)) == 0
    j = j + 1;
    mid(j) = res(i); 
   end
  end
  res = mid;
end

However, it failed.

QUESTION

  • How exact the real root in MATLAB?
  • Or is there a built-in like Chop in Mathematica?

enter image description here

Community
  • 1
  • 1
xyz
  • 413
  • 2
  • 19
  • This is how a question should be, with complete code, and with an explanation of what you have tried, what you expect and what you get! Kudos! =) – Stewie Griffin Sep 24 '15 at 08:19

1 Answers1

3
a=-65;
b=-312;
c=-582;
d=-488;
e=-153;


p = (8*a*c - 3*b^2)/(8*a^2);
q = (b^3 - 4*a*b*c + 8*a^2*d)/(8*a^3);

delta0 = c^2-3*b*d + 12*a*e;
delta1 = 2*c^3 - 9*b*c*d + 27*b^2*e + 27*a*d^2 - 72*a*c*e;
Q = ((delta1 + sqrt(delta1^2 - 4*delta0^3))*0.5)^(1/3);
S = 0.5*sqrt(-2/3*p + (Q + delta0/Q)/(3*a));

x1 = -b/(4*a) - S - 0.5*sqrt(-4*S^2-2*p + q/S);
x2 = -b/(4*a) - S + 0.5*sqrt(-4*S^2-2*p + q/S);
x3 = -b/(4*a) + S - 0.5*sqrt(-4*S^2-2*p - q/S);
x4 = -b/(4*a) + S + 0.5*sqrt(-4*S^2-2*p - q/S);
res = [x1; x2; x3; x4];

%exact the real roots  
realNumber = real(res((abs(imag(res)) <= 1e-10)))

Please try these codes. Because of the characteristics of floating point data type, there could be a little difference between imag(S) and imag(0.5*sqrt(-4*S^2-2*p + q/S)), even if two values are the same in the mathematics. In order to check the error, please type "format shortEng" on your matlab window, and then type imag(S) - imag(0.5*sqrt(-4*S^2-2*p + q/S)). The answer is 55.5112e-018i. The imagnary number of answer is definitely not 0, but 55.5112e-018. This phenomenon is caused by the nature of floating point data type.

PS. Thank you very much Stewie Griffin!! :-)

Stewie Griffin
  • 14,889
  • 11
  • 39
  • 70
KKS
  • 1,389
  • 1
  • 14
  • 33
  • 1
    No, no, NO! Don't use growing vectors/matrices in MATLAB. Skip the loop and use logical indexing instead! `realNumber = res(real(abs(imag(res)) <= 1e-10))`. It might be I made a mix up, I haven't tested it, but you get the picture. You can also do (quite common): `mask = real(abs(imag(res)) <= 1e-10); realNumber = res(mask);` I won't post this as an answer, as the rest of your answer is nice =) (Sorry if that was a bit harsh, it was not directed at you personally, but towards the very common mistake of having growing matrices (`A = [A, res]`). – Stewie Griffin Sep 24 '15 at 08:06
  • 1
    BTW! **Please** continue answering questions! Your contributions are very much appreciated. It's not my intention to scare you away... It's just a very common mistake, and your code will become _very_ slow if you use this approach, the MATLAB editor even warns you about it. You have my update, but I hope you take the time to edit the answer. =) I see I forgot at least one parentheses in my comment above, nut you figure it out I guess =) – Stewie Griffin Sep 24 '15 at 08:13
  • Many thanks for your comments. I will try to update my code, according to your advice. :-) Thank you~ – KKS Sep 24 '15 at 08:19