[Top][All Lists]

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

Re: How to sort roots?

From: Ben Abbott
Subject: Re: How to sort roots?
Date: Fri, 2 Nov 2007 19:26:29 -0400

On Nov 2, 2007, at 3:22 PM, Schirmacher, Rolf wrote:


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,


While not exactly compliant with your desires, you might find the sort function a good starting point.

p = [1 0 0 1 1];
r = roots(p)
r =

   0.72714 + 0.93410i
   0.72714 - 0.93410i
  -0.72714 + 0.43001i
  -0.72714 - 0.43001i

s = sort(r)
s =

  -0.72714 - 0.43001i
  -0.72714 + 0.43001i
   0.72714 - 0.93410i
   0.72714 + 0.93410i

You may also find the function mpoles useful.

[m,n] = mpoles(s);

ans =

   0.72714 + 0.93410i
   0.72714 - 0.93410i
  -0.72714 + 0.43001i
  -0.72714 - 0.43001i

I realize that each of these examples dance around the solution you desire, but perhaps they can be useful?


reply via email to

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