On Mon, 06 Jul 2009 16:53:42 +0200 Fabrice Silva <[email protected]> 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. > _______________________________________________ Numpy-discussion mailing list [email protected] http://mail.scipy.org/mailman/listinfo/numpy-discussion
