5
$\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 = {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}
]
$\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
5
$\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
  • 1
    $\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
4
$\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$
4
$\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 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.

$\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.