[Numpy-discussion] roots and high-order polynomial

[Charles R Harris]

On Mon, Jul 6, 2009 at 8:16 AM, Charles R Harris
<charlesr.harris at gmail.com>wrote:

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

I should mention that the computation of the eigenvalues of non-Hermitean
matrices is numerically difficult in itself, much more so than for Hermitean
matrices. The companion matrix approach is a convenient way to get the roots
but I wouldn't trust it very far. Pauli's suggestion of using Chebyshev
polynomials instead of power series might yield a better conditioned matrix.
I think his suggestion worth pursuing if you can make the adjustment.

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