|
2847
|
1 ## Copyright (C) 1996, 1997 John W. Eaton |
|
2313
|
2 ## |
|
|
3 ## This file is part of Octave. |
|
|
4 ## |
|
|
5 ## Octave is free software; you can redistribute it and/or modify it |
|
|
6 ## under the terms of the GNU General Public License as published by |
|
|
7 ## the Free Software Foundation; either version 2, or (at your option) |
|
|
8 ## any later version. |
|
|
9 ## |
|
|
10 ## Octave is distributed in the hope that it will be useful, but |
|
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
13 ## General Public License for more details. |
|
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
|
16 ## along with Octave; see the file COPYING. If not, write to the Free |
|
|
17 ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA |
|
|
18 ## 02111-1307, USA. |
|
904
|
19 |
|
3368
|
20 ## -*- texinfo -*- |
|
|
21 ## @deftypefn {Function File} {} residue (@var{b}, @var{a}, @var{tol}) |
|
|
22 ## If @var{b} and @var{a} are vectors of polynomial coefficients, then |
|
|
23 ## residue calculates the partial fraction expansion corresponding to the |
|
|
24 ## ratio of the two polynomials. |
|
|
25 ## @cindex partial fraction expansion |
|
3426
|
26 ## |
|
3368
|
27 ## The function @code{residue} returns @var{r}, @var{p}, @var{k}, and |
|
|
28 ## @var{e}, where the vector @var{r} contains the residue terms, @var{p} |
|
|
29 ## contains the pole values, @var{k} contains the coefficients of a direct |
|
|
30 ## polynomial term (if it exists) and @var{e} is a vector containing the |
|
|
31 ## powers of the denominators in the partial fraction terms. |
|
3426
|
32 ## |
|
3368
|
33 ## Assuming @var{b} and @var{a} represent polynomials |
|
|
34 ## @iftex |
|
|
35 ## @tex |
|
|
36 ## $P(s)$ and $Q(s)$ |
|
|
37 ## @end tex |
|
|
38 ## @end iftex |
|
|
39 ## @ifinfo |
|
|
40 ## P (s) and Q(s) |
|
|
41 ## @end ifinfo |
|
|
42 ## we have: |
|
|
43 ## @iftex |
|
|
44 ## @tex |
|
|
45 ## $$ |
|
|
46 ## {P(s)\over Q(s)} = \sum_{m=1}^M {r_m\over (s-p_m)^e_m} |
|
|
47 ## + \sum_{i=1}^N k_i s^{N-i}. |
|
|
48 ## $$ |
|
|
49 ## @end tex |
|
|
50 ## @end iftex |
|
|
51 ## @ifinfo |
|
3426
|
52 ## |
|
3368
|
53 ## @example |
|
2311
|
54 ## P(s) M r(m) N |
|
3368
|
55 ## ---- = SUM ------------- + SUM k(i)*s^(N-i) |
|
|
56 ## Q(s) m=1 (s-p(m))^e(m) i=1 |
|
|
57 ## @end example |
|
|
58 ## @end ifinfo |
|
3426
|
59 ## |
|
3368
|
60 ## @noindent |
|
|
61 ## where @var{M} is the number of poles (the length of the @var{r}, |
|
|
62 ## @var{p}, and @var{e} vectors) and @var{N} is the length of the @var{k} |
|
|
63 ## vector. |
|
3426
|
64 ## |
|
3368
|
65 ## The argument @var{tol} is optional, and if not specified, a default |
|
|
66 ## value of 0.001 is assumed. The tolerance value is used to determine |
|
|
67 ## whether poles with small imaginary components are declared real. It is |
|
|
68 ## also used to determine if two poles are distinct. If the ratio of the |
|
|
69 ## imaginary part of a pole to the real part is less than @var{tol}, the |
|
|
70 ## imaginary part is discarded. If two poles are farther apart than |
|
|
71 ## @var{tol} they are distinct. For example, |
|
3426
|
72 ## |
|
3368
|
73 ## @example |
|
|
74 ## @group |
|
|
75 ## b = [1, 1, 1]; |
|
2311
|
76 ## a = [1, -5, 8, -4]; |
|
3368
|
77 ## [r, p, k, e] = residue (b, a); |
|
|
78 ## @result{} r = [-2, 7, 3] |
|
|
79 ## @result{} p = [2, 2, 1] |
|
|
80 ## @result{} k = [](0x0) |
|
|
81 ## @result{} e = [1, 2, 1] |
|
|
82 ## @end group |
|
|
83 ## @end example |
|
3426
|
84 ## |
|
3368
|
85 ## @noindent |
|
|
86 ## which implies the following partial fraction expansion |
|
|
87 ## @iftex |
|
|
88 ## @tex |
|
|
89 ## $$ |
|
|
90 ## {s^2+s+1\over s^3-5s^2+8s-4} = {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} |
|
|
91 ## $$ |
|
|
92 ## @end tex |
|
|
93 ## @end iftex |
|
|
94 ## @ifinfo |
|
3426
|
95 ## |
|
3368
|
96 ## @example |
|
|
97 ## s^2 + s + 1 -2 7 3 |
|
2311
|
98 ## ------------------- = ----- + ------- + ----- |
|
|
99 ## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1) |
|
3368
|
100 ## @end example |
|
|
101 ## @end ifinfo |
|
|
102 ## @end deftypefn |
|
|
103 |
|
2311
|
104 ## SEE ALSO: poly, roots, conv, deconv, polyval, polyderiv, polyinteg |
|
1025
|
105 |
|
3202
|
106 ## Author: Tony Richardson <arichard@stark.cc.oh.us> |
|
2312
|
107 ## Created: June 1994 |
|
|
108 ## Adapted-By: jwe |
|
559
|
109 |
|
2312
|
110 function [r, p, k, e] = residue (b, a, toler) |
|
559
|
111 |
|
2303
|
112 ## Here's the method used to find the residues. |
|
|
113 ## The partial fraction expansion can be written as: |
|
|
114 ## |
|
|
115 ## |
|
|
116 ## P(s) D M(k) A(k,m) |
|
|
117 ## ---- = # # ------------- |
|
|
118 ## Q(s) k=1 m=1 (s - pr(k))^m |
|
|
119 ## |
|
|
120 ## (# is used to represent a summation) where D is the number of |
|
|
121 ## distinct roots, pr(k) is the kth distinct root, M(k) is the |
|
|
122 ## multiplicity of the root, and A(k,m) is the residue cooresponding |
|
|
123 ## to the kth distinct root with multiplicity m. For example, |
|
|
124 ## |
|
|
125 ## s^2 A(1,1) A(2,1) A(2,2) |
|
|
126 ## ------------------- = ------ + ------ + ------- |
|
|
127 ## s^3 + 4s^2 + 5s + 2 (s+2) (s+1) (s+1)^2 |
|
|
128 ## |
|
|
129 ## In this case there are two distinct roots (D=2 and pr = [-2 -1]), |
|
|
130 ## the first root has multiplicity one and the second multiplicity |
|
|
131 ## two (M = [1 2]) The residues are actually stored in vector format as |
|
|
132 ## r = [ A(1,1) A(2,1) A(2,2) ]. |
|
|
133 ## |
|
|
134 ## We then multiply both sides by Q(s). Continuing the example: |
|
|
135 ## |
|
|
136 ## s^2 = r(1)*(s+1)^2 + r(2)*(s+1)*(s+2) + r(3)*(s+2) |
|
|
137 ## |
|
|
138 ## or |
|
|
139 ## |
|
|
140 ## s^2 = r(1)*(s^2+2s+1) + r(2)*(s^2+3s+2) +r(3)*(s+2) |
|
|
141 ## |
|
|
142 ## The coefficients of the polynomials on the right are stored in a row |
|
|
143 ## vector called rhs, while the coefficients of the polynomial on the |
|
|
144 ## left is stored in a row vector called lhs. If the multiplicity of |
|
|
145 ## any root is greater than one we'll also need derivatives of this |
|
|
146 ## equation of order up to the maximum multiplicity minus one. The |
|
|
147 ## derivative coefficients are stored in successive rows of lhs and |
|
|
148 ## rhs. |
|
|
149 ## |
|
|
150 ## For our example lhs and rhs would be: |
|
|
151 ## |
|
|
152 ## | 1 0 0 | |
|
|
153 ## lhs = | | |
|
|
154 ## | 0 2 0 | |
|
|
155 ## |
|
|
156 ## | 1 2 1 1 3 2 0 1 2 | |
|
|
157 ## rhs = | | |
|
|
158 ## | 0 2 2 0 2 3 0 0 1 | |
|
|
159 ## |
|
|
160 ## We then form a vector B and a matrix A obtained by evaluating the |
|
|
161 ## polynomials in lhs and rhs at the pole values. If a pole has a |
|
|
162 ## multiplicity greater than one we also evaluate the derivative |
|
|
163 ## polynomials (successive rows) at the pole value. |
|
|
164 ## |
|
|
165 ## For our example we would have |
|
|
166 ## |
|
|
167 ## | 4| | 1 0 0 | | r(1) | |
|
|
168 ## | 1| = | 0 0 1 | * | r(2) | |
|
|
169 ## |-2| | 0 1 1 | | r(3) | |
|
|
170 ## |
|
|
171 ## We then solve for the residues using matrix division. |
|
559
|
172 |
|
1025
|
173 if (nargin < 2 || nargin > 3) |
|
3456
|
174 usage ("residue (b, a, toler)"); |
|
559
|
175 endif |
|
|
176 |
|
|
177 if (nargin == 2) |
|
|
178 toler = .001; |
|
|
179 endif |
|
|
180 |
|
2303
|
181 ## Make sure both polynomials are in reduced form. |
|
1025
|
182 |
|
|
183 a = polyreduce (a); |
|
|
184 b = polyreduce (b); |
|
559
|
185 |
|
1025
|
186 b = b / a(1); |
|
|
187 a = a / a(1); |
|
559
|
188 |
|
1025
|
189 la = length (a); |
|
|
190 lb = length (b); |
|
559
|
191 |
|
2303
|
192 ## Handle special cases here. |
|
1025
|
193 |
|
|
194 if (la == 0 || lb == 0) |
|
559
|
195 k = r = p = e = []; |
|
|
196 return; |
|
|
197 elseif (la == 1) |
|
1025
|
198 k = b / a; |
|
559
|
199 r = p = e = []; |
|
|
200 return; |
|
|
201 endif |
|
|
202 |
|
2303
|
203 ## Find the poles. |
|
1025
|
204 |
|
|
205 p = roots (a); |
|
|
206 lp = length (p); |
|
559
|
207 |
|
2303
|
208 ## Determine if the poles are (effectively) real. |
|
1025
|
209 |
|
|
210 index = find (abs (imag (p) ./ real (p)) < toler); |
|
|
211 if (length (index) != 0) |
|
|
212 p (index) = real (p (index)); |
|
559
|
213 endif |
|
|
214 |
|
2303
|
215 ## Find the direct term if there is one. |
|
1025
|
216 |
|
|
217 if (lb >= la) |
|
2303
|
218 ## Also returns the reduced numerator. |
|
1025
|
219 [k, b] = deconv (b, a); |
|
|
220 lb = length (b); |
|
559
|
221 else |
|
|
222 k = []; |
|
|
223 endif |
|
|
224 |
|
1025
|
225 if (lp == 1) |
|
|
226 r = polyval (b, p); |
|
559
|
227 e = 1; |
|
|
228 return; |
|
|
229 endif |
|
|
230 |
|
|
231 |
|
2303
|
232 ## We need to determine the number and multiplicity of the roots. |
|
|
233 ## |
|
|
234 ## D is the number of distinct roots. |
|
|
235 ## M is a vector of length D containing the multiplicity of each root. |
|
|
236 ## pr is a vector of length D containing only the distinct roots. |
|
|
237 ## e is a vector of length lp which indicates the power in the partial |
|
|
238 ## fraction expansion of each term in p. |
|
559
|
239 |
|
2303
|
240 ## Set initial values. We'll remove elements from pr as we find |
|
|
241 ## multiplicities. We'll shorten M afterwards. |
|
1025
|
242 |
|
|
243 e = ones (lp, 1); |
|
|
244 M = zeros (lp, 1); |
|
559
|
245 pr = p; |
|
1025
|
246 D = 1; |
|
|
247 M(1) = 1; |
|
559
|
248 |
|
1025
|
249 old_p_index = 1; |
|
|
250 new_p_index = 2; |
|
|
251 M_index = 1; |
|
|
252 pr_index = 2; |
|
|
253 |
|
|
254 while (new_p_index <= lp) |
|
|
255 if (abs (p (new_p_index) - p (old_p_index)) < toler) |
|
2303
|
256 ## We've found a multiple pole. |
|
1025
|
257 M (M_index) = M (M_index) + 1; |
|
|
258 e (new_p_index) = e (new_p_index-1) + 1; |
|
2303
|
259 ## Remove the pole from pr. |
|
1025
|
260 pr (pr_index) = []; |
|
559
|
261 else |
|
2303
|
262 ## It's a different pole. |
|
1025
|
263 D++; |
|
|
264 M_index++; |
|
|
265 M (M_index) = 1; |
|
|
266 old_p_index = new_p_index; |
|
|
267 pr_index++; |
|
559
|
268 endif |
|
|
269 new_p_index++; |
|
|
270 endwhile |
|
|
271 |
|
2303
|
272 ## Shorten M to it's proper length |
|
559
|
273 |
|
1025
|
274 M = M (1:D); |
|
|
275 |
|
2303
|
276 ## Now set up the polynomial matrices. |
|
559
|
277 |
|
|
278 MM = max(M); |
|
1025
|
279 |
|
2303
|
280 ## Left hand side polynomial |
|
1025
|
281 |
|
|
282 lhs = zeros (MM, lb); |
|
|
283 rhs = zeros (MM, lp*lp); |
|
|
284 lhs (1, :) = b; |
|
559
|
285 rhi = 1; |
|
|
286 dpi = 1; |
|
|
287 mpi = 1; |
|
1025
|
288 while (dpi <= D) |
|
559
|
289 for ind = 1:M(dpi) |
|
1025
|
290 if (mpi > 1 && (mpi+ind) <= lp) |
|
559
|
291 cp = [p(1:mpi-1); p(mpi+ind:lp)]; |
|
1025
|
292 elseif (mpi == 1) |
|
|
293 cp = p (mpi+ind:lp); |
|
559
|
294 else |
|
1025
|
295 cp = p (1:mpi-1); |
|
559
|
296 endif |
|
1025
|
297 rhs (1, rhi:rhi+lp-1) = prepad (poly (cp), lp); |
|
559
|
298 rhi = rhi + lp; |
|
|
299 endfor |
|
1025
|
300 mpi = mpi + M (dpi); |
|
559
|
301 dpi++; |
|
|
302 endwhile |
|
1025
|
303 if (MM > 1) |
|
559
|
304 for index = 2:MM |
|
1025
|
305 lhs (index, :) = prepad (polyderiv (lhs (index-1, :)), lb); |
|
559
|
306 ind = 1; |
|
|
307 for rhi = 1:lp |
|
1025
|
308 cp = rhs (index-1, ind:ind+lp-1); |
|
|
309 rhs (index, ind:ind+lp-1) = prepad (polyderiv (cp), lp); |
|
559
|
310 ind = ind + lp; |
|
|
311 endfor |
|
|
312 endfor |
|
|
313 endif |
|
|
314 |
|
2303
|
315 ## Now lhs contains the numerator polynomial and as many derivatives as |
|
|
316 ## are required. rhs is a matrix of polynomials, the first row |
|
|
317 ## contains the corresponding polynomial for each residue and |
|
|
318 ## successive rows are derivatives. |
|
559
|
319 |
|
2303
|
320 ## Now we need to evaluate the first row of lhs and rhs at each |
|
|
321 ## distinct pole value. If there are multiple poles we will also need |
|
|
322 ## to evaluate the derivatives at the pole value also. |
|
559
|
323 |
|
1025
|
324 B = zeros (lp, 1); |
|
|
325 A = zeros (lp, lp); |
|
559
|
326 |
|
|
327 dpi = 1; |
|
|
328 row = 1; |
|
1025
|
329 while (dpi <= D) |
|
559
|
330 for mi = 1:M(dpi) |
|
1025
|
331 B (row) = polyval (lhs (mi, :), pr (dpi)); |
|
559
|
332 ci = 1; |
|
|
333 for col = 1:lp |
|
1025
|
334 cp = rhs (mi, ci:ci+lp-1); |
|
|
335 A (row, col) = polyval (cp, pr(dpi)); |
|
559
|
336 ci = ci + lp; |
|
|
337 endfor |
|
|
338 row++; |
|
|
339 endfor |
|
|
340 dpi++; |
|
|
341 endwhile |
|
|
342 |
|
2303
|
343 ## Solve for the residues. |
|
1025
|
344 |
|
|
345 r = A \ B; |
|
559
|
346 |
|
|
347 endfunction |
