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 = {an, $$\dots$$, a2, a1}

coefficients of the polynomial b = {bm, $$\dots$$, b2, b1}

Three lists as outputs:

residues r = {rn, $$\dots$$ , r2, r1}

poles p = {pn, $$\dots$$ , p2, p1}

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 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 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 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 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 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 at 17:59

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


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 review queue and will at some point be available as ResourceFunction["IncompletePFD"]. It does (or at least tries to do) what the name implies, the so-called "incomplete" partial fraction decomposition.