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;