Mathematica function equivalent to MATLAB's residue function (partial fraction expansion)
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.
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}} *)
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