help-octave
[Top][All Lists]
Advanced

[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

Hello,

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
data

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
is 

  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)];
  endif 

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

  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,

Rolf



reply via email to

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