octave_packages/m/polynomial/residue.m

   1 ## Copyright (C) 1994-2012 John W. Eaton
   2 ## Copyright (C) 2007 Ben Abbott
   3 ##
   4 ## This file is part of Octave.
   5 ##
   6 ## Octave is free software; you can redistribute it and/or modify it
   7 ## under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation; either version 3 of the License, or (at
   9 ## your option) any later version.
  10 ##
  11 ## Octave is distributed in the hope that it will be useful, but
  12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 ## General Public License for more details.
  15 ##
  16 ## You should have received a copy of the GNU General Public License
  17 ## along with Octave; see the file COPYING.  If not, see
  18 ## <http://www.gnu.org/licenses/>.
  19
  20 ## -*- texinfo -*-
  21 ## @deftypefn  {Function File} {[@var{r}, @var{p}, @var{k}, @var{e}] =} residue (@var{b}, @var{a})
  22 ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k})
  23 ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}, @var{e})
  24 ## The first calling form computes the partial fraction expansion for the
  25 ## quotient of the polynomials, @var{b} and @var{a}.
  26 ## @tex
  27 ## $$
  28 ## {B(s)\over A(s)} = \sum_{m=1}^M {r_m\over (s-p_m)^e_m}
  29 ##   + \sum_{i=1}^N k_i s^{N-i}.
  30 ## $$
  31 ## @end tex
  32 ## @ifnottex
  33 ##
  34 ## @example
  35 ## @group
  36 ## B(s)    M       r(m)         N
  37 ## ---- = SUM -------------  + SUM k(i)*s^(N-i)
  38 ## A(s)   m=1 (s-p(m))^e(m)    i=1
  39 ## @end group
  40 ## @end example
  41 ##
  42 ## @end ifnottex
  43 ## @noindent
  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 m-th residue's pole.
  48 ##
  49 ## For example,
  50 ##
  51 ## @example
  52 ## @group
  53 ## b = [1, 1, 1];
  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]
  60 ## @end group
  61 ## @end example
  62 ##
  63 ## @noindent
  64 ## which represents the following partial fraction expansion
  65 ## @tex
  66 ## $$
  67 ## {s^2+s+1\over s^3-5s^2+8s-4} = {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1}
  68 ## $$
  69 ## @end tex
  70 ## @ifnottex
  71 ##
  72 ## @example
  73 ## @group
  74 ##         s^2 + s + 1       -2        7        3
  75 ##    ------------------- = ----- + ------- + -----
  76 ##    s^3 - 5s^2 + 8s - 4   (s-2)   (s-2)^2   (s-1)
  77 ## @end group
  78 ## @end example
  79 ##
  80 ## @end ifnottex
  81 ##
  82 ## The second calling form performs the inverse operation and computes
  83 ## the reconstituted quotient of polynomials, @var{b}(s)/@var{a}(s),
  84 ## from the partial fraction expansion; represented by the residues,
  85 ## poles, and a direct polynomial specified by @var{r}, @var{p} and
  86 ## @var{k}, and the pole multiplicity @var{e}.
  87 ##
  88 ## If the multiplicity, @var{e}, is not explicitly specified the multiplicity is
  89 ## determined by the function @code{mpoles}.
  90 ##
  91 ## For example:
  92 ##
  93 ## @example
  94 ## @group
  95 ## r = [-2; 7; 3];
  96 ## p = [2; 2; 1];
  97 ## k = [1, 0];
  98 ## [b, a] = residue (r, p, k)
  99 ##    @result{} b = [1, -5, 9, -3, 1]
 100 ##    @result{} a = [1, -5, 8, -4]
 101 ##
 102 ## where mpoles is used to determine e = [1; 2; 1]
 103 ## @end group
 104 ## @end example
 105 ##
 106 ## Alternatively the multiplicity may be defined explicitly, for example,
 107 ##
 108 ## @example
 109 ## @group
 110 ## r = [7; 3; -2];
 111 ## p = [2; 1; 2];
 112 ## k = [1, 0];
 113 ## e = [2; 1; 1];
 114 ## [b, a] = residue (r, p, k, e)
 115 ##    @result{} b = [1, -5, 9, -3, 1]
 116 ##    @result{} a = [1, -5, 8, -4]
 117 ## @end group
 118 ## @end example
 119 ##
 120 ## @noindent
 121 ## which represents the following partial fraction expansion
 122 ## @tex
 123 ## $$
 124 ## {-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}
 125 ## $$
 126 ## @end tex
 127 ## @ifnottex
 128 ##
 129 ## @example
 130 ## @group
 131 ##  -2        7        3         s^4 - 5s^3 + 9s^2 - 3s + 1
 132 ## ----- + ------- + ----- + s = --------------------------
 133 ## (s-2)   (s-2)^2   (s-1)          s^3 - 5s^2 + 8s - 4
 134 ## @end group
 135 ## @end example
 136 ##
 137 ## @end ifnottex
 138 ## @seealso{mpoles, poly, roots, conv, deconv}
 139 ## @end deftypefn
 140
 141 ## Author: Tony Richardson <arichard@stark.cc.oh.us>
 142 ## Author: Ben Abbott <bpabbott@mac.com>
 143 ## Created: June 1994
 144 ## Adapted-By: jwe
 145
 146 function [r, p, k, e] = residue (b, a, varargin)
 147
 148   if (nargin < 2 || nargin > 4)
 149     print_usage ();
 150   endif
 151
 152   toler = .001;
 153
 154   if (nargin >= 3)
 155     if (nargin >= 4)
 156       e = varargin{2};
 157     else
 158       e = [];
 159     endif
 160     ## The inputs are the residue, pole, and direct part. Solve for the
 161     ## corresponding numerator and denominator polynomials
 162     [r, p] = rresidue (b, a, varargin{1}, toler, e);
 163     return
 164   endif
 165
 166   ## Make sure both polynomials are in reduced form.
 167
 168   a = polyreduce (a);
 169   b = polyreduce (b);
 170
 171   b = b / a(1);
 172   a = a / a(1);
 173
 174   la = length (a);
 175   lb = length (b);
 176
 177   ## Handle special cases here.
 178
 179   if (la == 0 || lb == 0)
 180     k = r = p = e = [];
 181     return;
 182   elseif (la == 1)
 183     k = b / a;
 184     r = p = e = [];
 185     return;
 186   endif
 187
 188   ## Find the poles.
 189
 190   p = roots (a);
 191   lp = length (p);
 192
 193   ## Sort poles so that multiplicity loop will work.
 194
 195   [e, indx] = mpoles (p, toler, 1);
 196   p = p (indx);
 197
 198   ## For each group of pole multiplicity, set the value of each
 199   ## pole to the average of the group. This reduces the error in
 200   ## the resulting poles.
 201
 202   p_group = cumsum (e == 1);
 203   for ng = 1:p_group(end)
 204     m = find (p_group == ng);
 205     p(m) = mean (p(m));
 206   endfor
 207
 208   ## Find the direct term if there is one.
 209
 210   if (lb >= la)
 211     ## Also return the reduced numerator.
 212     [k, b] = deconv (b, a);
 213     lb = length (b);
 214   else
 215     k = [];
 216   endif
 217
 218   ## Determine if the poles are (effectively) zero.
 219
 220   small = max (abs (p));
 221   if (isa (a, "single") || isa (b, "single"))
 222     small = max ([small, 1]) * eps ("single") * 1e4 * (1 + numel (p))^2;
 223   else
 224     small = max ([small, 1]) * eps * 1e4 * (1 + numel (p))^2;
 225   endif
 226   p(abs (p) < small) = 0;
 227
 228   ## Determine if the poles are (effectively) real, or imaginary.
 229
 230   index = (abs (imag (p)) < small);
 231   p(index) = real (p(index));
 232   index = (abs (real (p)) < small);
 233   p(index) = 1i * imag (p(index));
 234
 235   ## The remainder determines the residues.  The case of one pole
 236   ## is trivial.
 237
 238   if (lp == 1)
 239     r = polyval (b, p);
 240     return;
 241   endif
 242
 243   ## Determine the order of the denominator and remaining numerator.
 244   ## With the direct term removed the potential order of the numerator
 245   ## is one less than the order of the denominator.
 246
 247   aorder = numel (a) - 1;
 248   border = aorder - 1;
 249
 250   ## Construct a system of equations relating the individual
 251   ## contributions from each residue to the complete numerator.
 252
 253   A = zeros (border+1, border+1);
 254   B = prepad (reshape (b, [numel(b), 1]), border+1, 0);
 255   for ip = 1:numel(p)
 256     ri = zeros (size (p));
 257     ri(ip) = 1;
 258     A(:,ip) = prepad (rresidue (ri, p, [], toler), border+1, 0).';
 259   endfor
 260
 261   ## Solve for the residues.
 262
 263   r = A \ B;
 264
 265 endfunction
 266
 267 function [pnum, pden, e] = rresidue (r, p, k, toler, e)
 268
 269   ## Reconstitute the numerator and denominator polynomials from the
 270   ## residues, poles, and direct term.
 271
 272   if (nargin < 2 || nargin > 5)
 273     print_usage ();
 274   endif
 275
 276   if (nargin < 5)
 277     e = [];
 278   endif
 279
 280   if (nargin < 4)
 281     toler = [];
 282   endif
 283
 284   if (nargin < 3)
 285     k = [];
 286   endif
 287
 288   if numel (e)
 289     indx = 1:numel(p);
 290   else
 291     [e, indx] = mpoles (p, toler, 0);
 292     p = p (indx);
 293     r = r (indx);
 294   endif
 295
 296   indx = 1:numel(p);
 297
 298   for n = indx
 299     pn = [1, -p(n)];
 300     if n == 1
 301       pden = pn;
 302     else
 303       pden = conv (pden, pn);
 304     endif
 305   endfor
 306
 307   ## D is the order of the denominator
 308   ## K is the order of the direct polynomial
 309   ## N is the order of the resulting numerator
 310   ## pnum(1:(N+1)) is the numerator's polynomial
 311   ## pden(1:(D+1)) is the denominator's polynomial
 312   ## pm is the multible pole for the nth residue
 313   ## pn is the numerator contribution for the nth residue
 314
 315   D = numel (pden) - 1;
 316   K = numel (k) - 1;
 317   N = K + D;
 318   pnum = zeros (1, N+1);
 319   for n = indx(abs (r) > 0)
 320     p1 = [1, -p(n)];
 321     for m = 1:e(n)
 322       if (m == 1)
 323         pm = p1;
 324       else
 325         pm = conv (pm, p1);
 326       endif
 327     endfor
 328     pn = deconv (pden, pm);
 329     pn = r(n) * pn;
 330     pnum = pnum + prepad (pn, N+1, 0, 2);
 331   endfor
 332
 333   ## Add the direct term.
 334
 335   if (numel (k))
 336     pnum = pnum + conv (pden, k);
 337   endif
 338
 339   ## Check for leading zeros and trim the polynomial coefficients.
 340   if (isa (r, "single") || isa (p, "single") || isa (k, "single"))
 341     small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps ("single");
 342   else
 343     small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps;
 344   endif
 345
 346   pnum(abs (pnum) < small) = 0;
 347   pden(abs (pden) < small) = 0;
 348
 349   pnum = polyreduce (pnum);
 350   pden = polyreduce (pden);
 351
 352 endfunction
 353
 354 %!test
 355 %! b = [1, 1, 1];
 356 %! a = [1, -5, 8, -4];
 357 %! [r, p, k, e] = residue (b, a);
 358 %! assert (abs (r - [-2; 7; 3]) < 1e-12
 359 %!   && abs (p - [2; 2; 1]) < 1e-12
 360 %!   && isempty (k)
 361 %!   && e == [1; 2; 1]);
 362 %! k = [1 0];
 363 %! b = conv (k, a) + prepad (b, numel (k) + numel (a) - 1, 0);
 364 %! a = a;
 365 %! [br, ar] = residue (r, p, k);
 366 %! assert ((abs (br - b) < 1e-12
 367 %!   && abs (ar - a) < 1e-12));
 368 %! [br, ar] = residue (r, p, k, e);
 369 %! assert ((abs (br - b) < 1e-12
 370 %!   && abs (ar - a) < 1e-12));
 371
 372 %!test
 373 %! b = [1, 0, 1];
 374 %! a = [1, 0, 18, 0, 81];
 375 %! [r, p, k, e] = residue (b, a);
 376 %! r1 = [-5i; 12; +5i; 12]/54;
 377 %! p1 = [+3i; +3i; -3i; -3i];
 378 %! assert (abs (r - r1) < 1e-12 && abs (p - p1) < 1e-12
 379 %!   && isempty (k)
 380 %!   && e == [1; 2; 1; 2]);
 381 %! [br, ar] = residue (r, p, k);
 382 %! assert ((abs (br - b) < 1e-12
 383 %!   && abs (ar - a) < 1e-12));
 384
 385 %!test
 386 %! r = [7; 3; -2];
 387 %! p = [2; 1; 2];
 388 %! k = [1 0];
 389 %! e = [2; 1; 1];
 390 %! [b, a] = residue (r, p, k, e);
 391 %! assert ((abs (b - [1, -5, 9, -3, 1]) < 1e-12
 392 %!   && abs (a - [1, -5, 8, -4]) < 1e-12));
 393 %! [rr, pr, kr, er] = residue (b, a);
 394 %! [jnk, n] = mpoles(p);
 395 %! assert ((abs (rr - r(n)) < 1e-12
 396 %!   && abs (pr - p(n)) < 1e-12
 397 %!   && abs (kr - k) < 1e-12
 398 %!   && abs (er - e(n)) < 1e-12));
 399
 400 %!test
 401 %! b = [1];
 402 %! a = [1, 10, 25];
 403 %! [r, p, k, e] = residue (b, a);
 404 %! r1 = [0; 1];
 405 %! p1 = [-5; -5];
 406 %! assert (abs (r - r1) < 1e-12 && abs (p - p1) < 1e-12
 407 %!   && isempty (k)
 408 %!   && e == [1; 2]);
 409 %! [br, ar] = residue (r, p, k);
 410 %! assert ((abs (br - b) < 1e-12
 411 %!   && abs (ar - a) < 1e-12));
 412
 413 ## The following test is due to Bernard Grung (bug #34266)
 414 %!xtest
 415 %! z1 =  7.0372976777e6;
 416 %! p1 = -3.1415926536e9;
 417 %! p2 = -4.9964813512e8;
 418 %! r1 = -(1 + z1/p1)/(1 - p1/p2)/p2/p1;
 419 %! r2 = -(1 + z1/p2)/(1 - p2/p1)/p2/p1;
 420 %! r3 = (1 + (p2 + p1)/p2/p1*z1)/p2/p1;
 421 %! r4 = z1/p2/p1;
 422 %! r = [r1; r2; r3; r4];
 423 %! p = [p1; p2; 0; 0];
 424 %! k = [];
 425 %! e = [1; 1; 1; 2];
 426 %! b = [1, z1];
 427 %! a = [1, -(p1 + p2), p1*p2, 0, 0];
 428 %! [br, ar] = residue (r, p, k, e);
 429 %! assert (br, b, 1e-8);
 430 %! assert (ar, a, 1e-8);