help-octave
[Top][All Lists]

## 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:

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

Rolf

```
```
```
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);

s(n)
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?
```
Ben

```