Code covered by the BSD License

Highlights from Multiple-root polynomial solved by partial fraction expansion

Be the first to rate this file! 4 Downloads (last 30 days) File Size: 1.62 KB File ID: #22375

Multiple-root polynomial solved by partial fraction expansion

Feng Cheng Chang (view profile)

10 Dec 2008 (Updated )

To find poles/residues of the rational function, instead of roots/multiplicities of the polynomial

File Information
Description

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

Acknowledgements

This file inspired Polynomials With Multiple Roots Solved.

MATLAB release MATLAB 6.5 (R13)
11 Dec 2008

16 Dec 2008

Update the source code. Add two numerical examples, including print-out the complete PRS and related polynomials.

15 Jan 2009

Update the m-file to include both monic-head and monic-tail for polynomial GCD computation.

05 Feb 2009

Update the Polynomial GCD routine.

27 Apr 2009

Update the m-file, to include the overall deviation error of the original polynomial.

29 Apr 2009

Correct typo in m-file