[Numpy-discussion] roots and high-order polynomial
Charles R Harris
Mon Jul 6 09:32:58 CDT 2009
On Mon, Jul 6, 2009 at 8:16 AM, Charles R Harris
> On Mon, Jul 6, 2009 at 3:44 AM, Fabrice Silva <firstname.lastname@example.org>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
>> 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
> 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.
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the NumPy-Discussion