[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

How to sort roots?

From: Schirmacher, Rolf
Subject: How to sort roots?
Date: Fri, 2 Nov 2007 20:22:03 +0100


I have a problem that includes the roots of a forth order polynomial, so
with r being the positive real and i being the positive imaginary part of
the solution, they shoud come out as

-r + i
-r - i
r + i 
r - i

For example have a look at 

octave.exe> roots([1 0 0 0 1])
ans =

  -0.707106781186548 + 0.707106781186548i
  -0.707106781186548 - 0.707106781186548i
   0.707106781186548 + 0.707106781186548i
   0.707106781186548 - 0.707106781186548i

for a simple case and 

octave.exe> roots([1 0 0 0 1 0 0 0 1e10])
ans =

  -16.42915099372338 +  6.80520121982268i
  -16.42915099372338 -  6.80520121982268i
   -6.80520121982269 + 16.42915099372340i
   -6.80520121982269 - 16.42915099372340i
    6.80520121982269 + 16.42915099372338i
    6.80520121982269 - 16.42915099372338i
   16.42915099372337 +  6.80520121982269i
   16.42915099372337 -  6.80520121982269i

showing a simple noisy example with two fourth order solutions.

As my polynomial is of higher order (currently about 40) and less well
conditioned, there is some numerical noise in the solution I get from roots
and it has the form of 

-r1 + i1
-r1 - i1
+r2 + i2
+r2 - i2

with r1, r2 and i1, i2 some slightly different positive numbers.

Now, I automatically want to extract the roots with the largest real part
and those with the largest imaginary part and I know that I will always have
second order solutions, i.e. +/- (cplx. solution), so I can really look for
the largest real / imag part and then have +/- the two distinct solutions.
Sometimes, the ploynomial degrades to the 4th order solution and then due to
the numerical noise, there are some difficulties:

Sorting for the largest real part is o.k., it gives the r2-i2 solution as
sort keeps the order:
c_n2(f_index,:) contains the 40 roots of a polynomial and

  [s,idx] = sort(- (real(c_n2(f_index,:))));  # sort idx according to the
real part, 
                                              # largest real part on top
  c_n2_sorted(f_index,:) = c_n2(f_index,idx); # assemble sorted original

Sorting for the largest imag part in the same way is o.k. IFF there are only
2nd order roots:

  [s,idx] = sort(- (imag(c_n2(f_index,:)))); # sort idx according to the
imag part 

BUT if I get a 4th order solution, then it depends on the numerical noise
whether i1 or i2 is larger and I get -r1+i1 or +r2+i2. I want to get +r2+i2
because that is the second, lineary independent solution while -r1+i2 is
just minus the first solution (except of the numerical noise). So what I do

  if(abs(c_n2(f_index,idx(1)) + c_n2_sorted(f_index,1)) < \
     (abs(c_n2(f_index,idx(1))) * eps)*20)
    ##    disp(" Index shift detected ");
    idx = [shift(idx(1:2),1) idx(3:end)];
    idx = [idx(1:end-2) shift(idx(end-1:end),1)];

to detect the "wrong", linear dependent solution and shift the index and

  c_n2_sorted_imag(f_index,:) = c_n2(f_index,idx); # assemble sorted
original data

20*eps*abs(quantity) seems to be a good estimate for the "numerical noise"
here, at least for my problem.
f_index is an index running from 1 to 1000 (typically).

Now, is there any more convenient / efficient way to handle such an issue? 
Is there any more robust method less depending on the (undetermined?)
ordering of roots from roots()? 

Thanks in advance,


reply via email to

[Prev in Thread] Current Thread [Next in Thread]