Geek Answers Handbook

How to find polynomial roots?

Consider a polynomial such as:

p = [1 -9 27 -27];

obviously the real root is 3:

 polyval(p,3) 0

When 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 check the correctness of the roots:

 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 correctly calculate the roots of a polynomial?

+6
source share
3 answers

This is due to floating point inaccuracies. Take a look at this post for details: Could floating point math fail?

One thing you can do is to round the answer to a few decimal places as follows:

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

You may need to work symbolically. For this you need the Symbolic Math Toolbox.

  • Define a polynomial as a symbolic function. You can (a) use poly2sym to create a character polynomial from its coefficients. Or (b) even better, define a symbolic function directly using a string. Thus, 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 any case, the result is symbolic:

 >> rs rs = 3 3 3

You can convert to numerical values ​​with

 r = double(rs);

In your example, this gives

 >> format long >> r r = 3 3 3
+6
source

This is very specific to your polynomial. In general, you should expect that the root of multiplicity m has a relative floating point error of mu^(1/m) , where mu=1e-15 is the accuracy of the machine. In this case, the multiplicity is m=3 , so the error is in the range of 10^(-5) . This is exactly the scale of the errors in your results.

The fact that this happens here with explicit integer coefficients is the result of using the matlab method, it calculates the eigenvalues ​​of the matrix of related devices, and the eigenvalue algorithm converts the integer matrix into its own floating point matrix with the corresponding rounding errors in the first step of the algorithm.

Other algorithms have empirical tests for multiplicities and related clusters of approximative roots and, thus, are able to correct this error. In this case, you could achieve this by replacing each root with an average of three roots.

Mathematically you have some polynomial

 p(x)=(xa)^m*q(x)

with a root in x=a multiplicity m . Due to floating point operations, the solver effectively "sees" the polynomial

 p(x)+e(x)

where the coefficients e(x) have a size that is the value of the coefficients p times mu . Near the root a this perturbed polynomial is effectively replaced by

 (xa)^m*q(a)+e(a) = 0 <==> (xa)^m = -e(a)/q(a)

so that the solutions form a m-oriented regular polygon or star centered at a with radius |e(a)/q(a)|^(1/m) , which should be in the region |a|*mu^(1/m)

0
source

More articles:


All Articles