[Numpy-discussion] roots and high-order polynomial

[Charles R Harris]

On Mon, Jul 6, 2009 at 3:44 AM, Fabrice Silva <silva at lma.cnrs-mrs.fr> wrote:

> Le vendredi 03 juillet 2009 à 10:00 -0600, Charles R Harris a écrit :
>
> > What do you mean by erratic? Are the computed roots different from
> > known roots? The connection between polynomial coefficients and
> > polynomial values becomes somewhat vague when the polynomial degree
> > becomes large, it is numerically ill conditioned.
>
> For an illustration of what I mean by 'erratic', see
> http://fsilva.perso.ec-marseille.fr/visible/script_als_nbmodes.png or
> http://fsilva.perso.ec-marseille.fr/visible/script_als_nbmodes.pdf
> with the real part of the roots on the left, and the imaginary part of
> the right:
> - for low orders (<26), a pair of complex conjugate roots appears when
> the order of polynomial is increased (with a step equal to 2). As
> expected in my physical application, their real parts are negative.
> - for higher order (>=26), there is no more 'hermitian
> symmetry' (complex conjugate solutions), and real parts become
> positive...
>

Double precision breaks down at about degree 25  if things are well scaled,
so that is suspicious in itself. Also, the companion matrix isn't Hermitean
so that property of the roots isn't preserved by the algorithm. If it were
Hermitean then eigh would be used instead of eig. That said, there are other
ways of computing polynomial roots that might preserve the Hermitean
property, but I don't know that that would solve your problem.


>
> The computation of the coefficients of the polynomial is out of topic, I
> already checked it and there is no errors.
>

The problem is floating point round off error in representing the
coefficients. Even if you know the coefficients exactly they can't generally
be represented exactly in double precision. Any computational roundoff error
just adds to that. If the coefficients were all integers I would have more
confidence in the no error claim.

Where do the coefficients come from? Perhaps there are options there.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20090706/8328dab9/attachment.html>