[Numpy-discussion] roots and high-order polynomial

[Fabrice Silva]

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.