Mercurial > octave-antonio
annotate scripts/polynomial/residue.m @ 8920:eb63fbe60fab
update copyright notices
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Sat, 07 Mar 2009 10:41:27 -0500 |
| parents | 81d6ab3ac93c |
| children | 1bf0ce0930be |
| rev | line source |
|---|---|
| 7017 | 1 ## Copyright (C) 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2004, 2005 |
| 8920 | 2 ## 2006, 2007, 2008, 2009 John W. Eaton |
| 6964 | 3 ## Copyright (C) 2007 Ben Abbott |
| 2313 | 4 ## |
| 5 ## This file is part of Octave. | |
| 6 ## | |
| 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 | |
| 7016 | 9 ## the Free Software Foundation; either version 3 of the License, or (at |
| 10 ## your option) any later version. | |
| 2313 | 11 ## |
| 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. | |
| 16 ## | |
| 17 ## You should have received a copy of the GNU General Public License | |
| 7016 | 18 ## along with Octave; see the file COPYING. If not, see |
| 19 ## <http://www.gnu.org/licenses/>. | |
| 904 | 20 |
| 3368 | 21 ## -*- texinfo -*- |
| 6978 | 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}. | |
| 3426 | 25 ## |
| 3368 | 26 ## @iftex |
| 27 ## @tex | |
| 28 ## $$ | |
| 6978 | 29 ## {B(s)\over A(s)} = \sum_{m=1}^M {r_m\over (s-p_m)^e_m} |
| 3368 | 30 ## + \sum_{i=1}^N k_i s^{N-i}. |
| 31 ## $$ | |
| 32 ## @end tex | |
| 33 ## @end iftex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
34 ## @ifnottex |
| 3426 | 35 ## |
| 3368 | 36 ## @example |
| 6978 | 37 ## B(s) M r(m) N |
| 3368 | 38 ## ---- = SUM ------------- + SUM k(i)*s^(N-i) |
| 6978 | 39 ## A(s) m=1 (s-p(m))^e(m) i=1 |
| 3368 | 40 ## @end example |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
41 ## @end ifnottex |
| 3426 | 42 ## |
| 3368 | 43 ## @noindent |
| 3499 | 44 ## where @math{M} is the number of poles (the length of the @var{r}, |
| 6978 | 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 | |
| 8494 | 47 ## the multiplicity of the m-th residue's pole. |
| 3426 | 48 ## |
| 6964 | 49 ## For example, |
| 3426 | 50 ## |
| 3368 | 51 ## @example |
| 52 ## @group | |
| 6964 | 53 ## b = [1, 1, 1]; |
| 54 ## a = [1, -5, 8, -4]; | |
| 7011 | 55 ## [r, p, k, e] = residue (b, a); |
| 56 ## @result{} r = [-2; 7; 3] | |
| 57 ## @result{} p = [2; 2; 1] | |
| 3368 | 58 ## @result{} k = [](0x0) |
| 7011 | 59 ## @result{} e = [1; 2; 1] |
| 3368 | 60 ## @end group |
| 61 ## @end example | |
| 3426 | 62 ## |
| 3368 | 63 ## @noindent |
| 6978 | 64 ## which represents the following partial fraction expansion |
| 3368 | 65 ## @iftex |
| 66 ## @tex | |
| 67 ## $$ | |
| 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} | |
| 69 ## $$ | |
| 70 ## @end tex | |
| 71 ## @end iftex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
72 ## @ifnottex |
| 3426 | 73 ## |
| 3368 | 74 ## @example |
| 75 ## s^2 + s + 1 -2 7 3 | |
| 2311 | 76 ## ------------------- = ----- + ------- + ----- |
| 77 ## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1) | |
| 3368 | 78 ## @end example |
| 6964 | 79 ## |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
80 ## @end ifnottex |
| 6978 | 81 ## |
| 82 ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}) | |
| 7011 | 83 ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}, @var{e}) |
| 6978 | 84 ## Compute the reconstituted quotient of polynomials, |
| 7398 | 85 ## @var{b}(s)/@var{a}(s), from the partial fraction expansion; |
| 6978 | 86 ## represented by the residues, poles, and a direct polynomial specified |
| 7011 | 87 ## by @var{r}, @var{p} and @var{k}, and the pole multiplicity @var{e}. |
| 88 ## | |
| 89 ## If the multiplicity, @var{e}, is not explicitly specified the multiplicity is | |
| 90 ## determined by the script mpoles.m. | |
| 6978 | 91 ## |
| 92 ## For example, | |
| 6964 | 93 ## |
| 94 ## @example | |
| 95 ## @group | |
| 7011 | 96 ## r = [-2; 7; 3]; |
| 97 ## p = [2; 2; 1]; | |
| 98 ## k = [1, 0]; | |
| 6964 | 99 ## [b, a] = residue (r, p, k); |
| 100 ## @result{} b = [1, -5, 9, -3, 1] | |
| 6978 | 101 ## @result{} a = [1, -5, 8, -4] |
| 7011 | 102 ## |
| 103 ## where mpoles.m is used to determine e = [1; 2; 1] | |
| 104 ## | |
| 105 ## @end group | |
| 106 ## @end example | |
| 107 ## | |
| 108 ## Alternatively the multiplicity may be defined explicitly, for example, | |
| 109 ## | |
| 110 ## @example | |
| 111 ## @group | |
| 112 ## r = [7; 3; -2]; | |
| 113 ## p = [2; 1; 2]; | |
| 114 ## k = [1, 0]; | |
| 115 ## e = [2; 1; 1]; | |
| 116 ## [b, a] = residue (r, p, k, e); | |
| 117 ## @result{} b = [1, -5, 9, -3, 1] | |
| 118 ## @result{} a = [1, -5, 8, -4] | |
| 6964 | 119 ## @end group |
| 120 ## @end example | |
| 121 ## | |
| 122 ## @noindent | |
| 6978 | 123 ## which represents the following partial fraction expansion |
| 6964 | 124 ## @iftex |
| 125 ## @tex | |
| 126 ## $$ | |
| 6978 | 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} |
| 6964 | 128 ## $$ |
| 129 ## @end tex | |
| 130 ## @end iftex | |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
131 ## @ifnottex |
| 6964 | 132 ## |
| 133 ## @example | |
| 6978 | 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 | |
| 6964 | 137 ## @end example |
|
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8494
diff
changeset
|
138 ## @end ifnottex |
| 6964 | 139 ## @seealso{poly, roots, conv, deconv, mpoles, polyval, polyderiv, polyinteg} |
| 3368 | 140 ## @end deftypefn |
| 1025 | 141 |
| 3202 | 142 ## Author: Tony Richardson <arichard@stark.cc.oh.us> |
| 6964 | 143 ## Author: Ben Abbott <bpabbott@mac.com> |
| 2312 | 144 ## Created: June 1994 |
| 145 ## Adapted-By: jwe | |
| 559 | 146 |
| 6978 | 147 function [r, p, k, e] = residue (b, a, varargin) |
| 559 | 148 |
| 7011 | 149 if (nargin < 2 || nargin > 4) |
| 6046 | 150 print_usage (); |
| 559 | 151 endif |
| 152 | |
| 6964 | 153 toler = .001; |
| 154 | |
| 7011 | 155 if (nargin >= 3) |
| 156 if (nargin >= 4) | |
| 157 e = varargin{2}; | |
| 158 else | |
| 159 e = []; | |
| 160 endif | |
| 6964 | 161 ## The inputs are the residue, pole, and direct part. Solve for the |
| 162 ## corresponding numerator and denominator polynomials | |
| 7011 | 163 [r, p] = rresidue (b, a, varargin{1}, toler, e); |
| 6964 | 164 return |
| 7011 | 165 endif |
| 559 | 166 |
| 2303 | 167 ## Make sure both polynomials are in reduced form. |
| 1025 | 168 |
| 169 a = polyreduce (a); | |
| 170 b = polyreduce (b); | |
| 559 | 171 |
| 1025 | 172 b = b / a(1); |
| 173 a = a / a(1); | |
| 559 | 174 |
| 1025 | 175 la = length (a); |
| 176 lb = length (b); | |
| 559 | 177 |
| 2303 | 178 ## Handle special cases here. |
| 1025 | 179 |
| 180 if (la == 0 || lb == 0) | |
| 559 | 181 k = r = p = e = []; |
| 182 return; | |
| 183 elseif (la == 1) | |
| 1025 | 184 k = b / a; |
| 559 | 185 r = p = e = []; |
| 186 return; | |
| 187 endif | |
| 188 | |
| 2303 | 189 ## Find the poles. |
| 1025 | 190 |
| 191 p = roots (a); | |
| 192 lp = length (p); | |
| 559 | 193 |
| 6964 | 194 ## Sort poles so that multiplicity loop will work. |
| 195 | |
| 196 [e, indx] = mpoles (p, toler, 1); | |
| 197 p = p (indx); | |
| 559 | 198 |
| 7398 | 199 ## For each group of pole multiplicity, set the value of each |
| 200 ## pole to the average of the group. This reduces the error in | |
| 201 ## the resulting poles. | |
| 202 | |
| 203 p_group = cumsum (e == 1); | |
| 204 for ng = 1:p_group(end) | |
| 205 m = find (p_group == ng); | |
| 206 p(m) = mean (p(m)); | |
| 207 endfor | |
| 208 | |
| 2303 | 209 ## Find the direct term if there is one. |
| 1025 | 210 |
| 211 if (lb >= la) | |
| 6964 | 212 ## Also return the reduced numerator. |
| 1025 | 213 [k, b] = deconv (b, a); |
| 214 lb = length (b); | |
| 559 | 215 else |
| 216 k = []; | |
| 217 endif | |
| 218 | |
| 7398 | 219 ## Determine if the poles are (effectively) zero. |
| 220 | |
| 221 small = max (abs (p)); | |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
222 if (isa (a, "single") || isa (b, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
223 small = max ([small, 1]) * eps ("single") * 1e4 * (1 + numel (p))^2; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
224 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
225 small = max ([small, 1]) * eps * 1e4 * (1 + numel (p))^2; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
226 endif |
| 7398 | 227 p(abs (p) < small) = 0; |
| 228 | |
| 229 ## Determine if the poles are (effectively) real, or imaginary. | |
| 230 | |
| 231 index = (abs (imag (p)) < small); | |
| 232 p(index) = real (p(index)); | |
| 233 index = (abs (real (p)) < small); | |
| 234 p(index) = 1i * imag (p(index)); | |
| 235 | |
| 236 ## The remainder determines the residues. The case of one pole | |
| 237 ## is trivial. | |
| 238 | |
| 1025 | 239 if (lp == 1) |
| 240 r = polyval (b, p); | |
| 559 | 241 return; |
| 242 endif | |
| 243 | |
| 6964 | 244 ## Determine the order of the denominator and remaining numerator. |
| 245 ## With the direct term removed the potential order of the numerator | |
| 246 ## is one less than the order of the denominator. | |
| 1025 | 247 |
| 6964 | 248 aorder = numel (a) - 1; |
| 249 border = aorder - 1; | |
| 1025 | 250 |
| 6964 | 251 ## Construct a system of equations relating the individual |
| 252 ## contributions from each residue to the complete numerator. | |
| 559 | 253 |
| 6964 | 254 A = zeros (border+1, border+1); |
| 255 B = prepad (reshape (b, [numel(b), 1]), border+1, 0); | |
| 256 for ip = 1:numel(p) | |
| 257 ri = zeros (size (p)); | |
| 258 ri(ip) = 1; | |
| 259 A(:,ip) = prepad (rresidue (ri, p, [], toler), border+1, 0).'; | |
| 260 endfor | |
| 559 | 261 |
| 2303 | 262 ## Solve for the residues. |
| 1025 | 263 |
| 264 r = A \ B; | |
| 559 | 265 |
| 266 endfunction | |
| 6964 | 267 |
| 7011 | 268 function [pnum, pden, e] = rresidue (r, p, k, toler, e) |
| 6964 | 269 |
| 270 ## Reconstitute the numerator and denominator polynomials from the | |
| 271 ## residues, poles, and direct term. | |
| 272 | |
| 7011 | 273 if (nargin < 2 || nargin > 5) |
| 6964 | 274 print_usage (); |
| 275 endif | |
| 276 | |
| 7011 | 277 if (nargin < 5) |
| 278 e = []; | |
| 279 endif | |
| 280 | |
| 6964 | 281 if (nargin < 4) |
| 282 toler = []; | |
| 283 endif | |
| 284 | |
| 285 if (nargin < 3) | |
| 286 k = []; | |
| 287 endif | |
| 7011 | 288 |
| 289 if numel (e) | |
| 290 indx = 1:numel(p); | |
| 291 else | |
| 292 [e, indx] = mpoles (p, toler, 0); | |
| 293 p = p (indx); | |
| 294 r = r (indx); | |
| 295 endif | |
| 6964 | 296 |
| 297 indx = 1:numel(p); | |
| 298 | |
| 299 for n = indx | |
| 300 pn = [1, -p(n)]; | |
| 301 if n == 1 | |
| 302 pden = pn; | |
| 303 else | |
| 304 pden = conv (pden, pn); | |
| 305 endif | |
| 306 endfor | |
| 307 | |
| 308 ## D is the order of the denominator | |
| 309 ## K is the order of the direct polynomial | |
| 310 ## N is the order of the resulting numerator | |
| 311 ## pnum(1:(N+1)) is the numerator's polynomial | |
| 312 ## pden(1:(D+1)) is the denominator's polynomial | |
| 313 ## pm is the multible pole for the nth residue | |
| 314 ## pn is the numerator contribution for the nth residue | |
| 315 | |
| 316 D = numel (pden) - 1; | |
| 317 K = numel (k) - 1; | |
| 318 N = K + D; | |
| 319 pnum = zeros (1, N+1); | |
| 7011 | 320 for n = indx(abs (r) > 0) |
| 6964 | 321 p1 = [1, -p(n)]; |
| 7011 | 322 for m = 1:e(n) |
| 323 if (m == 1) | |
| 6964 | 324 pm = p1; |
| 325 else | |
| 326 pm = conv (pm, p1); | |
| 327 endif | |
| 328 endfor | |
| 329 pn = deconv (pden, pm); | |
| 330 pn = r(n) * pn; | |
| 7183 | 331 pnum = pnum + prepad (pn, N+1, 0, 2); |
| 6964 | 332 endfor |
| 333 | |
| 334 ## Add the direct term. | |
| 335 | |
| 336 if (numel (k)) | |
| 337 pnum = pnum + conv (pden, k); | |
| 338 endif | |
| 339 | |
| 340 ## Check for leading zeros and trim the polynomial coefficients. | |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
341 if (isa (r, "single") || isa (p, "single") || isa (k, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
342 small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps ("single"); |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
343 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
344 small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7398
diff
changeset
|
345 endif |
| 6964 | 346 |
| 7011 | 347 pnum(abs (pnum) < small) = 0; |
| 348 pden(abs (pden) < small) = 0; | |
| 6964 | 349 |
| 350 pnum = polyreduce (pnum); | |
| 351 pden = polyreduce (pden); | |
| 352 | |
| 353 endfunction | |
| 6968 | 354 |
| 355 %!test | |
| 356 %! b = [1, 1, 1]; | |
| 357 %! a = [1, -5, 8, -4]; | |
| 358 %! [r, p, k, e] = residue (b, a); | |
| 7398 | 359 %! assert (abs (r - [-2; 7; 3]) < 1e-12 |
| 360 %! && abs (p - [2; 2; 1]) < 1e-12 | |
| 6998 | 361 %! && isempty (k) |
| 362 %! && e == [1; 2; 1]); | |
| 6994 | 363 %! k = [1 0]; |
| 7011 | 364 %! b = conv (k, a) + prepad (b, numel (k) + numel (a) - 1, 0); |
| 365 %! a = a; | |
| 366 %! [br, ar] = residue (r, p, k); | |
| 367 %! assert ((abs (br - b) < 1e-12 | |
| 368 %! && abs (ar - a) < 1e-12)); | |
| 369 %! [br, ar] = residue (r, p, k, e); | |
| 370 %! assert ((abs (br - b) < 1e-12 | |
| 371 %! && abs (ar - a) < 1e-12)); | |
| 6994 | 372 |
| 373 %!test | |
| 374 %! b = [1, 0, 1]; | |
| 375 %! a = [1, 0, 18, 0, 81]; | |
| 7398 | 376 %! [r, p, k, e] = residue (b, a); |
| 6994 | 377 %! r1 = [-5i; 12; +5i; 12]/54; |
| 378 %! p1 = [+3i; +3i; -3i; -3i]; | |
| 7398 | 379 %! assert (abs (r - r1) < 1e-12 && abs (p - p1) < 1e-12 |
| 6998 | 380 %! && isempty (k) |
| 381 %! && e == [1; 2; 1; 2]); | |
| 6994 | 382 %! [br, ar] = residue (r, p, k); |
| 383 %! assert ((abs (br - b) < 1e-12 | |
| 6998 | 384 %! && abs (ar - a) < 1e-12)); |
| 7011 | 385 |
| 386 %!test | |
| 387 %! r = [7; 3; -2]; | |
| 388 %! p = [2; 1; 2]; | |
| 389 %! k = [1 0]; | |
| 390 %! e = [2; 1; 1]; | |
| 391 %! [b, a] = residue (r, p, k, e); | |
| 392 %! assert ((abs (b - [1, -5, 9, -3, 1]) < 1e-12 | |
| 393 %! && abs (a - [1, -5, 8, -4]) < 1e-12)); | |
| 394 %! [rr, pr, kr, er] = residue (b, a); | |
| 395 %! [jnk, n] = mpoles(p); | |
| 7398 | 396 %! assert ((abs (rr - r(n)) < 1e-12 |
| 397 %! && abs (pr - p(n)) < 1e-12 | |
| 7011 | 398 %! && abs (kr - k) < 1e-12 |
| 399 %! && abs (er - e(n)) < 1e-12)); | |
| 400 | |
| 7188 | 401 %!test |
| 402 %! b = [1]; | |
| 403 %! a = [1, 10, 25]; | |
| 7398 | 404 %! [r, p, k, e] = residue (b, a); |
| 7188 | 405 %! r1 = [0; 1]; |
| 406 %! p1 = [-5; -5]; | |
| 7398 | 407 %! assert (abs (r - r1) < 1e-12 && abs (p - p1) < 1e-12 |
| 7188 | 408 %! && isempty (k) |
| 409 %! && e == [1; 2]); | |
| 410 %! [br, ar] = residue (r, p, k); | |
| 411 %! assert ((abs (br - b) < 1e-12 | |
| 412 %! && abs (ar - a) < 1e-12)); |