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