Mercurial > octave
changeset 31916:b35cd94027b0
maint: merge stable to default
author | Rik <rik@octave.org> |
---|---|
date | Tue, 21 Mar 2023 17:38:55 -0700 |
parents | 034c8ab3d7d3 (current diff) 8eeb49ec16d8 (diff) |
children | 655e757c7522 |
files | |
diffstat | 1 files changed, 80 insertions(+), 42 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/polynomial/mpoles.m Mon Mar 20 22:54:54 2023 -0400 +++ b/scripts/polynomial/mpoles.m Tue Mar 21 17:38:55 2023 -0700 @@ -29,13 +29,21 @@ ## @deftypefnx {} {[@var{multp}, @var{idxp}] =} mpoles (@var{p}, @var{tol}, @var{reorder}) ## Identify unique poles in @var{p} and their associated multiplicity. ## -## The output is ordered from pole with largest magnitude to smallest -## magnitude. +## By default, the output is ordered from the pole with the largest magnitude +## to the smallest magnitude. +## +## Two poles are considered to be multiples if the difference between them +## is less than the relative tolerance @var{tol}. ## -## If the relative difference of two poles is less than @var{tol} then they are -## considered to be multiples. The default value for @var{tol} is 0.001. +## @example +## abs (@var{p1} - @var{p0}) / abs (@var{p0}) < @var{tol} +## @end example ## -## If the optional parameter @var{reorder} is zero, poles are not sorted. +## If the pole is 0 then no scaling is done and @var{tol} is interpreted as an +## absolute tolerance. The default value for @var{tol} is 0.001. +## +## If the optional parameter @var{reorder} is false/zero, poles are not +## sorted. ## ## The output @var{multp} is a vector specifying the multiplicity of the poles. ## @code{@var{multp}(n)} refers to the multiplicity of the Nth pole @@ -56,71 +64,101 @@ ## @seealso{residue, poly, roots, conv, deconv} ## @end deftypefn -function [multp, indx] = mpoles (p, tol, reorder) +function [multp, idxp] = mpoles (p, tol, reorder) if (nargin < 1) print_usage (); endif - if (nargin < 2 || isempty (tol)) - tol = 0.001; - endif + if (! isfloat (p)) + error ("mpoles: P must be a single or double floating point vector"); + endif - if (nargin < 3 || isempty (reorder)) - reorder = true; - endif + if (nargin < 2 || isempty (tol)) + tol = 0.001; + elseif (! (isscalar (tol) && isreal (tol) && tol > 0)) + error ("mpoles: TOL must be a real scalar greater than 0"); + endif + + if (nargin < 3 || isempty (reorder)) + reorder = true; + elseif (! (isscalar (reorder) && isreal (reorder))) + error ("mpoles: REORDER must be a numeric or logical scalar"); + endif Np = numel (p); - - ## force poles to be a column vector - - p = p(:); + p = p(:); # force poles to be a column vector if (reorder) ## sort with largest magnitude first - [~, ordr] = sort (abs (p), "descend"); - p = p(ordr); + [~, order] = sort (abs (p), "descend"); + p = p(order); else - ordr = (1:Np).'; + order = (1:Np).'; endif - ## find pole multiplicity by comparing relative difference of poles + ## Create vector of tolerances for use in algorithm. + vtol = zeros (Np, 1, class (p)); + p_nz = (p != 0); # non-zero poles + vtol(! p_nz) = tol; # use absolute tolerance for zero poles - multp = zeros (Np, 1); - indx = []; + ## Find pole multiplicity by comparing relative difference of poles. + multp = zeros (Np, 1, class (p)); + idxp = []; n = find (multp == 0, 1); while (n) - dp = abs (p-p(n)); - if (p(n) == 0.0) - if (any (abs (p) > 0 & isfinite (p))) - p0 = mean (abs (p(abs (p) > 0 & isfinite (p)))); - else - p0 = 1; - endif - else - p0 = abs (p(n)); - endif - k = find (dp < tol * p0); + dp = abs (p - p(n)); + vtol(p_nz) = tol * abs (p(n)); + k = find (dp < vtol); ## Poles can only be members of one multiplicity group. - if (numel (indx)) - k = k(! ismember (k, indx)); + if (numel (idxp)) + k = k(! ismember (k, idxp)); endif m = 1:numel (k); multp(k) = m; - indx = [indx; k]; + idxp = [idxp; k]; n = find (multp == 0, 1); endwhile - multp = multp(indx); - indx = ordr(indx); + multp = multp(idxp); + idxp = order(idxp); endfunction %!test -%! [mp, n] = mpoles ([0 0], 0.01); +%! [mp, ip] = mpoles ([0 0], 0.01); %! assert (mp, [1; 2]); %!test -%! [mp, n] = mpoles ([-1e4, -0.1, 0]); -%! assert (mp, ones (3, 1)); -%! assert (n, [1; 2; 3]); +%! [mp, ip] = mpoles ([-1e4, -0.1, 0]); +%! assert (mp, [1; 1; 1]); +%! assert (ip, [1; 2; 3]); + +## Test single inputs +%!test +%! [mp, ip] = mpoles (single ([-1e4, -0.1, 0])); +%! assert (mp, single ([1; 1; 1])); +%! assert (ip, [1; 2; 3]); + +## Test relative tolerance criteria +%!test +%! [mp, ip] = mpoles ([1, 1.1, 1.3], .1/1.1); +%! assert (mp, [1; 1; 1]); +%! [mp, ip] = mpoles ([1, 1.1, 1.3], .1/1.1 + eps); +%! assert (mp, [1; 1; 2]); + +## Test absolute tolerance criteria with a zero pole +%!test +%! [mp, ip] = mpoles ([0, -0.1, 0.3], .1); +%! assert (mp, [1; 1; 1]); +%! [mp, ip] = mpoles ([0, -0.1, 0.3], .1 + eps); +%! assert (mp, [1; 1; 2]); + +## Test input validation +%!error <Invalid call> mpoles () +%!error <P must be a single or double floating point vector> mpoles (uint8 (1)) +%!error <TOL must be a real scalar greater than 0> mpoles (1, [1, 2]) +%!error <TOL must be a real scalar greater than 0> mpoles (1, 1i) +%!error <TOL must be a real scalar greater than 0> mpoles (1, 0) +%!error <REORDER must be a numeric or logical scalar> mpoles (1, 1, [1, 2]) +%!error <REORDER must be a numeric or logical scalar> mpoles (1, 1, {1})