<
Partial fraction decomposition

To build the definite integral or to do a Laplace transformation (quite a topic in Electronic science) of a truly broken rational function with real coefficients like:


PartialFract


where m >= n, it’s quite handy to decompose this function into its partial fractions like:



PartialFract


With s1 … sk the roots of the denominator polynomial and A1 … Ak as single real values. That’s an interesting task as s1 … sk can be real or complex :-)

When I started to occupy myself with this task, I found some descriptions containing quite extensive operations with different cases for single or multiple real or complex roots. While studying these approaches and implementing some small algorithms to perform these operations I suddenly realized that all these operations are not necessary. In a computer it’s enough to stuff all the parts into a matrix equation and solve that. That yields all the needed values at once. The only crux is to stuff the things the right way (and of course first to determine the roots of the denominator polynomial. But to this we get later). To find how to stuff the right tings we first have to have a look into the basic approach and there are 4 cases:



Single real roots

Single real roots means all occurring roots of the denominator polynomial are real (not complex) and have different values:


PartialFract


And s1 … sk, the roots, are all real and have different values. In this case we have to find A1...Ak of the expression:


PartialFract


To solve this, the approach is to eliminate the denominator, to select some useful values for s and calculate the A values with these. In this case useful values for s are the roots of the denominator polynomial. With these values for each root all the Ax fractions become 0 except the one containing the root in its denominator.

The equation without the denominatro is:


PartialFract


And if we set s = s1 we get:


PartialFractPartialFract


With s = s2 we get:


PartialFractPartialFract


And so on.


This is basically the same as setting up a matrix equation with the upper formulation and s1..sk. The same matrix equation could be set up with just some random values for s as well. With s = s1..sk.in the main matrix all elements except the main diagonal elements are 0. With just some random s-values that would not be the case. That’s just to be mentioned :-)




Multiple real roots

Multiple real roots means the same real root occurs more than once and we have


PartialFract


A multiple root of the order n is decomposed into n fractions of the order n, n-1, n-2,…,1 like:


PartialFract


To handle this we first have a look at the case where there is only one triple root and nothing else:


PartialFract


This root is decomposed into 3 fractions and the partial fraction decomposition is carried out as:


PartialFract


The classic approach to solve this is again to eliminate the denominator polynomial. For the above sample that would look like:


PartialFractPartialFract



But now we cannot just set s = sk. That would not work. Here the approach is in, a first step, to differentiate the equation (in this case 2 times) to get rid of A2 and A3. With this A1 can be calculated. In a second step with A1 inserted into the equation and it is differentiated once. Then A2 can be determined and with A1 and A2 inserted A3 can be calculated. That’s quite easy for a triple root. But do that once for a sextuple root or higher :-)

If a computer should do this partial fraction decomposition, there is no need for any differentiation. Using the above formulation we can set up a matrix equation for this situation to. Similar to what I did for the first case. Only here we are not free in the selection of the values for s. The root should be avoided as value for s here. Otherwise the matrix equation could not be solved. 




Single complex root

As the coefficients of our broken rational function are real, complex roots will occur as conjugate complex pairs only. These conjugate complex pairs can be added like


PartialFract


And the sum of them is real.

And with this a partial fraction is built like:


PartialFract



The solution is the same as above: Eliminate the denominator polynomial. But then build a matrix equation with one row and column for each Ci and Cjs at once.



Multiple complex roots


For multiple complex roots the procedure is the same as for multiple real roots. Only the denominators are conjugate complex and the numerators are polynomials like.


PartialFractPartialFract


For a triple complex root.


The solution is the same again: Eliminate the denominator polynomial and build a matrix equation with one row and column for each Ci and each Cjs.

Now there is one more complication for setting up the matrix equation. If there are some real roots beside multiple complex roots, we cannot use the values of the real roots as s values for equation rows of the matrix equation. Otherwise the columns of the complex roots in the matrix equation become 0. Therefore the values of all real roots should be collected in a list and these values should be avoided while setting up the matrix equation.




If all these 4 cases are combined we have to build one big matrix equation containing all fractions.



The algorithm

When I started with this algorithm, my first thoughts about the design were, how to manage to have to implement as few decisions as possible in a sequence.

This is quite a simple task for single real roots. We just have to multiply the enumerator of each fraction by the denominator of all other fractions and remove the enumerator on the left side of the equation to set up the equation row for the matrix. For multiple real (or complex) roots things become a bit more complicate as we cannot just multiply an enumerator by all denominators. In the above sample with the tripple real root the matrix equation row looks like:


PartialFractPartialFract


The left side of the equation must be extended by the multiple root only once and A2 and A3 in the above sample are extended by (s – sk) and (s – sk)2. That requires some more decisions in the algorithm.

My idea to handle this was to implement a class for one fraction containing a denominator and a enumerator polynomial and to introduce 2 additional polynomials: An expander and a multiplier. The expander I use to expand the fractions in order to eliminate the denominator polynomial of the origin equation. A multiple root I like:


PartialFract


I virtually modify a little and threat it as follows


PartialFract


In this 3 fractions the first one is the parent of the other two’s. With this parent set in the other fractions I know that I don’t have to expand these one's by the expander of the parent but just multiply them by their multiplier polynomial. All the other fraction I just expand by all expander polynomials and that’s only once (s - sk)3 as the other 2 expanders are 1.

 For the first occurrence of a root I just copy the polynomial of this root into the expander and the denominator of this fraction. If I find a second root of the same value I first copy the expander of the first fraction into the multiplier of the new fraction and then multiply the expander of the first fraction by the polynomial of the second root and set the parent of the second root to the index of the first one. The expander of the second fraction I leave as 1. Then I increase the power variable of the parent fraction. If I find a third root of the same value I do that same procedure again and get the above formulation.

Exactly the same I do for single complex roots and multiple complex roots and with this I can fill my matrix equation.

When all the fractions are determined I set the power variables of all the multiple roots by walking through all elements and checking whether the power of the actual element is >1. If this is the case, I know that the following elements with this element as parent are child elements of this first parent element and set the power of these elements.

The denominator and the enumerator polynomial and the power finally are used for displaying the result only.

To compute the roots of the denominator polynomial I use my algorithm that computes the roots by the use of eigenvalues computation (see Computing roots of a polynomial by the use of eigenvalues for a detailed description) . It performs a QR decomposition and so the instance name of it is "qr". O.k. the QR algorithm is quite a complex algorithm but it yields all real and complex roots at once.



My partial fraction class with all its member variables looks like:


public class PartialFraction
{
    public double[] expander;
    public double[] enumerator;
    public double[] denominator;
    public POLETYPE type;
    public int parent;
    public int power;
    public double[] multiplier;
 
public PartialFraction(int size)
{
expander = new double[size];    // polynomial to expand each other polynomial
denominator = new double[size]; // denominator polynomial
enumerator = new double[1];     // enumerator polynomial
enumerator[0] = 1;
multiplier = new double[1];     // polynomial to multiply multiple and complex poles
multiplier[0] = 1;
type = POLETYPE.Single;
parent = -1;
power = 1;
}
 
public void InitExpander(double[] source)
{
CopyPoly(source, ref expander);
}
 
public void InitDenominator(double[] source)
{
CopyPoly(source, ref denominator);
}
 
private void CopyPoly(double[] source, ref double[] dest)
{
int i;
Array.Resize(ref dest, source.Length);
for (i = 0; i < dest.Length; i++)
{
dest[i] = source[i];
}
}
}




And the loop to fill the roots into the data sequence is



foreach (double[] root in qr.roots)
{
found = false;
if (i > 0)
{
for (k = 0; k < i; k++)
{
double[] tempRoot = new double[2];
CopyPoly(root, ref tempRoot);
if (AreEqual(partialFract[k].denominator, tempRoot) && !found)
{
// already one root of the same value found -> multiple root
if (!found)
{
// allready one pole found => multiple real pole
partialFract[i] = new PartialFraction(1);
partialFract[i].expander[0] = 1.0;
partialFract[i].InitDenominator(root);
partialFract[i].multiplier = MultPoly(partialFract[i].multiplier, partialFract[k].expander);
partialFract[i].parent = k;
partialFract[i].type = POLETYPE.Multiple;
partialFract[k].type = POLETYPE.Multiple;
 
if (root.Length > 2) // complex root
{
partialFract[i].multiplier[0] = 1.0;
Array.Resize(ref partialFract, partialFract.Length + 1);
i++;
 
partialFract[i] = new PartialFraction(1);
partialFract[i].expander[0] = 1.0;
partialFract[i].InitDenominator(root);
Array.Resize(ref partialFract[i].enumerator, 2);
partialFract[i].enumerator[1] = 0;
Array.Resize(ref partialFract[i].multiplier, 2);
partialFract[i].multiplier[0] = 1.0;
partialFract[i].multiplier[1] = 0;
partialFract[i].multiplier = MultPoly(partialFract[i].multiplier, partialFract[k].expander);
partialFract[i].type = POLETYPE.Multiple;
partialFract[i].parent = k;
partialFract[k + 1].power++;
}
else
{
multipleRoots.Add(Math.Round(root[1], 2));
}
partialFract[k].expander = MultPoly(partialFract[k].expander, tempRoot);
partialFract[k].power++;
found = true;
i++;
}
}
}
}
if (!found)
{
partialFract[i] = new PartialFraction(root.Length);
partialFract[i].InitExpander(root);
partialFract[i].InitDenominator(root);
if (root.Length > 2) // complex root
{
// conjugate complex pole pair
partialFract[i].type = POLETYPE.Complex;
partialFract[i].multiplier[0] = 1.0;
Array.Resize(ref partialFract, partialFract.Length + 1);
i++;
 
partialFract[i] = new PartialFraction(1);
partialFract[i].expander[0] = 1.0;
partialFract[i].InitDenominator(root);
Array.Resize(ref partialFract[i].enumerator, 2);
partialFract[i].enumerator[1] = 0;
Array.Resize(ref partialFract[i].multiplier, 2);
partialFract[i].multiplier[0] = 1.0;
partialFract[i].multiplier[1] = 0;
partialFract[i].type = POLETYPE.Complex;
partialFract[i].parent = i - 1;
}
else
{
// real pol
partialFract[i].type = POLETYPE.Single;
}
i++;
}
}



Now there is the value of the power for the child roots (the ones that have a parent > -1). These values are not set jet. Therefor I implemented the function


private void SetPowerValues()
{
     int i, j;
     SPower[] pwrValues = new SPower[0];
     for (i = 0; i < partialFract.Length; i++)
     {
         if (partialFract[i].power > 1)
         {
              Array.Resize(ref pwrValues, pwrValues.Length + 1);
              pwrValues[pwrValues.Length-1].index = i;
              pwrValues[pwrValues.Length - 1].power = partialFract[i].power;
         }
         else
         {
              if (pwrValues.Length >= 1)
              {
                   for(j = 0; j < pwrValues.Length; j++)
                   {
                        if (pwrValues[j].index == partialFract[i].parent)
                        {
                            pwrValues[j].power--;
                            if (partialFract[i].denominator.Length > 2)
                            {
                                 partialFract[i].power = pwrValues[j].power;
                                 partialFract[i + 1].power = pwrValues[j].power;
                                 i++;
                            }
                            else
                                 partialFract[i].power = pwrValues[j].power;
                        }
                   }
              }
         }
     }
}



That sets all the power values.


As input I get the found roots in qr.roots. This is where my Eigenvalue algorithm puts its results and qr is the instance of this QR algorithm.



To fill the matrix equation and solve it I implemented the function


private void BuildMatrixAndSolve(int numberOfFractions)


With the sequence


for (k = 0; k < numberOfFractions; k++)
{
for (i = 0; i < numberOfFractions; i++)
{
nom = 0;
Array.Resize(ref tempDenom, actPosition + 1);
tempDenom[0] = 1;
for (j = 0; j < numberOfFractions; j++)
{
if ((j != i) && !((partialFract[i].type > POLETYPE.Single) && (j == partialFract[i].parent)))
tempDenom = MultPoly(tempDenom, partialFract[j].expander);
}
if (partialFract[i].type > POLETYPE.Single)
{
tempDenom = MultPoly(tempDenom, partialFract[i].multiplier);
}
for (j = 0; j < tempDenom.Length; j++)
{
nom = nom + tempDenom[j] * Math.Pow(equationValue[k], tempDenom.Length - j - 1);
}
cell = dGMatrix[i, k];
cell.Value = Convert.ToString(Convert.ToString(nom));
m.a[k, i] = nom;
}
res = 0;
for (j = 0; j < enumeratorPolynom.Length; j++)
{
res = res + enumeratorPolynom[j] * Math.Pow(equationValue[k], enumeratorPolynom.Length - j - 1);
}
cell = dgResult[0, k];
cell.Value = Convert.ToString(Convert.ToString(res));
m.y[k] = res;
}
gauss.Eliminate();
gauss.Solve();



I use the list


public List<double> realRoots = new List<double>();



to store the values of all real roots in order to avoid these values while filling the matrix equation.

Therefore I determine the values to calculate the matrix rows in the loop


for (i = 0; i < orderDenom; i++)
{
while (realRoots.Contains(s))
s++;
equationValue[i] = s;
s++;
}



And use the array


double[] equationValue = new double[orderDenom];



to keep the values that can be used for the matrix equation.


With this approach I ended with a quite manageable algorithm (at least I think so :-) )



Sample project


In my sample project there is one main window


PartialFract



There the top data grid contains the broken rational function with its enumerator in the top row and the denominator in the bottom row. In the combobox on the right side the order of the denominator polynomial can be set. If the enumerator shall have a lower order, just leave the last fields empty. In the middle the set up matrix equation is displayed and in the bottom data grid the fractions will be displayed.

To use it, first set the denominator order. Then enter the denominator and the enumerator with the highest power starting on the left side. Pressing <calc> computes the roots, fills the matrix equation, solves it and displays the result in the bottom data grid.


On start-up the window contains the sample


PartialFract


which has one single real pole at


PartialFract


And a double complex at


PartialFract



C# Demo Project Partial fraction decomposition
  • PartialFraction.zip