[Scilab-users] partial fractions

[Samuel Gougeon]

Hello Federico,

Le 16/07/2021 à 21:11, Federico Miyara a écrit :
>
> Dear All,
>
> I'm wonderig whether there is a standard method to obtain the partial 
> fraction expansion of a rational function.
>
> Searching, I've found the function pfss(), which purportedly is meant 
> for that, but I find it doesn't give a complete solution. For instance:
>
> // Poles and zeros
> p = [-8:-1];
> z = [0 0 0 0];
>
> // Numerator and denominator
> N = prod(%s - z)
> D = prod(%s - p)
>
> // Rational function
> H = N/D
>
> // Partial fractions
> HH = pfss(H)
>
> The result is
>
> HH  =
>        HH(1)
>
>    -8 -5.7428571s -3.9515873s² +0.2960317s³ -0.0946429s⁴ -0.006746s⁵ 
> -0.0001984s⁶
> ------------------------------------------------------------------------------ 
>
>               40320 +69264s +48860s² +18424s³ +4025s⁴ +511s⁵ +35s⁶ +s⁷
>
>        HH(2)
>
>    0.0001984
>    ---------
>      1 +s
>
> The second component in the list is indeed one of the partial 
> fractions, but the first one is not an irreductible version. 


The pfss() rmax option must be used, and set to a big value, for 
instance 1/%eps.
https://help.scilab.org/docs/6.1.0/en_US/pfss.html
Then:

--> D = list2vec(pfss(H,1/%eps))'
  D  =

    -0.8126984  3.3347222  -5.4  4.3402778  -1.7777778  0.3375 
-0.0222222  0.0001984
    ----------  ---------  ----  ---------  ----------  ------ 
----------  ---------
       8 +s       7 +s     6 +s    5 +s        4 +s      3 +s      2 
+s       1 +s

--> clean(sum(D)-H)
  ans  =

    0
    -
    1

> Computing HH(1) + HH(2) yields H with good approximation.
>
> I guess I could proceed iteratively with HH(1), but it would be nice 
> to have a function performing this automatically. I couldn't find such 
> a function


In addition, partfrac() is available on FileExchange since 2015 @ 
https://fileexchange.scilab.org/toolboxes/451000
It manages poles multiplicity. With it, your example (with only simple 
poles) gives the following:

--> H = N/D
  H  =
s⁴
-------------------------------------------------------------------
    40320 +109584s +118124s² +67284s³ +22449s⁴ +4536s⁵ +546s⁶ +36s⁷ +s⁸

--> [I, R, F, T] = partfrac(H);

--> I  // integer part of the fraction
  I  =
    0.

--> R   // Undecomposed remainder of the fraction
  R  =
s⁴
-------------------------------------------------------------------
    40320 +109584s +118124s² +67284s³ +22449s⁴ +4536s⁵ +546s⁶ +36s⁷ +s⁸

--> F   // Decomposition: row #1 = numerators ; row #2 = poles value ; 
row #3 = poles multiplicity
  F  =

   -0.8126984   3.3347222  -5.4   4.3402778  -1.7777778   0.3375 
-0.0222222   0.0001984
   -8.         -7.         -6.   -5.         -4.         -3. -2.         
-1.
    1.          1.          1.    1.          1. 1.       1.          1.

--> T  // Decomposition output as text
  T  =
   "    -0.8126984   3.3347222   -5.4   4.3402778   -1.7777778 0.3375   
-0.0222222   0.0001984"
   "0 + ---------- + --------- + ---- + --------- + ---------- + ------ 
+ ---------- + ---------"
   "       8+%s        7+%s      6+%s     5+%s         4+%s 3+%s       
2+%s        1+%s   "


--> D = F(1,:)./((%s-F(2,:)).^F(3,:))
  D  =

    -0.8126984  3.3347222  -5.4  4.3402778 -1.7777778  0.3375  
-0.0222222  0.0001984
    ----------  ---------  ----  ---------  ----------  ------ 
----------  ---------
       8 +s       7 +s     6 +s    5 +s        4 +s      3 +s      2 
+s       1 +s

-->  clean(R-sum(D))
  ans  =

    0
    -
    1

Regards
Samuel
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20211229/51f5d770/attachment.htm>