Code covered by the BSD License

# An Easy Efficient Pure Algebraic Method for Partial Fraction Expansion

### MA Youneng (view profile)

A matlab code for partial expansion of rational functions with multiple high-order poles.

pfe.m
```%_______________________________________
% This a matlab code for partial fraction expansion of rational functions
% The related article can be found at
%______________________________________
%input:
%S : array of poles
%m : array of multiplitities of poles
%D : polynominal coefficients of numerator in ascending order
%s0: a constant coefficient of the numerator

%output:
%res   matrix residues
%E:      coefficients of residueal polynomial in ascending order
%Author: Younenng MA, Fudan University
%data  : June 24, 2013

function [res,E]=pfe(S,m,D,s0)

M=length(S);  %number of poles
m_max=max(m); %largest multiplicity
if M<2 ||M~=length(m);
res=0;
E=0;
disp('wrong input !');
return;
end
%A is connected with bionominal coefficients
A=zeros(m_max);
A(:,1)=1;
A(1,:)=1;
for i=2:m_max
for j=2:m_max
A(i,j)=A(i-1,j)+A(i,j-1);
end
end
B=zeros(M);
for i=2:M
for j=1:i-1
B(i,j)=1/(S(i)-S(j));
B(j,i)=-B(i,j);
end
end
%
%______________________________________
C_0=zeros(M+2,m_max+2);
%____________
%Exanding R_0^2(0)
for j=1:m(1)
C_0(1,j)=(-1)^(m(1)-j)*A(m(2),m(1)-j+1)*B(1,2)^(m(1)+m(2)-j);
end
for j=1:m(2)
C_0(2,j)=(-1)^(m(2)-j)*A(m(1),m(2)-j+1)*B(2,1)^(m(1)+m(2)-j);
end
%__________
%Exanding R_0(0)
for n=3:M
C_1=zeros(M+2,m_max+2);
%
for i=1:n-1 %
for L=1:m(i)
temp=0;
for j=L:m(i)
temp=temp+(-1)^(j-L) *C_0(i,j)*A(m(n),j-L+1)*B(i,n)^j;
end
C_1(i,L)=C_1(i,L)+B(i,n)^(m(n)-L)*temp;
end
end

%
for L=1:m(n)
for i=1:n-1
temp=0;
for j=1:m(i)
temp=temp+C_0(i,j)*A(j,m(n)-L+1)*B(n,i)^j;
end
C_1(n,L)=C_1(n,L)+B(i,n)^(m(n)-L)*temp;
end
end
C_0=C_1;
end

%__________________________________________
%___________________________________________
%Calculating coefficients of R(s)

N=length(D)-1;    %degree of numerator
K=sum(m);         %degree of denumerator
E=zeros(N-K+1,1); %coefficients of
if N>=K
E(1)=D(K+1);
e_n=E;
e_n(1)=1;
end
C_1=zeros(M+2,m_max+2);  % c_n
res=D(1)*C_0;
for n=2:N+1
%
%compute residues of R_n; store it in C_1;add it to res
for i=1:M
for j=1:m(i)
C_1(i,j)=(S(i)-s0)*C_0(i,j)+C_0(i,j+1);
res(i,j)=res(i,j)+D(n)*C_1(i,j);
end
end
%_______________________
%computing coefficients of residue polynomial of R_n; store it in
%e_n;add it to E_n
if n>K+1
for j=(n-K):-1:2
e_n(j)=e_n(j-1);
E(j)=E(j)+e_n(j)*D(n);
end
e_n(1)=0;
for i=1:M
for j=1:m(i)
if(s0==S(i))
continue;
end
e_n(1)=e_n(1)+C_1(i,j)/(s0-S(i))^j;
end
end
e_n(1)=-e_n(1);
E(1)=E(1)+e_n(1)*D(n);
end
C_0=C_1;
end
res(1:M,m_max+2)=S;
res(M+2,1:m_max)=1:m_max;

```