Re: Finding only some roots of a polynomial

[Jaroslav Hajek]

On Fri, Oct 9, 2009 at 1:58 AM, Pablo <address@hidden> wrote:
> Hi everyone!
>
> I need to find the roots of some quite big polynomials (order ~800). The
> root function works fine (and surprisingly fast!), but I don't think
> it's the best tool forwhat I need. I'm only interested in those roots
> with high absolute value, so after "roots" works hard to find the 800
> complex values, I discard almost all of them.
>
> Do you have any suggestion for an algorithm to do that?
>
> I did find a book on the subject (McNamee - Numerical Methods for Roots
> of Polynomials, Vol1 - 2007), and it suggest (I didn't read it
> thoroughly, yet) that the companion method finds the roots one by one,
> starting with the one with the largest magnitude. However, there are a
> multitude of similar methods in the book.
>
> Is Octave doing that?
>
> In that case I might try to modify the function to stop after a few
> roots. I peeked in the source code, but I don't know Fortran, so it's
> kind of hard to understand what is going on.
>
> For reference, my current script is
>  %lenght(c) is ~800
>  r=roots(c)'; %This takes quite long
>  %I only keep roots with abs(r)>=0 and angle(r)>=0
>  P=abs(r);
>  I=find(P>=1);
>  BB=log(P(I))/deltat;
>  WW=angle(r(I))/deltat;
>  I=find(WW>=0);
>  B=BB(I);
>  W=WW(I);
>  %Now lenght(W) is just around 10
>
> Thank you very much for your answers, and for maintaining this great
> program!
>  Arnoques

Octave takes a simple approach: it constructs a companion matrix of
the polynomial whose eigenvalues match the polynomial roots, then uses
an eigenvalue solver.
You can construct the companion matrix using "compan".
So, if you have Octave linked with ARPACK (so that you have the eigs
function available), you can do sth like:

C = compan (c);
eigs (C, 3) # find 3 largest roots.

hth

-- 
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz