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?

Answer 1

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

1. 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.

2. 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

Answer 2

This is due to floating point inaccuracies. Have a look at this post for details: Is floating point math broken?

One thing you can do is to round off the answer/s upto some decimal places like this:

q = round(roots([1 -9 27 -27]), 4) % rounding off to 4 decimal places