On Mon, Jul 6, 2009 at 8:16 AM, Charles R Harris <[email protected]>wrote:
> > > On Mon, Jul 6, 2009 at 3:44 AM, Fabrice Silva <[email protected]>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
_______________________________________________ Numpy-discussion mailing list [email protected] http://mail.scipy.org/mailman/listinfo/numpy-discussion
