File Exchange

image thumbnail

Polynomials with multiple roots solved

version (1.63 KB) by Feng Cheng Chang
Solving multiple roots polynomial, using simple elementary arithematic operations mostly.


Updated 15 May 2014

View Version History

View License

A given polynomial having multiple roots is solved by the routine
Z = poly_roots(p)
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.

Cite As

Feng Cheng Chang (2021). Polynomials with multiple roots solved (, MATLAB Central File Exchange. Retrieved .

Comments and Ratings (9)

Hans Werner Borchers

When applying poly_roots() to a polynomial with a threefold root in 0.1 and simple roots in 0.5, 0.6, and 0.7 (one of the Jenkins-Traub test examples), the result indicates 0.1 as a simple root. What is going wrong?

>> p = poly([0.1 0.1 0.1 0.5 0.6 0.7])
>> poly_roots(p)
ans =
0.100000 1.000000
0.500018 1.000000
0.599943 1.000000
0.700043 1.000000

Feng Cheng Chang

Reply to pippo: In the 7x2 matrix:
1st column --- polynomial roots
2nd column --- root's multiplicities


Hi, I have a polynomial with 34 orders, but this function gives me roots of 7x2 matrix. Why is that? I copied my polynomial coefficients as below. Could anyone help me out? Thanks!

1 -26.7479290179817 346.353790281681 -2891.79278087362 17494.1298230978 -81695.4463236101 306371.295805817 -947699.866103653 2464717.23435414 -5466000.43064478 10447379.3286537 -17350483.7757060 25192808.9914633 -32131544.3997781 36120714.4912932 -35872219.7954520 31515357.4741505 -24504114.0520695 16854492.5963053 -10241648.1386483 5485377.96325360 -2580928.16918672 1061956.10387350 -379845.324952006 117197.882178068 -30882.3382140184 6860.29267091761 -1262.93321056908 188.263121608582 -21.9982275506094 1.92027369835658 -0.115917764563237 0.00419847428401219 -6.48118467521201e-05

Bor Plestenjak

The correction works for the reported problem, but, there is more. For instance, for p=[3e-04 1e-09 1e-09 5e-08 2e+04] results are wrong and in addition multiplicities appear as complex numbers.

Feng Cheng Chang

Reply to Bor Plestenjak:

Thank you for your pinpoint the bug.

It should be OK if the file is correct as follows: In line 15, change
" if npo == 1, return, end;"
" if npo == 0, return, end;"

I will update this file later. But please try again now!

Bor Plestenjak

There seems to be a bug when it comes to linear polynomials (which of course can not have multiple zeros). poly_roots([1 2]) does not return a zero and poly_roots([1 1 0 0 0]) returns just the triple zero 0.

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.

MATLAB Release Compatibility
Created with R13
Compatible with any release
Platform Compatibility
Windows macOS Linux

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!