[Numpy-discussion] roots and high-order polynomial

Nils Wagner nwagner@iam.uni-stuttgart...
Mon Jul 6 10:13:10 CDT 2009

On Mon, 06 Jul 2009 16:53:42 +0200
  Fabrice Silva <silva@lma.cnrs-mrs.fr> wrote:
> 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). 

IIRC, the coefficients of your polynomial are complex.
So, you cannot guarantee that the roots are complex 
conjugate pairs.


> 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. 
>> computational roundoff error just adds to that. If the 
>> 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 
> 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.

More information about the NumPy-Discussion mailing list