How to find polynomial roots correctly?

Question

Consider a polynomial such as:

p = [1 -9 27 -27];

obviously the real root is 3:

polyval(p,3)

0

While using the roots function

q = roots([1 -9 27 -27]);

with format short:

q =

   3.0000 + 0.0000i
   3.0000 + 0.0000i
   3.0000 - 0.0000i

and to check if the the roots are real:

bsxfun(@eq,ones(size(q)),isreal(q))

0
0
0

And even worse with format long I get:

roots([1 -9 27 -27])

ans =

  3.000019414068325 + 0.000000000000000i
  2.999990292965843 + 0.000016813349886i
  2.999990292965843 - 0.000016813349886i

How can I calculate roots of a polynomial correctly?

Minor note: your check to see if the roots are real is not correct. isreal(q) gives false if the array q is complex. But some entries may have zero imaginary part. In fact, isreal(q) gives false, whereas for x = q(:).', isreal(x), end gives true, false, false. The first entry of q is real, the others are not, and q as a whole is not real — Luis Mendo

Luis Mendo Luis Mendo · Accepted Answer · 2016-09-01T09:59:10

You may have to work symbolically. You need the Symbolic Math Toolbox for that.

Define the polynomial as a symbolic function. You can (a) use poly2sym to generate the symbolic polynomial from its coefficients. Or (b) better yet, define the symbolic function directly using a string. That way you avoid the loss of accuracy that may result from representing the coefficients as double.
Use solve, which symbolically solves algebraic equations.

Code with option (a):

p = [1 -9 27 -27];
ps = poly2sym(p);
rs = solve(ps);

Code with option (b):

ps = sym('x^3-9*x^2+27*x-27');
rs = solve(ps);

In either case, the result is symbolic:

>> rs
rs =
 3
 3
 3

You may want to convert to numeric values using

r = double(rs);

In your example, this gives

>> format long
>> r
r =
     3
     3
     3

How to find polynomial roots correctly?

3 Answers