[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.
Nils
> 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.
>
More information about the NumPy-Discussion
mailing list