Mercurial > hg > octave-terminal
changeset 11617:ca0cbc46abce release-3-0-x
[3-0-0-branch @ 2008-01-18 05:47:55 by jwe]
author | jwe |
---|---|
date | Fri, 18 Jan 2008 05:47:56 +0000 |
parents | ad944c3cc888 |
children | 3257cccb9ed7 |
files | scripts/ChangeLog scripts/polynomial/residue.m |
diffstat | 2 files changed, 40 insertions(+), 22 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,8 @@ +2008-01-18 Ben Abbott <bpabbott@mac.com> + + * polynomial/residue.m: For each group of pole multiplicity, set + the poles of the group to the value of the group's average. + 2008-01-17 Tetsuro KURITA <tkurita@mac.com> * plot/print.m: Handle PDF output.
--- a/scripts/polynomial/residue.m +++ b/scripts/polynomial/residue.m @@ -82,7 +82,7 @@ ## @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 +## @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}, and the pole multiplicity @var{e}. ## @@ -191,24 +191,21 @@ p = roots (a); lp = length (p); - ## Determine if the poles are (effectively) zero. - - small = max (abs (p)); - small = max ([small, 1] ) * 1e-8 * (1 + numel (p))^2; - p(abs (p) < small) = 0; - - ## Determine if the poles are (effectively) real, or imaginary. - - index = (abs (imag (p)) < small); - p(index) = real (p(index)); - index = (abs (real (p)) < small); - p(index) = 1i * imag (p(index)); - ## Sort poles so that multiplicity loop will work. [e, indx] = mpoles (p, toler, 1); p = p (indx); + ## For each group of pole multiplicity, set the value of each + ## pole to the average of the group. This reduces the error in + ## the resulting poles. + + p_group = cumsum (e == 1); + for ng = 1:p_group(end) + m = find (p_group == ng); + p(m) = mean (p(m)); + endfor + ## Find the direct term if there is one. if (lb >= la) @@ -219,6 +216,22 @@ k = []; endif + ## Determine if the poles are (effectively) zero. + + small = max (abs (p)); + small = max ([small, 1]) * eps*1e4 * (1 + numel (p))^2; + p(abs (p) < small) = 0; + + ## Determine if the poles are (effectively) real, or imaginary. + + index = (abs (imag (p)) < small); + p(index) = real (p(index)); + index = (abs (real (p)) < small); + p(index) = 1i * imag (p(index)); + + ## The remainder determines the residues. The case of one pole + ## is trivial. + if (lp == 1) r = polyval (b, p); return; @@ -336,8 +349,8 @@ %! b = [1, 1, 1]; %! a = [1, -5, 8, -4]; %! [r, p, k, e] = residue (b, a); -%! assert (abs (r - [-2; 7; 3]) < 1e-5 -%! && abs (p - [2; 2; 1]) < 1e-7 +%! assert (abs (r - [-2; 7; 3]) < 1e-12 +%! && abs (p - [2; 2; 1]) < 1e-12 %! && isempty (k) %! && e == [1; 2; 1]); %! k = [1 0]; @@ -353,10 +366,10 @@ %!test %! b = [1, 0, 1]; %! a = [1, 0, 18, 0, 81]; -%! [r, p, k, e] = residue(b, a); +%! [r, p, k, e] = residue (b, a); %! r1 = [-5i; 12; +5i; 12]/54; %! p1 = [+3i; +3i; -3i; -3i]; -%! assert (abs (r - r1) < 1e-7 && abs (p - p1) < 1e-7 +%! assert (abs (r - r1) < 1e-12 && abs (p - p1) < 1e-12 %! && isempty (k) %! && e == [1; 2; 1; 2]); %! [br, ar] = residue (r, p, k); @@ -373,18 +386,18 @@ %! && 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 +%! assert ((abs (rr - r(n)) < 1e-12 +%! && abs (pr - p(n)) < 1e-12 %! && abs (kr - k) < 1e-12 %! && abs (er - e(n)) < 1e-12)); %!test %! b = [1]; %! a = [1, 10, 25]; -%! [r, p, k, e] = residue(b, a); +%! [r, p, k, e] = residue (b, a); %! r1 = [0; 1]; %! p1 = [-5; -5]; -%! assert (abs (r - r1) < 1e-7 && abs (p - p1) < 1e-7 +%! assert (abs (r - r1) < 1e-12 && abs (p - p1) < 1e-12 %! && isempty (k) %! && e == [1; 2]); %! [br, ar] = residue (r, p, k);