Le lundi 06 juillet 2009 à 08:16 -0600, Charles R Harris a écrit : > 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 Hermitian 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 think there is a misunderstanding: I was referring to the fact the solution had to be real or complex conjugate, and not that the companion matrix would be a hermitian matrix (which it isn't due to its construction). Something I forgot to tell is that the roots are complex eigenfrequencies of a physical system, the real and imaginary parts expressing the damping and the frequency, respectively. If a positive frequency is solution then the opposite negative is solution too but with the «same-signed» damping. > 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. Here is the construction: the coefficients are obtained from a modal analysis of a subsystem of a bigger system. A quantity called impedance depending of a variable X is the result of the combination of several terms: Z(X) = C1/(X-X1)+C1*/(X-X1*)+...+CN/(X-XN)+CN*/(X-XN*) where * denotes the complex conjugate. I have to get the solutions X of an equation Ze(X)=N(X)/D(X) with N and D are very-low order polynomial (orders 1 and 2). An alternative of using iterative root finding (for example with scipy.optimize.fsolve) is to expand the equation as a polynomial of order close to 2N. I do understand this approach might seem non-obvious but it is rather efficient for a low value of N (<10)... -- Fabrice Silva Laboratory of Mechanics and Acoustics - CNRS 31 chemin Joseph Aiguier, 13402 Marseille, France. _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion