This is residue.m in view mode; [Download] [Up]
## Copyright (C) 1996 John W. Eaton
##
## 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, 59 Temple Place - Suite 330, Boston, MA
## 02111-1307, USA.
## usage: [r, p, k, e] = residue (b, a)
##
## If b and a are vectors of polynomial coefficients, then residue
## calculates the partial fraction expansion corresponding to the
## ratio of the two polynomials. The vector r contains the residue
## terms, p contains the pole values, k contains the coefficients of
## a direct polynomial term (if it exists) and e is a vector containing
## the powers of the denominators in the partial fraction terms.
## Assuming b and a represent polynomials P(s) and Q(s) we have:
##
## P(s) M r(m) N
## ---- = # ------------- + # k(n)*s^(N-n)
## Q(s) m=1 (s-p(m))^e(m) n=1
##
## (# represents summation) where M is the number of poles (the length of
## the r, p, and e vectors) and N is the length of the k vector.
##
## [r p k e] = residue(b,a,tol)
##
## This form of the function call may be used to set a tolerance value.
## The default value is 0.001. 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 tol, the
## imaginary part is discarded. If two poles are farther apart than tol
## they are distinct.
##
## Example:
## b = [1, 1, 1];
## a = [1, -5, 8, -4];
##
## [r, p, k, e] = residue (b, a)
##
## returns
##
## r = [-2, 7, 3]; p = [2, 2, 1]; k = []; e = [1, 2, 1];
##
## which implies the following partial fraction expansion
##
## s^2 + s + 1 -2 7 3
## ------------------- = ----- + ------- + -----
## s^3 - 5s^2 + 8s - 4 (s-2) (s-2)^2 (s-1)
##
## SEE ALSO: poly, roots, conv, deconv, polyval, polyderiv, polyinteg
## Author: Tony Richardson <amr@mpl.ucsd.edu>
## 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.
if (nargin < 2 || nargin > 3)
usage ("residue (b, a [, toler])");
endif
if (nargin == 2)
toler = .001;
endif
## Make sure both polynomials are in reduced form.
a = polyreduce (a);
b = polyreduce (b);
b = b / a(1);
a = 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);
## Determine if the poles are (effectively) real.
index = find (abs (imag (p) ./ real (p)) < toler);
if (length (index) != 0)
p (index) = real (p (index));
endif
## Find the direct term if there is one.
if (lb >= la)
## Also returns the reduced numerator.
[k, b] = deconv (b, a);
lb = length (b);
else
k = [];
endif
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;
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
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);
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);
ind = 1;
for rhi = 1:lp
cp = rhs (index-1, ind:ind+lp-1);
rhs (index, ind:ind+lp-1) = prepad (polyderiv (cp), lp);
ind = ind + lp;
endfor
endfor
endif
## 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
## Solve for the residues.
r = A \ B;
endfunction
These are the contents of the former NiCE NeXT User Group NeXTSTEP/OpenStep software archive, currently hosted by Netfuture.ch.