File Exchange

image thumbnail

Multiple-root polynomial solved by partial fraction expansion

version 1.11.0.0 (1.62 KB) by Feng Cheng Chang
To find poles/residues of the rational function, instead of roots/multiplicities of the polynomial

2.4K Downloads

Updated 29 Apr 2009

View License

A given polynomial p(x) is transformed into a rational function r(x). The poles and residues of the derived rational function are found to be equivalent to the roots and multiplicities of the original polynomial.
p(x) = Given polynomial
= PROD[k=1:K]{(x - z_k)^m_k}
d(x) = (d/dx)p(x)
g(x) = GCD(p(x),d(x))
u(x) = p(x)/g(x)
w(x) = (d/dx)u(x)
v(x) = d(x)/g(x)
r(x) = v(x)/u(x)
= SUM[k=1:K]{m_k/(x - z_k)}
Thus, the roots z_k are computed from solving the simple-root polynomial u(x)=0, instead of the original multiple-root polynomial p(x)=0; and the multiplicities m_k are determined as the partial fraction expansion coefficients of the derived rational function r(x)=v(x)/u(x),
z_k = Roots(u(x)), k=1,K
m_k = v(z_k)/w(z_k), k=1,K
In addition, re-constructing a polynomial pz(x) from the computed z_k and m_k, the overall deviation error of the original polynomial p(x) is calculated,
er = Norm(pz - p)/Norm(p)
The polynomial GCD is calculated from "Monic polynomial subtraction" derived from the longhand polynomial division in classical Euclidean GCD algorithm. It requirs only simple algebric operations without any high mathematics.
The source code contains total of only 43 lines, using merely basic built-in MATLAB functions, and applying only existing double precision. Amazingly, it gives the expected results of test polynomials of very high degree , such as
p(x) = (x - 123456789)^30
p(x) = (x + 100)^20 * (100x-1)^10
p(x) = (x+1)^40 * (x-2)^30 * (x+3)^20 * (x-4)^10
p(x) = (x + 1)^1000

______________________________________________________

The code is list here for reader's convenince. (only 43 lines)

function [zm,er] = polyroots(p)
% *** A polymonial with multiple roots ***
% Solved via partial fraction expansion
d = polyder(p);
g = polygcd(p,d);
u = deconv(p,g);
v = deconv(d,g);
w = polyder(u);
z = roots(u);
m = round(abs(polyval(v,z)./polyval(w,z)));
zm = [z,m]; % p,d,g,u,v,w,z,m,zm
pz = polyget([m,-z,ones(length(z),1)])*p(1);
er = norm(pz-p)/norm(p); % pz,er
function g = polygcd(p,q)
% *** GCD of a pair of polynomials ***
% by "Monic polynomial subtraction"
n = length(p)-1; nc = max(find(p))-1;
m = length(q)-1; mc = max(find(q))-1;
nz = min(n-nc,m-mc);
if nc*mc == 0, g = [1,zeros(1,nz)]; return, end;
p2 = [p(1:nc+1)];
p3 = [q(1:mc+1)];
for k = 1:nc+nc,
p3 = [p3(min(find(abs(p3)>1.e-6)): max(find(abs(p3)>1.e-6)))];
p1 = [p2/p2(1)]; % k,p1,
p2 = [p3/p3(1)];
p3 = [p2,zeros(1,length(p1)-length(p2))]-[p1,zeros(1,length(p2)-length(p1))];
if norm(p3)/norm(p2) < 1.e-3, break; end;
end;
g = [p1,zeros(1,nz)];
function p = polyget(A)
% *** A polynomial coefficient vector from sub-polynomial factors ***
p = 1;
for i = 1:length(A(:,1)),
q = 1;
for j = 1:A(i,1),
q = conv(q,A(i,max(find(A(i,:))):-1:2));
end;
p = conv(p,q);
end;

_____________________________________________________

Typical Numerical Example:

>> % Contruct a test polynomial:
>> p = poly([ 1 1 1 1 1 1 1 -1 -1 -1 -1+2i -1+2i -1+2i -1-2i -1-2i -1-2i 2 2 3 3 +i +i +i -i -i -i -3 0 0 0 0 0 ])
p =
1 -5 2 -6 76 140 -802 954 -4251 13663 -18740 28472 -53504 45776 5212 -77580 185243 -220631 104794 52458 -193356 248612 -146266 9202 65791 -87555 55800 -13500 0 0 0 0 0
>> % Roots and multiplicities for the polynomial are computed.
>> zm = polyroots(p)
zm =
3.0000 -------------- : 2.0000
-3.0000 -------------- : 1.0000
-1.0000 + 2.0000i ---: 3.0000
-1.0000 - 2.0000i ---: 3.0000
2.0000 ---------------: 2.0000
-1.0000 ---------------: 3.0000
-0.0000 + 1.0000i ---: 3.0000
-0.0000 - 1.0000i ---: 3.0000
1.0000 ---------------: 7.0000
0.0000 ---------------: 5.0000

Cite As

Feng Cheng Chang (2021). Multiple-root polynomial solved by partial fraction expansion (https://www.mathworks.com/matlabcentral/fileexchange/22375-multiple-root-polynomial-solved-by-partial-fraction-expansion), MATLAB Central File Exchange. Retrieved .

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!