Re: problem with consistency of root finding and sorting

[Ben Abbott]

On Oct 9, 2007, at 8:23 PM, David Bateman wrote:

> Ben Abbott wrote:
> > I have Octave 2.9.14 on Mac OS X (intel & PPC) installed via Fink and
> > am concerned that my installations are giving different results than
> > others.
> >
> > During the process of modifying residue.m so that it would function
> > in a manner consistent with Matlab's I've become suspicious that my
> > Octave installation  does not work the same as others.
> >
> > The difference is small, but is troublesome.
> >
> > Given the following short example,
> >
> >         r = sort(roots([1 0 18 0 81]))
> >         r_err = r - round(r)
> >
> > I obtain the results below
> >
> > r = sort(roots([1 0 18 0 81]))
> > r =
> >
> >     0.0000 - 3.0000i
> >     0.0000 + 3.0000i
> >    -0.0000 - 3.0000i
> >    -0.0000 + 3.0000i
> >
> > r_err = r - round(r)
> > r_err =
> >
> >    2.5814e-08 + 1.6920e-08i
> >    2.5814e-08 - 1.6920e-08i
> >    -2.5814e-08 - 1.6920e-08i
> >    -2.5814e-08 + 1.6920e-08i
>
>
> Under Octave 2.9.14 compiled with g++ 4.0.1 on a Mandriva 2006  
> machine I get
>
> octave:3> r = sort(roots([1 0 18 0 81]))
> r =
>
>    0.0000 - 3.0000i
>    0.0000 + 3.0000i
>   -0.0000 - 3.0000i
>   -0.0000 + 3.0000i
>
> octave:4> r_err = r - round(r)
> r_err =
>
>   1.5973e-08 + 6.5742e-09i
>   1.5973e-08 - 6.5742e-09i
>   -1.5973e-08 - 6.5742e-09i
>   -1.5973e-08 + 6.5742e-09i
>
> D.

The order of your roots is the same as mine! I've attached two files,  
residue.m and mpoles.m. Perhaps the small differences in r_err are  
due to slight differences in BLAS? <shrugs>

David, if you have the time, please download /scripts/polynomial/ 
residue.m and /scripts/polynomial/mpoles.m from the cvs, and try  
another test for me.

b = [1 0 1];
a = [1 0 18 0 81];
[r, p, k, e] = residue(b, a)

r =

  -0.00000 - 0.09259i
   0.22222 - 0.08333i
   0.00000 + 0.09259i
   0.22222 + 0.08333i

p =

   0.0000 + 3.0000i
   0.0000 + 3.0000i
  -0.0000 - 3.0000i
  -0.0000 - 3.0000i

k = [](0x0)
e =

   1
   2
   1
   2