Re: roots accuracy

[Nicholas Jankowski]

On Oct 31, 2015 9:28 AM, "Doug Stewart" <address@hidden> wrote:
>
>  Recently these was a user who tried to find the roots of:
> a=poly([ -2 -2 -2 ])
> a =
>
>     1    6   12    8
>
> and our normal root finding code gives:
>
> >> roots(a)
> ans =
>
>   -1.99998027897029 + 0.00000000000000i
>   -2.00000986051485 + 0.00001707909463i
>   -2.00000986051485 - 0.00001707909463i
>
> Matlab also gives the same numbers. So we are compatible with Matlab.
>
> But I was sure that we could do better and JordiGH suggested:
> http://math.cts.nthu.edu.tw/Mathematics/RANMEP%20Slides/Zhonggang%20Zeng.pdf
>
> After browsing through that article (but understanding only a small part of it)
> I could see that my ideas were similar.
>
> The idea is that
>    1)  repeated roots are always difficult.
>    2)  when you use the eigen value method to find the repeated roots
>         it will find a set that are equally spaced, in a circle, around the actual roots 
>         in the complex plane.
>    3)  the mean of these equally spaced answer is always much closer
>         to the correct answer.
>
> My method:
>  1) use roots to find the roots as we do now.
>  2) use mpoles to find if there are repeated roots
>  3)  if there are repeated roots then
>     4)  factor out all roots with a multiplicity of 1
>     5) find the group of roots with the highest multiplicity (M)
>     6) find the mean of this group and us that M times in the list of roots
>     7) factor out these roots
>     8) repeat steps 5 -7 until all roots have been found.
>  9) put all the different roots together in a vector.
>
> Using this method the results are:
>
> a=poly([ -22 -22 -22 -4 -6 -2.01 -2.01]);
> # using old method
> roots(a)
> ans =
>
>   -21.99990711577122 +  0.00016087892205i
>   -21.99990711577122 -  0.00016087892205i
>   -22.00018576845768 +  0.00000000000000i
>    -5.99999999999993 +  0.00000000000000i
>    -4.00000000000002 +  0.00000000000000i
>    -2.01000015603606 +  0.00000000000000i
>    -2.00999984396397 +  0.00000000000000i
>
> # now my new method
> newroots(a)
> ans =
>
>   -4.00000000000002
>   -5.99999999999995
>   -22.00000000000001
>   -22.00000000000001
>   -22.00000000000001
>    -2.01000000000000
>    -2.01000000000000
>
>
> I have tried many different polynomials with similar results.
>
> I would like the Octave teem's feedback  on where to go with this.
>  1) We could make a separate file and advise users that  if they 
>       have multiple roots to use newroots.
>  2) We could have roots itself advise this, if it found a multiplicity greater than 1
>  3) We can modify roots to automatically do this
>  4) We can add an option to roots to do this if it is selected.
>
>  Some of these options keep compatibility with matlab and others don't 
>
> I know that I would always use this new way, but then my answers will not match Matlab's answers.
>
> Feedback, suggestions are most welcome.
> Doug
>
> PS I have a proof of concept file running, but there is much work yet to finish this code.
>
>
>  
> --
> DAS
>

Personally I like the optional setting. Add a ...'Method','NewRoots') or similar. One less function to track and document. Does roots for octave or MATLAB currently take options? Thinking of 2-way compatibility, do MATLAB functions throw a usage error if called with extra arguments?