[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: roots accuracy
From: |
Nicholas Jankowski |
Subject: |
Re: roots accuracy |
Date: |
Sat, 31 Oct 2015 10:04:03 -0400 |
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?
- roots accuracy, Doug Stewart, 2015/10/31
- Re: roots accuracy,
Nicholas Jankowski <=