1
## Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2004, 2005
2
## 2006, 2007 John W. Eaton
3
## Copyright (C) 2007 Ben Abbott
5
## This file is part of Octave.
7
## Octave is free software; you can redistribute it and/or modify it
8
## under the terms of the GNU General Public License as published by
9
## the Free Software Foundation; either version 3 of the License, or (at
10
## your option) any later version.
12
## Octave is distributed in the hope that it will be useful, but
13
## WITHOUT ANY WARRANTY; without even the implied warranty of
14
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15
## General Public License for more details.
17
## You should have received a copy of the GNU General Public License
18
## along with Octave; see the file COPYING. If not, see
19
## <http://www.gnu.org/licenses/>.
22
## @deftypefn {Function File} {[@var{r}, @var{p}, @var{k}, @var{e}] =} residue (@var{b}, @var{a})
23
## Compute the partial fraction expansion for the quotient of the
24
## polynomials, @var{b} and @var{a}.
29
## {B(s)\over A(s)} = \sum_{m=1}^M {r_m\over (s-p_m)^e_m}
30
## + \sum_{i=1}^N k_i s^{N-i}.
38
## ---- = SUM ------------- + SUM k(i)*s^(N-i)
39
## A(s) m=1 (s-p(m))^e(m) i=1
44
## where @math{M} is the number of poles (the length of the @var{r},
45
## @var{p}, and @var{e}), the @var{k} vector is a polynomial of order @math{N-1}
46
## representing the direct contribution, and the @var{e} vector specifies
47
## the multiplicity of the mth residue's pole.
54
## a = [1, -5, 8, -4];
55
## [r, p, k, e] = residue (b, a);
56
## @result{} r = [-2; 7; 3]
57
## @result{} p = [2; 2; 1]
58
## @result{} k = [](0x0)
59
## @result{} e = [1; 2; 1]
64
## which represents the following partial fraction expansion
68
## {s^2+s+1\over s^3-5s^2+8s-4} = {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1}
76
## ------------------- = ----- + ------- + -----
77
## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1)
82
## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k})
83
## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}, @var{e})
84
## Compute the reconstituted quotient of polynomials,
85
## @var{b}(s)/@var{a}(s), from the partial fraction expansion
86
## represented by the residues, poles, and a direct polynomial specified
87
## by @var{r}, @var{p} and @var{k}, and the pole multiplicity @var{e}.
89
## If the multiplicity, @var{e}, is not explicitly specified the multiplicity is
90
## determined by the script mpoles.m.
99
## [b, a] = residue (r, p, k);
100
## @result{} b = [1, -5, 9, -3, 1]
101
## @result{} a = [1, -5, 8, -4]
103
## where mpoles.m is used to determine e = [1; 2; 1]
108
## Alternatively the multiplicity may be defined explicitly, for example,
116
## [b, a] = residue (r, p, k, e);
117
## @result{} b = [1, -5, 9, -3, 1]
118
## @result{} a = [1, -5, 8, -4]
123
## which represents the following partial fraction expansion
127
## {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} + s = {s^4-5s^3+9s^2-3s+1\over s^3-5s^2+8s-4}
134
## -2 7 3 s^4 - 5s^3 + 9s^2 - 3s + 1
135
## ----- + ------- + ----- + s = --------------------------
136
## (s-2) (s-2)^2 (s-1) s^3 - 5s^2 + 8s - 4
139
## @seealso{poly, roots, conv, deconv, mpoles, polyval, polyderiv, polyinteg}
142
## Author: Tony Richardson <arichard@stark.cc.oh.us>
143
## Author: Ben Abbott <bpabbott@mac.com>
144
## Created: June 1994
147
function [r, p, k, e] = residue (b, a, varargin)
149
if (nargin < 2 || nargin > 4)
161
## The inputs are the residue, pole, and direct part. Solve for the
162
## corresponding numerator and denominator polynomials
163
[r, p] = rresidue (b, a, varargin{1}, toler, e);
167
## Make sure both polynomials are in reduced form.
178
## Handle special cases here.
180
if (la == 0 || lb == 0)
194
## Determine if the poles are (effectively) zero.
196
small = max (abs (p));
197
small = max ([small, 1] ) * 1e-8 * (1 + numel (p))^2;
198
p(abs (p) < small) = 0;
200
## Determine if the poles are (effectively) real, or imaginary.
202
index = (abs (imag (p)) < small);
203
p(index) = real (p(index));
204
index = (abs (real (p)) < small);
205
p(index) = 1i * imag (p(index));
207
## Sort poles so that multiplicity loop will work.
209
[e, indx] = mpoles (p, toler, 1);
212
## Find the direct term if there is one.
215
## Also return the reduced numerator.
216
[k, b] = deconv (b, a);
227
## Determine the order of the denominator and remaining numerator.
228
## With the direct term removed the potential order of the numerator
229
## is one less than the order of the denominator.
231
aorder = numel (a) - 1;
234
## Construct a system of equations relating the individual
235
## contributions from each residue to the complete numerator.
237
A = zeros (border+1, border+1);
238
B = prepad (reshape (b, [numel(b), 1]), border+1, 0);
240
ri = zeros (size (p));
242
A(:,ip) = prepad (rresidue (ri, p, [], toler), border+1, 0).';
245
## Solve for the residues.
251
function [pnum, pden, e] = rresidue (r, p, k, toler, e)
253
## Reconstitute the numerator and denominator polynomials from the
254
## residues, poles, and direct term.
256
if (nargin < 2 || nargin > 5)
275
[e, indx] = mpoles (p, toler, 0);
287
pden = conv (pden, pn);
291
## D is the order of the denominator
292
## K is the order of the direct polynomial
293
## N is the order of the resulting numerator
294
## pnum(1:(N+1)) is the numerator's polynomial
295
## pden(1:(D+1)) is the denominator's polynomial
296
## pm is the multible pole for the nth residue
297
## pn is the numerator contribution for the nth residue
299
D = numel (pden) - 1;
302
pnum = zeros (1, N+1);
303
for n = indx(abs (r) > 0)
312
pn = deconv (pden, pm);
314
pnum = pnum + prepad (pn, N+1, 0, 2);
317
## Add the direct term.
320
pnum = pnum + conv (pden, k);
323
## Check for leading zeros and trim the polynomial coefficients.
325
small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps;
327
pnum(abs (pnum) < small) = 0;
328
pden(abs (pden) < small) = 0;
330
pnum = polyreduce (pnum);
331
pden = polyreduce (pden);
337
%! a = [1, -5, 8, -4];
338
%! [r, p, k, e] = residue (b, a);
339
%! assert (abs (r - [-2; 7; 3]) < 1e-5
340
%! && abs (p - [2; 2; 1]) < 1e-7
342
%! && e == [1; 2; 1]);
344
%! b = conv (k, a) + prepad (b, numel (k) + numel (a) - 1, 0);
346
%! [br, ar] = residue (r, p, k);
347
%! assert ((abs (br - b) < 1e-12
348
%! && abs (ar - a) < 1e-12));
349
%! [br, ar] = residue (r, p, k, e);
350
%! assert ((abs (br - b) < 1e-12
351
%! && abs (ar - a) < 1e-12));
355
%! a = [1, 0, 18, 0, 81];
356
%! [r, p, k, e] = residue(b, a);
357
%! r1 = [-5i; 12; +5i; 12]/54;
358
%! p1 = [+3i; +3i; -3i; -3i];
359
%! assert (abs (r - r1) < 1e-7 && abs (p - p1) < 1e-7
361
%! && e == [1; 2; 1; 2]);
362
%! [br, ar] = residue (r, p, k);
363
%! assert ((abs (br - b) < 1e-12
364
%! && abs (ar - a) < 1e-12));
371
%! [b, a] = residue (r, p, k, e);
372
%! assert ((abs (b - [1, -5, 9, -3, 1]) < 1e-12
373
%! && abs (a - [1, -5, 8, -4]) < 1e-12));
374
%! [rr, pr, kr, er] = residue (b, a);
375
%! [jnk, n] = mpoles(p);
376
%! assert ((abs (rr - r(n)) < 1e-5
377
%! && abs (pr - p(n)) < 1e-7
378
%! && abs (kr - k) < 1e-12
379
%! && abs (er - e(n)) < 1e-12));
384
%! [r, p, k, e] = residue(b, a);
387
%! assert (abs (r - r1) < 1e-7 && abs (p - p1) < 1e-7
390
%! [br, ar] = residue (r, p, k);
391
%! assert ((abs (br - b) < 1e-12
392
%! && abs (ar - a) < 1e-12));
added
removed
scripts/polynomial/residue.m
