<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<div class="moz-cite-prefix">Hello Federico,</div>
<div class="moz-cite-prefix"><br>
</div>
<div class="moz-cite-prefix">Le 16/07/2021 à 21:11, Federico Miyara
a écrit :<br>
</div>
<blockquote type="cite"
cite="mid:286c21c3-10b7-9610-0314-83770c7ecd07@fceia.unr.edu.ar">
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
<br>
<font face="Courier New">Dear All,<br>
<br>
I'm wonderig whether there is a standard method to obtain the
partial fraction expansion of a rational function.<br>
<br>
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:<br>
<br>
// Poles and zeros<br>
p = [-8:-1];<br>
z = [0 0 0 0];<br>
<br>
// Numerator and denominator<br>
N = prod(%s - z)<br>
D = prod(%s - p)<br>
<br>
// Rational function<br>
H = N/D<br>
<br>
// Partial fractions<br>
HH = pfss(H)<br>
<br>
The result is <br>
<br>
HH = <br>
HH(1)<br>
<br>
-8 -5.7428571s -3.9515873s² +0.2960317s³ -0.0946429s⁴
-0.006746s⁵ -0.0001984s⁶ <br>
------------------------------------------------------------------------------
<br>
40320 +69264s +48860s² +18424s³ +4025s⁴ +511s⁵
+35s⁶ +s⁷ <br>
<br>
HH(2)<br>
<br>
0.0001984 <br>
--------- <br>
1 +s <br>
<br>
The second component in the list is indeed one of the partial
fractions, but the first one is not an irreductible version.
</font></blockquote>
<p><br>
The pfss() rmax option must be used, and set to a big value, for
instance 1/%eps. <br>
<a class="moz-txt-link-freetext" href="https://help.scilab.org/docs/6.1.0/en_US/pfss.html">https://help.scilab.org/docs/6.1.0/en_US/pfss.html</a><br>
Then:<br>
<font face="monospace"><br>
--> D = list2vec(pfss(H,1/%eps))'<br>
D = <br>
<br>
-0.8126984 3.3347222 -5.4 4.3402778 -1.7777778 0.3375
-0.0222222 0.0001984 <br>
---------- --------- ---- --------- ---------- ------
---------- --------- <br>
8 +s 7 +s 6 +s 5 +s 4 +s 3
+s 2 +s 1 +s <br>
<br>
--> clean(sum(D)-H)<br>
ans =<br>
<br>
0 <br>
- <br>
1 <br>
</font><br>
</p>
<blockquote type="cite"
cite="mid:286c21c3-10b7-9610-0314-83770c7ecd07@fceia.unr.edu.ar"><font
face="Courier New"> Computing HH(1) + HH(2) yields H with good
approximation.<br>
<br>
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 <br>
</font></blockquote>
<p><br>
In addition, partfrac() is available on FileExchange since 2015 @
<a class="moz-txt-link-freetext" href="https://fileexchange.scilab.org/toolboxes/451000">https://fileexchange.scilab.org/toolboxes/451000</a><br>
It manages poles multiplicity. With it, your example (with only
simple poles) gives the following:<br>
<br>
<font face="monospace">--> H = N/D<br>
H = <br>
s⁴ <br>
-------------------------------------------------------------------
<br>
40320 +109584s +118124s² +67284s³ +22449s⁴ +4536s⁵ +546s⁶
+36s⁷ +s⁸ <br>
<br>
--> [I, R, F, T] = partfrac(H);<br>
<br>
--> I // integer part of the fraction<br>
I = <br>
0.<br>
<br>
--> R // Undecomposed remainder of the fraction<br>
R = <br>
s⁴ <br>
-------------------------------------------------------------------
<br>
40320 +109584s +118124s² +67284s³ +22449s⁴ +4536s⁵ +546s⁶
+36s⁷ +s⁸ <br>
<br>
--> F // Decomposition: row #1 = numerators ; row #2 =
poles value ; row #3 = poles multiplicity<br>
F = <br>
<br>
-0.8126984 3.3347222 -5.4 4.3402778 -1.7777778 0.3375
-0.0222222 0.0001984<br>
-8. -7. -6. -5. -4. -3.
-2. -1. <br>
1. 1. 1. 1. 1.
1. 1. 1. <br>
<br>
--> T // Decomposition output as text<br>
T = <br>
" -0.8126984 3.3347222 -5.4 4.3402778 -1.7777778
0.3375 -0.0222222 0.0001984"<br>
"0 + ---------- + --------- + ---- + --------- + ---------- +
------ + ---------- + ---------"<br>
" 8+%s 7+%s 6+%s 5+%s 4+%s
3+%s 2+%s 1+%s "</font></p>
<p><font face="monospace"><br>
</font></p>
<p><font face="monospace">--> D = F(1,:)./((%s-F(2,:)).^F(3,:))<br>
D = <br>
</font></p>
<p><font face="monospace"> -0.8126984 3.3347222 -5.4 4.3402778
-1.7777778 0.3375 -0.0222222 0.0001984 <br>
---------- --------- ---- --------- ---------- ------
---------- --------- <br>
8 +s 7 +s 6 +s 5 +s 4 +s 3
+s 2 +s 1 +s <br>
<br>
--> clean(R-sum(D))<br>
ans =<br>
<br>
0 <br>
- <br>
1 </font></p>
Regards<br>
Samuel<br>
<br>
</body>
</html>