10
$\begingroup$

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}
]
$\endgroup$
  • 1
    $\begingroup$ @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. $\endgroup$ – jacobi16 Jul 15 at 11:59
6
$\begingroup$

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}}            *)
$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – jacobi16 Jul 15 at 17:16
  • $\begingroup$ 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 = [ ] $\endgroup$ – jacobi16 Jul 15 at 17:17
  • 1
    $\begingroup$ 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? $\endgroup$ – Roman Jul 15 at 17:22
  • 2
    $\begingroup$ 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? $\endgroup$ – Roman Jul 15 at 17:31
  • $\begingroup$ 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). $\endgroup$ – jacobi16 Jul 15 at 17:59
8
$\begingroup$

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.

$\endgroup$
5
$\begingroup$

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];

Mathematica graphics

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];

Mathematica graphics

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];

Mathematica graphics

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
$\endgroup$
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.