|
"Grzegorz Knor" wrote in message <iomh7h$j76$1@fred.mathworks.com>...
> Hi!
>
> Is there a way to avoid this type of errors in matlab:
>
> p = poly([1 1 1 1 1 1 1])
> p =
> 1 -7 21 -35 35 -21 7 -1
> r = roots(p)
> r =
> 1.0090
> 1.0056 + 0.0070i
> 1.0056 - 0.0070i
> 0.9980 + 0.0088i
> 0.9980 - 0.0088i
> 0.9919 + 0.0039i
> 0.9919 - 0.0039i
>
> ?
>
> Grzegorz
- - - - - - - - - -
From a numerical analysis point of view you have given 'roots' a difficult problem, Grzegorz. Your polynomial is y = (x-1)^7 which means that if y differs from 0 by only 1.0e-16, the corresponding x would differ from 1 by a whopping big 0.005, which as you see is of the order of error you are experiencing. Since 'roots' must use the polynomial form to compute y, such errors in y are inevitable.
Suppose your roots were all to be pi as in y = (x-pi)^7. If you expand this polynomial to powers of x, and round each of the resulting coefficients to the 53-bit accuracy of which your computer is capable, then I believe you will find that the theoretical roots of infinite precision for such a "rounded" polynomial would depart from pi by very substantial amounts, of the same order of magnitude that you have observed in your trial - that is, more or less in the third decimal place. It is a problem that is inherent in the mathematical nature of finding repeated polynomial roots.
Of course you might argue that in your particular case with all roots equal to a common integer, the polynomial coefficients would necessarily all be exact integers also so there should be no reason to make any errors whatever. I am afraid 'roots' is not so sophisticated as to be able to take advantage of such a special situation as that.
Roger Stafford
|
Feed for this Thread