# HG changeset patch # User jwe # Date 1191645078 0 # Node ID 33f20a41aeea0a07725d6b4c05c4aedca057536c # Parent 642f481d2d50a19323c23fafd5023ab029757187 [project @ 2007-10-06 04:31:18 by jwe] diff -r 642f481d2d50 -r 33f20a41aeea scripts/ChangeLog --- a/scripts/ChangeLog Fri Oct 05 22:40:59 2007 +0000 +++ b/scripts/ChangeLog Sat Oct 06 04:31:18 2007 +0000 @@ -1,3 +1,13 @@ +2007-10-05 Ben Abbott + + * polynomial/mpoles.m: New function. + * polynomial/residue.m: Modified to behave in reciprocal + manner. No longer compute 4th output, "e". No longer accept + tolerance input. Explicitly set tolerance parameter to 0.001. + Respect maximum relative difference in poles when determining + their multiplicity. Use mpoles to determine the multiplicity of + poles. + 2007-10-05 Peter A. Gustafson * plot/__go_draw_axes__.m: Add cbrange to the plot stream diff -r 642f481d2d50 -r 33f20a41aeea scripts/polynomial/Makefile.in --- a/scripts/polynomial/Makefile.in Fri Oct 05 22:40:59 2007 +0000 +++ b/scripts/polynomial/Makefile.in Sat Oct 06 04:31:18 2007 +0000 @@ -20,7 +20,7 @@ INSTALL_PROGRAM = @INSTALL_PROGRAM@ INSTALL_DATA = @INSTALL_DATA@ -SOURCES = compan.m conv.m deconv.m mkpp.m pchip.m poly.m \ +SOURCES = compan.m conv.m deconv.m mkpp.m mpoles.m pchip.m poly.m \ polyder.m polyderiv.m polyfit.m polygcd.m polyint.m \ polyout.m polyreduce.m polyval.m polyvalm.m ppval.m residue.m \ roots.m spline.m unmkpp.m diff -r 642f481d2d50 -r 33f20a41aeea scripts/polynomial/mpoles.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/polynomial/mpoles.m Sat Oct 06 04:31:18 2007 +0000 @@ -0,0 +1,112 @@ +## Copyright (C) 2007 Ben Abbott +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{multp}, @var{indx}] =} mpoles (@var{p}) +## @deftypefnx {Function File} {[@var{multp}, @var{indx}] =} mpoles (@var{p}, @var{tol}) +## @deftypefnx {Function File} {[@var{multp}, @var{indx}] =} mpoles (@var{p}, @var{tol}, @var{reorder}) +## Identifiy unique poles in @var{p} and associates their multiplicity, +## ordering them from largest to smallest. +## +## If the relative difference of the poles is less than @var{tol}, then +## they are considered to be multiples. The default value for @var{tol} +## is 0.001. +## +## If the optional parameter @var{reorder} is zero, poles are not sorted. +## +## The value @var{multp} is a vector specifying the multiplicity of the +## poles. @var{multp}(:) refers to mulitplicity of @var{p}(@var{indx}(:)). +## +## For example, +## +## @example +## @group +## p = [2 3 1 1 2]; +## [m, n] = mpoles(p); +## @result{} m = [1; 1; 2; 1; 2] +## @result{} n = [2; 5; 1; 4; 3] +## @result{} p(n) = [3, 2, 2, 1, 1] +## @end group +## @end example +## +## @seealso{poly, roots, conv, deconv, polyval, polyderiv, polyinteg, residue} +## @end deftypefn + +## Author: Ben Abbott +## Created: Sept 30, 2007 + +function [multp, indx] = mpoles (p, tol, reorder) + + if (nargin < 1 || nargin > 3) + print_usage (); + endif + + if (nargin < 2 || isempty (tol)) + tol = 0.001; + endif + + if (nargin < 3 || isempty (reorder)) + reorder = true; + endif + + Np = numel (p); + + ## Force the poles to be a column vector. + + p = p(:); + + ## Sort the poles according to their magnitidues, largest first. + + if (reorder) + ## Sort with smallest magnitude first. + [p, ordr] = sort (p); + ## Reverse order, largest maginitude first. + n = Np:-1:1; + p = p(n); + ordr = ordr(n); + else + ordr = 1:Np; + endif + + ## Find pole multiplicty by comparing the relative differnce in the + ## poles. + + multp = zeros (Np, 1); + indx = []; + n = find (multp == 0, 1); + while (n) + dp = abs (p-p(n)); + if (p(n) == 0.0) + p0 = mean (abs (p(find (abs (p) > 0)))); + if (isempty (p0)) + p0 = 1; + end + else + p0 = abs (p(n)); + endif + k = find (dp < tol * p0); + m = 1:numel (k); + multp(k) = m; + indx = [indx; k]; + n = find (multp == 0, 1); + endwhile + multp = multp(indx); + indx = indx(ordr); + +endfunction diff -r 642f481d2d50 -r 33f20a41aeea scripts/polynomial/residue.m --- a/scripts/polynomial/residue.m Fri Oct 05 22:40:59 2007 +0000 +++ b/scripts/polynomial/residue.m Sat Oct 06 04:31:18 2007 +0000 @@ -1,4 +1,5 @@ ## Copyright (C) 1996, 1997 John W. Eaton +## Copyright (C) 2007 Ben Abbott ## ## This file is part of Octave. ## @@ -18,17 +19,22 @@ ## 02110-1301, USA. ## -*- texinfo -*- -## @deftypefn {Function File} {} residue (@var{b}, @var{a}, @var{tol}) -## If @var{b} and @var{a} are vectors of polynomial coefficients, then -## residue calculates the partial fraction expansion corresponding to the -## ratio of the two polynomials. +## @deftypefn {Function File} {[@var{r}, @var{p}, @var{k}] =} residue (@var{b}, @var{a}) +## @deftypefnx {Function File} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}) +## In the instance of two inputs, they are assumed to correspond to vectors +## of polynomial coefficients, @var{b} and @var{a}. From these polynomial +## coefficients, the function residue calculates the partial fraction +## expansion corresponding to the ratio of the two polynomials. ## @cindex partial fraction expansion ## -## The function @code{residue} returns @var{r}, @var{p}, @var{k}, and -## @var{e}, where the vector @var{r} contains the residue terms, @var{p} -## contains the pole values, @var{k} contains the coefficients of a direct -## polynomial term (if it exists) and @var{e} is a vector containing the -## powers of the denominators in the partial fraction terms. +## In this instance, the function @code{residue} returns @var{r}, @var{p}, +## and @var{k}, where the vector @var{r} contains the residue terms, +## @var{p} contains the pole values, @var{k} contains the coefficients of a +## direct polynomial term (if it exists). +## +## In the instance of three inputs, the function 'residue' performs the +## reciprocal task. Meaning the partial fraction expansion is reconstituted +## into the corresponding pair of polynomial coefficients. ## ## Assuming @var{b} and @var{a} represent polynomials ## @iftex @@ -60,25 +66,19 @@ ## @noindent ## where @math{M} is the number of poles (the length of the @var{r}, ## @var{p}, and @var{e} vectors) and @math{N} is the length of the -## @var{k} vector. +## @var{k} vector. The @var{e} vector specifies the multiplicity of the +## mth residue's pole. ## -## The argument @var{tol} is optional, and if not specified, a default -## value of 0.001 is assumed. The tolerance value is used to determine -## whether poles with small imaginary components are declared real. It is -## also used to determine if two poles are distinct. If the ratio of the -## imaginary part of a pole to the real part is less than @var{tol}, the -## imaginary part is discarded. If two poles are farther apart than -## @var{tol} they are distinct. For example, +## For example, ## ## @example ## @group -## b = [1, 1, 1]; -## a = [1, -5, 8, -4]; -## [r, p, k, e] = residue (b, a); +## 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] ## @result{} k = [](0x0) -## @result{} e = [1, 2, 1] ## @end group ## @end example ## @@ -98,84 +98,61 @@ ## ------------------- = ----- + ------- + ----- ## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1) ## @end example +## ## @end ifinfo -## @seealso{poly, roots, conv, deconv, polyval, polyderiv, and polyinteg} +## A similar, but reciprocal example, where the fraction's polynomials are +## reconstituted from the residues, poles, and direct term is +## +## @example +## @group +## 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] +## @end group +## @end example +## +## @noindent +## which implies the following partial fraction expansion +## @iftex +## @tex +## $$ +## {s^4-5s^3+9s^2-3s+1\over s^3-5s^2+8s-4} = {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} + s +## $$ +## @end tex +## @end iftex +## @ifinfo +## +## @example +## s^4 - 5s^3 + 9s^2 - 3s + 1 -2 7 3 +## -------------------------- = ----- + ------- + ----- + s +## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1) +## @end example +## @end ifinfo +## @seealso{poly, roots, conv, deconv, mpoles, polyval, polyderiv, polyinteg} ## @end deftypefn ## Author: Tony Richardson +## Author: Ben Abbott ## Created: June 1994 ## Adapted-By: jwe -function [r, p, k, e] = residue (b, a, toler) - - ## Here's the method used to find the residues. - ## The partial fraction expansion can be written as: - ## - ## - ## P(s) D M(k) A(k,m) - ## ---- = # # ------------- - ## Q(s) k=1 m=1 (s - pr(k))^m - ## - ## (# is used to represent a summation) where D is the number of - ## distinct roots, pr(k) is the kth distinct root, M(k) is the - ## multiplicity of the root, and A(k,m) is the residue cooresponding - ## to the kth distinct root with multiplicity m. For example, - ## - ## s^2 A(1,1) A(2,1) A(2,2) - ## ------------------- = ------ + ------ + ------- - ## s^3 + 4s^2 + 5s + 2 (s+2) (s+1) (s+1)^2 - ## - ## In this case there are two distinct roots (D=2 and pr = [-2 -1]), - ## the first root has multiplicity one and the second multiplicity - ## two (M = [1 2]) The residues are actually stored in vector format as - ## r = [ A(1,1) A(2,1) A(2,2) ]. - ## - ## We then multiply both sides by Q(s). Continuing the example: - ## - ## s^2 = r(1)*(s+1)^2 + r(2)*(s+1)*(s+2) + r(3)*(s+2) - ## - ## or - ## - ## s^2 = r(1)*(s^2+2s+1) + r(2)*(s^2+3s+2) +r(3)*(s+2) - ## - ## The coefficients of the polynomials on the right are stored in a row - ## vector called rhs, while the coefficients of the polynomial on the - ## left is stored in a row vector called lhs. If the multiplicity of - ## any root is greater than one we'll also need derivatives of this - ## equation of order up to the maximum multiplicity minus one. The - ## derivative coefficients are stored in successive rows of lhs and - ## rhs. - ## - ## For our example lhs and rhs would be: - ## - ## | 1 0 0 | - ## lhs = | | - ## | 0 2 0 | - ## - ## | 1 2 1 1 3 2 0 1 2 | - ## rhs = | | - ## | 0 2 2 0 2 3 0 0 1 | - ## - ## We then form a vector B and a matrix A obtained by evaluating the - ## polynomials in lhs and rhs at the pole values. If a pole has a - ## multiplicity greater than one we also evaluate the derivative - ## polynomials (successive rows) at the pole value. - ## - ## For our example we would have - ## - ## | 4| | 1 0 0 | | r(1) | - ## | 1| = | 0 0 1 | * | r(2) | - ## |-2| | 0 1 1 | | r(3) | - ## - ## We then solve for the residues using matrix division. +function [r, p, k] = residue (b, a, varargin) if (nargin < 2 || nargin > 3) print_usage (); endif - if (nargin == 2) - toler = .001; - endif + toler = .001; + + if (nargin == 3) + ## 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); + return + end ## Make sure both polynomials are in reduced form. @@ -206,17 +183,26 @@ ## Determine if the poles are (effectively) zero. - p(abs (p) < toler) = 0; + 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. - ## Determine if the poles are (effectively) real. + index = (abs (imag (p)) < small); + p(index) = real (p(index)); + index = (abs (real (p)) < small); + p(index) = 1i * imag (p(index)); - index = (abs (p) >= toler & (abs (imag (p)) ./ abs (p) < toler)); - p(index) = real (p(index)); + ## Sort poles so that multiplicity loop will work. + + [e, indx] = mpoles (p, toler, 1); + p = p (indx); ## Find the direct term if there is one. if (lb >= la) - ## Also returns the reduced numerator. + ## Also return the reduced numerator. [k, b] = deconv (b, a); lb = length (b); else @@ -225,123 +211,106 @@ if (lp == 1) r = polyval (b, p); - e = 1; return; endif - ## We need to determine the number and multiplicity of the roots. - ## - ## D is the number of distinct roots. - ## M is a vector of length D containing the multiplicity of each root. - ## pr is a vector of length D containing only the distinct roots. - ## e is a vector of length lp which indicates the power in the partial - ## fraction expansion of each term in p. - - ## Set initial values. We'll remove elements from pr as we find - ## multiplicities. We'll shorten M afterwards. - - e = ones (lp, 1); - M = zeros (lp, 1); - pr = p; - D = 1; - M(1) = 1; - - old_p_index = 1; - new_p_index = 2; - M_index = 1; - pr_index = 2; + ## Determine the order of the denominator and remaining numerator. + ## With the direct term removed the potential order of the numerator + ## is one less than the order of the denominator. - while (new_p_index <= lp) - if (abs (p (new_p_index) - p (old_p_index)) < toler) - ## We've found a multiple pole. - M (M_index) = M (M_index) + 1; - e (new_p_index) = e (new_p_index-1) + 1; - ## Remove the pole from pr. - pr (pr_index) = []; - else - ## It's a different pole. - D++; - M_index++; - M (M_index) = 1; - old_p_index = new_p_index; - pr_index++; - endif - new_p_index++; - endwhile - - ## Shorten M to it's proper length - - M = M (1:D); - - ## Now set up the polynomial matrices. - - MM = max(M); - - ## Left hand side polynomial + aorder = numel (a) - 1; + border = aorder - 1; - lhs = zeros (MM, lb); - rhs = zeros (MM, lp*lp); - lhs (1, :) = b; - rhi = 1; - dpi = 1; - mpi = 1; - while (dpi <= D) - for ind = 1:M(dpi) - if (mpi > 1 && (mpi+ind) <= lp) - cp = [p(1:mpi-1); p(mpi+ind:lp)]; - elseif (mpi == 1) - cp = p (mpi+ind:lp); - else - cp = p (1:mpi-1); - endif - rhs (1, rhi:rhi+lp-1) = prepad (poly (cp), lp, 0, 2); - rhi = rhi + lp; - endfor - mpi = mpi + M (dpi); - dpi++; - endwhile - if (MM > 1) - for index = 2:MM - lhs (index, :) = prepad (polyderiv (lhs (index-1, :)), lb, 0, 2); - ind = 1; - for rhi = 1:lp - cp = rhs (index-1, ind:ind+lp-1); - rhs (index, ind:ind+lp-1) = prepad (polyderiv (cp), lp, 0, 2); - ind = ind + lp; - endfor - endfor - endif + ## Construct a system of equations relating the individual + ## contributions from each residue to the complete numerator. - ## Now lhs contains the numerator polynomial and as many derivatives as - ## are required. rhs is a matrix of polynomials, the first row - ## contains the corresponding polynomial for each residue and - ## successive rows are derivatives. - - ## Now we need to evaluate the first row of lhs and rhs at each - ## distinct pole value. If there are multiple poles we will also need - ## to evaluate the derivatives at the pole value also. - - B = zeros (lp, 1); - A = zeros (lp, lp); - - dpi = 1; - row = 1; - while (dpi <= D) - for mi = 1:M(dpi) - B (row) = polyval (lhs (mi, :), pr (dpi)); - ci = 1; - for col = 1:lp - cp = rhs (mi, ci:ci+lp-1); - A (row, col) = polyval (cp, pr(dpi)); - ci = ci + lp; - endfor - row++; - endfor - dpi++; - endwhile + A = zeros (border+1, border+1); + B = prepad (reshape (b, [numel(b), 1]), border+1, 0); + for ip = 1:numel(p) + ri = zeros (size (p)); + ri(ip) = 1; + A(:,ip) = prepad (rresidue (ri, p, [], toler), border+1, 0).'; + endfor ## Solve for the residues. r = A \ B; endfunction + +function [pnum, pden] = rresidue (r, p, k, toler) + + ## Reconstitute the numerator and denominator polynomials from the + ## residues, poles, and direct term. + + if (nargin < 2 || nargin > 4) + print_usage (); + endif + + if (nargin < 4) + toler = []; + endif + + if (nargin < 3) + k = []; + endif + + [multp, indx] = mpoles (p, toler, 0); + + p = p (indx); + r = r (indx); + + indx = 1:numel(p); + + for n = indx + pn = [1, -p(n)]; + if n == 1 + pden = pn; + else + pden = conv (pden, pn); + endif + endfor + + ## D is the order of the denominator + ## K is the order of the direct polynomial + ## N is the order of the resulting numerator + ## pnum(1:(N+1)) is the numerator's polynomial + ## pden(1:(D+1)) is the denominator's polynomial + ## pm is the multible pole for the nth residue + ## pn is the numerator contribution for the nth residue + + D = numel (pden) - 1; + K = numel (k) - 1; + N = K + D; + pnum = zeros (1, N+1); + for n = indx(abs(r)>0) + p1 = [1, -p(n)]; + for m = 1:multp(n) + if m == 1 + pm = p1; + else + pm = conv (pm, p1); + endif + endfor + pn = deconv (pden, pm); + pn = r(n) * pn; + pnum = pnum + prepad ( pn, N+1, 0); + endfor + + ## Add the direct term. + + if (numel (k)) + pnum = pnum + conv (pden, k); + endif + + ## Check for leading zeros and trim the polynomial coefficients. + + small = max ([max(abs(pden)), max(abs(pnum)), 1]) * eps; + + pnum (abs (pnum) < small) = 0; + pden (abs (pden) < small) = 0; + + pnum = polyreduce (pnum); + pden = polyreduce (pden); + +endfunction