Mercurial > octave
view scripts/polynomial/residue.m @ 29358:0a5b15007766 stable
update Octave Project Developers copyright for the new year
In files that have the "Octave Project Developers" copyright notice,
update for 2021.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 10 Feb 2021 09:52:15 -0500 |
parents | 9f9ac219896d |
children | 7854d5752dd2 6dfc06f55cd2 |
line wrap: on
line source
######################################################################## ## ## Copyright (C) 1994-2021 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. ## ## 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 3 of the License, 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, see ## <https://www.gnu.org/licenses/>. ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {[@var{r}, @var{p}, @var{k}, @var{e}] =} residue (@var{b}, @var{a}) ## @deftypefnx {} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}) ## @deftypefnx {} {[@var{b}, @var{a}] =} residue (@var{r}, @var{p}, @var{k}, @var{e}) ## The first calling form computes the partial fraction expansion for the ## quotient of the polynomials, @var{b} and @var{a}. ## ## The quotient is defined as ## @tex ## $$ ## {B(s)\over A(s)} = \sum_{m=1}^M {r_m\over (s-p_m)^e_m} ## + \sum_{i=1}^N k_i s^{N-i}. ## $$ ## @end tex ## @ifnottex ## ## @example ## @group ## B(s) M r(m) N ## ---- = SUM ------------- + SUM k(i)*s^(N-i) ## A(s) m=1 (s-p(m))^e(m) i=1 ## @end group ## @end example ## ## @end ifnottex ## @noindent ## where @math{M} is the number of poles (the length of the @var{r}, @var{p}, ## and @var{e}), the @var{k} vector is a polynomial of order @math{N-1} ## representing the direct contribution, and the @var{e} vector specifies the ## multiplicity of the m-th residue's pole. ## ## For example, ## ## @example ## @group ## b = [1, 1, 1]; ## a = [1, -5, 8, -4]; ## [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 ## ## @noindent ## which represents the following partial fraction expansion ## @tex ## $$ ## {s^2+s+1\over s^3-5s^2+8s-4} = {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1} ## $$ ## @end tex ## @ifnottex ## ## @example ## @group ## s^2 + s + 1 -2 7 3 ## ------------------- = ----- + ------- + ----- ## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1) ## @end group ## @end example ## ## @end ifnottex ## ## The second calling form performs the inverse operation and computes 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}, and the pole ## multiplicity @var{e}. ## ## If the multiplicity, @var{e}, is not explicitly specified the multiplicity ## is determined by the function @code{mpoles}. ## ## For example: ## ## @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] ## ## where mpoles 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 ## ## @noindent ## which represents the following partial fraction expansion ## @tex ## $$ ## {-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} ## $$ ## @end tex ## @ifnottex ## ## @example ## @group ## -2 7 3 s^4 - 5s^3 + 9s^2 - 3s + 1 ## ----- + ------- + ----- + s = -------------------------- ## (s-2) (s-2)^2 (s-1) s^3 - 5s^2 + 8s - 4 ## @end group ## @end example ## ## @end ifnottex ## @seealso{mpoles, poly, roots, conv, deconv} ## @end deftypefn function [r, p, k, e] = residue (b, a, varargin) if (nargin < 2 || nargin > 4) print_usage (); endif tol = .001; 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}, tol, e); return; endif ## Make sure both polynomials are in reduced form, and scaled. a = polyreduce (a); b = polyreduce (b); b /= a(1); a /= a(1); la = length (a); lb = length (b); ## Handle special cases here. if (la == 0 || lb == 0) k = r = p = e = []; return; elseif (la == 1) k = b / a; r = p = e = []; return; endif ## Find the poles. p = roots (a); lp = length (p); ## Sort poles so that multiplicity loop will work. [e, idx] = mpoles (p, tol, 1); p = p(idx); ## 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) ## Also return the reduced numerator. [k, b] = deconv (b, a); lb = length (b); else k = []; endif ## Determine if the poles are (effectively) zero. small = max (abs (p)); if (isa (a, "single") || isa (b, "single")) small = max ([small, 1]) * eps ("single") * 1e4 * (1 + numel (p))^2; else small = max ([small, 1]) * eps * 1e4 * (1 + numel (p))^2; endif p(abs (p) < small) = 0; ## Determine if the poles are (effectively) real, or imaginary. idx = (abs (imag (p)) < small); p(idx) = real (p(idx)); idx = (abs (real (p)) < small); p(idx) = 1i * imag (p(idx)); ## The remainder determines the residues. The case of one pole is trivial. if (lp == 1) r = polyval (b, p); return; endif ## 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. aorder = numel (a) - 1; border = aorder - 1; ## Construct a system of equations relating the individual ## contributions from each residue to the complete numerator. 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, [], tol), border+1, 0).'; endfor ## Solve for the residues. ## FIXME: Use a pre-conditioner d to make A \ B work better (bug #53869). ## It would be better to construct A and B so they are not close to ## singular in the first place. d = max (abs (A), [], 2); r = (diag (d) \ A) \ (B ./ d); endfunction ## Reconstitute the numerator and denominator polynomials ## from the residues, poles, and direct term. function [pnum, pden, e] = rresidue (r, p, k = [], tol = [], e = []) if (! isempty (e)) idx = 1:numel (p); else [e, idx] = mpoles (p, tol, 0); p = p(idx); r = r(idx); endif idx = 1:numel (p); for n = idx 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 multiple 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 = idx(abs (r) > 0) p1 = [1, -p(n)]; pn = 1; for j = 1:n - 1 pn = conv (pn, [1, -p(j)]); end for j = n + 1:numel (p) pn = conv (pn, [1, -p(j)]); end for j = 1:e(n) - 1 pn = deconv (pn, p1); end pn = r(n) * pn; pnum += prepad (pn, N+1, 0, 2); endfor ## Add the direct term. if (numel (k)) pnum += conv (pden, k); endif pnum = polyreduce (pnum); pden = polyreduce (pden); endfunction %!test %! b = [1, 1, 1]; %! a = [1, -5, 8, -4]; %! [r, p, k, e] = residue (b, a); %! assert (r, [-2; 7; 3], 1e-12); %! assert (p, [2; 2; 1], 1e-12); %! assert (isempty (k)); %! assert (e, [1; 2; 1]); %! k = [1 0]; %! b = conv (k, a) + prepad (b, numel (k) + numel (a) - 1, 0); %! a = a; %! [br, ar] = residue (r, p, k); %! assert (br, b, 1e-12); %! assert (ar, a, 1e-12); %! [br, ar] = residue (r, p, k, e); %! assert (br, b, 1e-12); %! assert (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; %! p1 = [+3i; +3i; -3i; -3i]; %! assert (r, r1, 1e-12); %! assert (p, p1, 1e-12); %! assert (isempty (k)); %! assert (e, [1; 2; 1; 2]); %! [br, ar] = residue (r, p, k); %! assert (br, b, 1e-12); %! assert (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 (b, [1, -5, 9, -3, 1], 1e-12); %! assert (a, [1, -5, 8, -4], 1e-12); %! [rr, pr, kr, er] = residue (b, a); %! [jnk, n] = mpoles (p); %! assert (rr, r(n), 1e-12); %! assert (pr, p(n), 1e-12); %! assert (kr, k, 1e-12); %! assert (er, e(n), 1e-12); %!test %! b = [1]; %! a = [1, 10, 25]; %! [r, p, k, e] = residue (b, a); %! r1 = [0; 1]; %! p1 = [-5; -5]; %! assert (r, r1, 1e-12); %! assert (p, p1, 1e-12); %! assert (isempty (k)); %! assert (e, [1; 2]); %! [br, ar] = residue (r, p, k); %! assert (br, b, 1e-12); %! assert (ar, a, 1e-12); ## The following test is due to Bernard Grung %!test <*34266> %! z1 = 7.0372976777e6; %! p1 = -3.1415926536e9; %! p2 = -4.9964813512e8; %! r1 = -(1 + z1/p1)/(1 - p1/p2)/p2/p1; %! r2 = -(1 + z1/p2)/(1 - p2/p1)/p2/p1; %! r3 = (1 + (p2 + p1)/p2/p1*z1)/p2/p1; %! r4 = z1/p2/p1; %! r = [r1; r2; r3; r4]; %! p = [p1; p2; 0; 0]; %! k = []; %! e = [1; 1; 1; 2]; %! b = [1, z1]; %! a = [1, -(p1 + p2), p1*p2, 0, 0]; %! [br, ar] = residue (r, p, k, e); %! assert (br, [0,0,b], 1e-7); %! assert (ar, a, 1e-8); %!test <*49291> %! rf = [1e3, 2e3, 1e3, 2e3]; %! cf = [316.2e-9, 50e-9, 31.6e-9, 5e-9]; %! [num, den] = residue (1./cf,-1./(rf.*cf),0); %! assert (numel (num), 4); %! assert (numel (den), 5); %! assert (den(1), 1); %!test <*51148> %! r = [1.0000e+18, 3.5714e+12, 2.2222e+11, 2.1739e+10]; %! pin = [-1.9231e+15, -1.6234e+09, -4.1152e+07, -1.8116e+06]; %! k = 0; %! [p, q] = residue (r, pin, k); %! assert (p(4), 4.6828e+42, -1e-5);