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