help-octave
[Top][All Lists]

## 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()?