# Mathematica function equivalent to MATLAB's residue function (partial fraction expansion)

I am looking for a Mathematica function equivalent to MATLAB's residue function.

If there is no Mathematica equivalent, I would like to write a function that given the coefficients of two polynomials $$a$$ and $$b$$, outputs the residues $$r$$, the poles $$p$$ and the term $$k$$ of the partial fraction expansion of $$\frac{b}{a}$$:

$$\frac{b(x)}{a(x)} = \frac{b_m x^m + b_{m-1}x^{m-1} + \dots + b_1 x + b_0}{a_n x^n + a_{n-1}x^{n-1} + \dots + a_1 x + a_0} = \frac{r_n}{x-p_n} + \dots + \frac{r_2}{x-p_2} + \frac{r_1}{x-p_1} + k(x) \tag{1}$$

Two lists as inputs:

coefficients of the polynomial $$a = \{a_n, \dots, a_2, a_1\}$$

coefficients of the polynomial $$b = \{b_m, \dots, b_2, b_1\}$$

Three lists as outputs:

residues $$r = \{r_n, \dots, r_2, r_1\}$$

poles $$p = \{p_n, \dots, p_2, p_1\}$$

polynomial $$k = \{\dots\}$$

I know how to create the polynomials $$a(x)$$ and $$b(x)$$ from the input lists $$a$$ and $$b$$ and how to do a partial fraction expansion. However, I haven't been able to figure out how to compute $$r$$, $$p$$ and $$k$$

Residues[a_,b_] :=
Module[{apoly, bpoly,f,x},
(* build polynomials *)
bpoly = FromDigits[b, x];
apoly = FromDigits[a, x];
(* partial fraction expansion *)
f = Apart[Expand[bpoly]/Expand[apoly]];

r = ???;
p = ???;
k = ???;

{r,p,k}
]

• @HenrikSchumacher I do not want to simply do a partial fraction expansion (which can be done using Apart) but rather extract the residues, poles and constant term from this expansion. – jacobi16 Jul 15 '19 at 11:59

You get the poles with

p = x /. Solve[apoly == 0, x]


and the residues with

r = RootReduce[Residue[f, {x, #}]] & /@ p


I don't know how to get k directly though, except to do a difference (very inefficient):

k = f - Total[r/(x - p)] // Together // FullSimplify


All together in one function:

residue[num_, denom_] := Module[{apoly, bpoly, f, p, r, k},
bpoly = FromDigits[num, x];
apoly = FromDigits[denom, x];
f = bpoly/apoly;
p = x /. Solve[apoly == 0, x];
r = RootReduce[Residue[f, {x, #}]] & /@ p;
k = f - Total[r/(x - p)] // Together // FullSimplify;
{r, p, CoefficientList[k, x]}]


Let's go through Nasser's examples:

residue[{2, 1, 0, 0}, {1, 0, 1, 1}]
(*    {{-0.0708358, 0.535418 - 1.03899 I, 0.535418 + 1.03899 I},
{-0.682328, 0.341164 - 1.16154 I, 0.341164 + 1.16154 I},
{2}}                                                         *)

residue[{-4, 8}, {1, 6, 8}]
(*    {{-12, 8},
{-4, -2},
{}}          *)

residue[{2, 0, 0, 1, 0}, {1, 0, 1}]
(*    {{1/2 + I, 1/2 - I},
{-I, I},
{-2, 0, 2}}            *)

• When I run residue[{1},{-3.,2.,0.,0.,0.}] I get a message CoefficientList::poly: 0. +0.5/x^3+0.75/x^2-2.25/x is not a polynomial. and incorrect results for r and k – jacobi16 Jul 15 '19 at 17:16
• For the same example in Matlab, [r,p,k]=residue([1],[-3,2,0,0,0]) I obtain r = [-1.1250 1.1250 0.7500 0.5000], p = [0.6667 0 0 0], k = [ ] – jacobi16 Jul 15 '19 at 17:17
• Yes, there's a problem in this code for higher-order roots (in this case, $x=0$ is a third-order root). The problem lies in the calculation of $k$, it needs to be adapted and I don't know how. Maybe ask at math.SE for a robust algorithm? – Roman Jul 15 '19 at 17:22
• I don't think the Matlab answer is a correct decomposition: $-\frac98/(x - \frac23) + \frac98/x + \frac34/x + \frac12/x$ is not the same as your input polynomial $1/(2 x^3 - 3 x^4)$. What is the correct answer for this case? – Roman Jul 15 '19 at 17:31
• You are indeed correct. The Matlab decomposition is wrong. (sorry for not checking, I assumed Matlab was doing it properly). Matlab's documentation states Numerically,the partial fraction expansion of a ratio of polynomials represents an ill-posed problem.If the denominator polynomial,a(x), is near a polynomial with multiple roots,then small changes in the data,including roundoff errors,can result in arbitrarily large changes in the resulting poles and residues. I guess residue([1],[-3,2,0,0,0]) is such a case, i.e.,$\frac{1}{-3x^4 + 2x^3}$ cannot be written in the form of eq. (1). – jacobi16 Jul 15 '19 at 17:59

This is cheating, but:

num = FromDigits[{2, 1, 0, 0}, x]
den = FromDigits[{1., 0, 1, 1}, x]

(* Out[7]= x^2 (1 + 2 x)

Out[8]= 1 + x + 1. x^3 *)

In[10]:= IncompletePFD[num, den]

(* Out[10]= (2. + 0. I) + (
0.535417905939 + 1.03899170871 I)/((-0.341163901914 - 1.1615414 I) +
x) + (0.535417905939 -
1.03899170871 I)/((-0.341163901914 + 1.1615414 I) + x) - (
0.0708358118772 + 0. I)/(0.682327803828 + x) *)


IncompletePFD is in the Wolfram Function Repository and is available as ResourceFunction["IncompletePFD"]. It does (or at least tries to do) what the name implies, the so-called "incomplete" partial fraction decomposition.

I just saw the other answer uses some build in function to find the residues which I did not know about and which might be more efficient that what I did below. So this is something that can be improved in the function below if needed.

residue[numer0_, denom0_] :=
Module[{x, result, k, roots, numer, denom, factors, n, p1, p2, p3},
numer = FromDigits[numer0, x];
denom = FromDigits[denom0, x];
If[Length[numer0] >= Length[denom0],
k = PolynomialQuotient[numer, denom, x];
k = CoefficientList[k, x],
k = {}
];
roots = x /. NSolve[denom == 0, x];
factors = Flatten@Last@Reap@Do[
If[n == 1,
p1 = roots[[2 ;; -1]],
p1 = roots[[Join[Range[1, n - 1], Range[n + 1, Length[roots]]]]]
];
p2 = (Times @@ ((x - #) & /@ p1)) /. x -> roots[[n]];
p3 = numer /. x -> roots[[n]];
Sow[Simplify[p3/p2]]
,
{n, 1, Length[roots]}
];
{Reverse@Chop@factors, Reverse@Chop@roots, k}
]


Example 1

numer = {2, 1, 0, 0};
denom = {1, 0, 1, 1};
{r, p, k} = residue[numer, denom];


Matlab:

>> numerator = [2 1 0 0];
>> denominator = [1 0 1 1];
>> [r,p,k] = residue(numerator,denominator)

r =
0.5354 + 1.0390i
0.5354 - 1.0390i
-0.0708 + 0.0000i
p =
0.3412 + 1.1615i
0.3412 - 1.1615i
-0.6823 + 0.0000i

k =
2


Example 2

numer = {-4, 8};
denom = {1, 6, 8};
{r, p, k} = residue[numer, denom];


Matlab:

>> numerator = [-4 8];
denomenator = [1 6 8];
[r,p,k] = residue(numerator ,denomenator)

r =
-12
8
p =
-4
-2
k =
[]


example 3

numer = {2, 0, 0, 1, 0};
denom = {1, 0, 1};
{r, p, k} = residue[numer, denom];


Matlab

>> b = [2 0 0 1 0];
a = [1 0 1];
[r,p,k] = residue(b,a)

r =
0.5000 - 1.0000i
0.5000 + 1.0000i

p =
0.0000 + 1.0000i
0.0000 - 1.0000i

k =
2     0    -2


It has been mentioned in other answers that the case of the denominator having multiple roots is problematic. One should then recall what Residue[] actually does, and instead use the more general functionality of Series[] to get the necessary components of the partial fraction decomposition.

Let us demonstrate on two of Nasser's examples. Here's an easy one with complex residues:

rat = FromDigits[{2, 0, 0, 1, 0}, x]/FromDigits[{1, 0, 1}, x];

{num, den} = {Numerator[rat], Denominator[rat]};
{pp, rem} = PolynomialQuotientRemainder[num, den, x];

pfd = pp + Total[Normal[FullSimplify[Series[rem/den, {x, #, -1}]] & /@
DeleteDuplicates[x /. Solve[den == 0, x]]]]
-2 + 2 x^2 + (1/2 - I)/(-I + x) + (1/2 + I)/(I + x)


and a slightly harder case with an irreducible denominator:

rat = FromDigits[{2, 1, 0, 0}, x]/FromDigits[{1, 0, 1, 1}, x]

{num, den} = {Numerator[rat], Denominator[rat]};
{pp, rem} = PolynomialQuotientRemainder[num, den, x];

pfd = pp + Total[Normal[FullSimplify[Series[rem/den, {x, #, -1}]] & /@
DeleteDuplicates[x /. Solve[den == 0, x,
Cubics -> False]]]]
2 + Root[3 + 40 #1 - 31 #1^2 + 31 #1^3 &, 1]/(x - Root[1 + #1 + #1^3 &, 1]) +
Root[3 + 40 #1 - 31 #1^2 + 31 #1^3 &, 2]/(x - Root[1 + #1 + #1^3 &, 2]) +
Root[3 + 40 #1 - 31 #1^2 + 31 #1^3 &, 3]/(x - Root[1 + #1 + #1^3 &, 3])


where Cubics -> False was necessary to prevent needless complexity.

Now for examples with multiple roots. Here is a relatively simple one:

rat = (37 - 72 x + 20 x^2 + 25 x^3 + x^4 + x^5 - 7 x^6 - x^7 + x^8)/
(12 - 28 x + 11 x^2 + 14 x^3 - 8 x^4 - 2 x^5 + x^6);

{num, den} = {Numerator[rat], Denominator[rat]};
{pp, rem} = PolynomialQuotientRemainder[num, den, x];

pfd = pp + Total[Normal[FullSimplify[Series[rem/den, {x, #, -1}]] & /@
DeleteDuplicates[x /. Solve[den == 0, x]]]]
3 + 271/(200 (-3 + x)) - 5/(18 (-1 + x)^3) - 61/(108 (-1 + x)^2) - 149/(216 (-1 + x)) +
x + x^2 - 19/(135 (2 + x)^2) + 226/(675 (2 + x))


Check:

Simplify[Together[pfd - rat]]
0


Here is a more difficult case:

rat = (-17 - 39 x - 26 x^2 + 30 x^3 + 45 x^4 + 15 x^5 - 17 x^6 - 9 x^7 -
4 x^8 + 2 x^9 + x^10)/
(-6 - 11 x - 3 x^2 + 15 x^3 + 11 x^4 - 4 x^5 - 8 x^6 + x^7 + x^8);


Perform the procedure once more:

{num, den} = {Numerator[rat], Denominator[rat]};
{pp, rem} = PolynomialQuotientRemainder[num, den, x];

pfd = pp + Total[Normal[FullSimplify[Series[rem/den, {x, #, -1}]] & /@
DeleteDuplicates[x /. Solve[den == 0, x]]]];


This time, I have mercifully suppressed the output. Nevertheless, here is how to check the result:

Collect[Numerator[Together[pfd - rat]], x, RootReduce]
0


Note that this procedure is intended for polynomials with exact coefficients; involving inexact numbers is a much thornier problem, involving the decision of when to take a cluster of roots as a single multiple root. Encapsulation of this method into a routine is left as an exercise.