Mercurial > hg > octave-lyh
changeset 7011:4a682c7b2bd6
[project @ 2007-10-11 17:54:48 by jwe]
author | jwe |
---|---|
date | Thu, 11 Oct 2007 17:54:48 +0000 |
parents | f7d2f54008f5 |
children | fa4b43705e37 |
files | scripts/ChangeLog scripts/polynomial/residue.m |
diffstat | 2 files changed, 88 insertions(+), 28 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,8 @@ +2007-10-11 Ben Abbott <bpabbott@mac.com> + + * polynomial/residue.m: New optional input for pole multiplicity. + Doc fix. Fix tests. + 2007-10-11 Thomas Treichl <Thomas.Treichl@gmx.net> * toplev.cc (Foctave_config_info): Add field "mac".
--- a/scripts/polynomial/residue.m +++ b/scripts/polynomial/residue.m @@ -52,10 +52,11 @@ ## @group ## b = [1, 1, 1]; ## a = [1, -5, 8, -4]; -## [r, p, k] = residue (b, a); -## @result{} r = [-2, 7, 3] -## @result{} p = [2, 2, 1] +## [r, p, k, e] = residue (b, a); +## @result{} r = [-2; 7; 3] +## @result{} p = [2; 2; 1] ## @result{} k = [](0x0) +## @result{} e = [1; 2; 1] ## @end group ## @end example ## @@ -79,21 +80,42 @@ ## @end ifinfo ## ## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}) +## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}, @var{e}) ## Compute the reconstituted quotient of polynomials, ## @var{b}(s)/@var{a}(s), from the partial fraction expansion ## represented by the residues, poles, and a direct polynomial specified -## by @var{r}, @var{p} and @var{k}. +## by @var{r}, @var{p} and @var{k}, and the pole multiplicity @var{e}. +## +## If the multiplicity, @var{e}, is not explicitly specified the multiplicity is +## determined by the script mpoles.m. ## ## For example, ## ## @example ## @group -## r = [-2, 7, 3]; -## p = [2, 2, 1]; -## k = [1 0]; +## r = [-2; 7; 3]; +## p = [2; 2; 1]; +## k = [1, 0]; ## [b, a] = residue (r, p, k); ## @result{} b = [1, -5, 9, -3, 1] ## @result{} a = [1, -5, 8, -4] +## +## where mpoles.m is used to determine e = [1; 2; 1] +## +## @end group +## @end example +## +## Alternatively the multiplicity may be defined explicitly, for example, +## +## @example +## @group +## r = [7; 3; -2]; +## p = [2; 1; 2]; +## k = [1, 0]; +## e = [2; 1; 1]; +## [b, a] = residue (r, p, k, e); +## @result{} b = [1, -5, 9, -3, 1] +## @result{} a = [1, -5, 8, -4] ## @end group ## @end example ## @@ -124,18 +146,23 @@ function [r, p, k, e] = residue (b, a, varargin) - if (nargin < 2 || nargin > 3) + if (nargin < 2 || nargin > 4) print_usage (); endif toler = .001; - if (nargin == 3) + if (nargin >= 3) + if (nargin >= 4) + e = varargin{2}; + else + e = []; + endif ## The inputs are the residue, pole, and direct part. Solve for the ## corresponding numerator and denominator polynomials - [r, p] = rresidue (b, a, varargin{1}, toler); + [r, p] = rresidue (b, a, varargin{1}, toler, e); return - end + endif ## Make sure both polynomials are in reduced form. @@ -221,15 +248,19 @@ endfunction -function [pnum, pden, multp] = rresidue (r, p, k, toler) +function [pnum, pden, e] = rresidue (r, p, k, toler, e) ## Reconstitute the numerator and denominator polynomials from the ## residues, poles, and direct term. - if (nargin < 2 || nargin > 4) + if (nargin < 2 || nargin > 5) print_usage (); endif + if (nargin < 5) + e = []; + endif + if (nargin < 4) toler = []; endif @@ -237,11 +268,14 @@ if (nargin < 3) k = []; endif - - [multp, indx] = mpoles (p, toler, 0); - - p = p (indx); - r = r (indx); + + if numel (e) + indx = 1:numel(p); + else + [e, indx] = mpoles (p, toler, 0); + p = p (indx); + r = r (indx); + endif indx = 1:numel(p); @@ -266,10 +300,10 @@ K = numel (k) - 1; N = K + D; pnum = zeros (1, N+1); - for n = indx(abs(r)>0) + for n = indx(abs (r) > 0) p1 = [1, -p(n)]; - for m = 1:multp(n) - if m == 1 + for m = 1:e(n) + if (m == 1) pm = p1; else pm = conv (pm, p1); @@ -277,7 +311,7 @@ endfor pn = deconv (pden, pm); pn = r(n) * pn; - pnum = pnum + prepad ( pn, N+1, 0); + pnum = pnum + prepad (pn, N+1, 0); endfor ## Add the direct term. @@ -290,8 +324,8 @@ small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps; - pnum (abs (pnum) < small) = 0; - pden (abs (pden) < small) = 0; + pnum(abs (pnum) < small) = 0; + pden(abs (pden) < small) = 0; pnum = polyreduce (pnum); pden = polyreduce (pden); @@ -307,16 +341,20 @@ %! && isempty (k) %! && e == [1; 2; 1]); %! k = [1 0]; -%! [b, a] = residue (r, p, k); -%! assert ((abs (b - [1, -5, 9, -3, 1]) < 1e-12 -%! && abs (a - [1, -5, 8, -4]) < 1e-12)); +%! b = conv (k, a) + prepad (b, numel (k) + numel (a) - 1, 0); +%! a = a; +%! [br, ar] = residue (r, p, k); +%! assert ((abs (br - b) < 1e-12 +%! && abs (ar - a) < 1e-12)); +%! [br, ar] = residue (r, p, k, e); +%! assert ((abs (br - b) < 1e-12 +%! && abs (ar - a) < 1e-12)); %!test %! b = [1, 0, 1]; %! a = [1, 0, 18, 0, 81]; %! [r, p, k, e] = residue(b, a); %! r1 = [-5i; 12; +5i; 12]/54; -%! r2 = conj(r1); %! p1 = [+3i; +3i; -3i; -3i]; %! assert (abs (r - r1) < 1e-7 && abs (p - p1) < 1e-7 %! && isempty (k) @@ -324,3 +362,20 @@ %! [br, ar] = residue (r, p, k); %! assert ((abs (br - b) < 1e-12 %! && abs (ar - a) < 1e-12)); + +%!test +%! r = [7; 3; -2]; +%! p = [2; 1; 2]; +%! k = [1 0]; +%! e = [2; 1; 1]; +%! [b, a] = residue (r, p, k, e); +%! assert ((abs (b - [1, -5, 9, -3, 1]) < 1e-12 +%! && abs (a - [1, -5, 8, -4]) < 1e-12)); +%! [rr, pr, kr, er] = residue (b, a); +%! [jnk, n] = mpoles(p); +%! assert ((abs (rr - r(n)) < 1e-5 +%! && abs (pr - p(n)) < 1e-7 +%! && abs (kr - k) < 1e-12 +%! && abs (er - e(n)) < 1e-12)); + +