Polynomials with multiple roots solved

Solving multiple roots polynomial, using simple elementary arithematic operations mostly.

Description

A given polynomial having multiple roots is solved by the routine

              Z = poly_roots(p)
where
            Input p: polynomial coefficient vector
            Output Z: root-multiplicity pairs

The MATLAB source code is very simple and compact (fewer then 50 lines) and amazingly gives the expected results for any test polynomials of very high degree and multiplicities.

In this simple routine 'poly_roots.m', except for a MATLAB built-in function 'roots.m', all the required computations involve only simple elementary arithematic operations (such
as additions, subtractions, multiplications, divisions, and integer exponations), without applying any high mathematics. The only routine 'roots.m' is used here mainly to solve a simple roots polynomial (not any multiple roots polynomials).

Most of the calculations in this routine involve elementary arithematic operations, therefore, it is fairly easy to improve the expected results from the existing double precision to vpi' multiple precision.
 
For detail description, please refer to
 
F.C. Chang, "Solving multple-root polynomials" IEEE Antennas and Propagation Magazine, Vol.51, No.6, pp. 151-155, Dec. 2009.

Acknowledgements

Gcd Of Polynomials, Multiple Root Polynomial Solved By Partial Fraction Expansion, and Polynomial Coefficient Vector Derived From Sub Polynomial Factors inspired this file.

This file inspired Solving Multiple Root Polynomials.

• MATLAB release: MATLAB 6.5 (R13)

[John D'Errico]

Well done.

[Andrew Mendelsohn]

Nice little file finds multiple roots much more accurately than built in routine. Used it in an inverse Laplace transform application.

[Feng Cheng Chang]

Self comment:

Amazingly, the revised routine can find the desired roots of test polynomials p(x) -- expanded, such as
p(x) = (x+987654321)^N,
where N=any positive integers, such as 10, 100, 1000, 10000, ....

For example, find the roots of the following polynomial expanded
p(x) = (x+987654321)^12345.

>> format long g
>> p = poly([-987654321*ones(1,12345)]);
>> Z = poly_roots(p)
Z =
-987654321 12345

It takes about 2 min total time. Most is to create polynomial coefficient vector p, and only 0.5 sec to get the desired root-multiplicity pairs Z.

Of course, any numerical algorithm is not fool prove, neither is mine.

Foe example,it does work for p(x)=(9876x+12345)^70, but fails for p(x)=(9876x+12345)^80.

I hope someone would like to try it for some other existing test polynomials, and give me any valuable comments.