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