On Sun, Nov 1, 2015 at 12:23 PM, Lukas Reichlin <address@hidden> wrote:
@Lukas That bug I think is more about how rlocus chooses a range for k.Hi Doug,
> On 31.10.2015, at 15:04, Nicholas Jankowski <address@hidden> wrote:
>
> 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?
>
I’m interested in your concept and would prefer the optional setting as well. I’ve just seen bug #46325 <http://savannah.gnu.org/bugs/?46325> and thought that your concept could be a part of the solution.
Best regards,
Lukas
In this case of 3 repeated roots and no other roots it chooses a small range for k.