home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts / polynomial / residue.m < prev    next >
Text File  |  1998-11-10  |  9KB  |  307 lines
  1. ## Copyright (C) 1996, 1997 John W. Eaton
  2. ##
  3. ## This file is part of Octave.
  4. ##
  5. ## Octave is free software; you can redistribute it and/or modify it
  6. ## under the terms of the GNU General Public License as published by
  7. ## the Free Software Foundation; either version 2, or (at your option)
  8. ## any later version.
  9. ##
  10. ## Octave is distributed in the hope that it will be useful, but
  11. ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12. ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13. ## General Public License for more details.
  14. ##
  15. ## You should have received a copy of the GNU General Public License
  16. ## along with Octave; see the file COPYING.  If not, write to the Free
  17. ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
  18. ## 02111-1307, USA.
  19.  
  20. ## usage: [r, p, k, e] = residue (b, a)
  21. ##
  22. ## If b and a are vectors of polynomial coefficients, then residue
  23. ## calculates the partial fraction expansion corresponding to the
  24. ## ratio of the two polynomials. The vector r contains the residue
  25. ## terms, p contains the pole values, k contains the coefficients of
  26. ## a direct polynomial term (if it exists) and e is a vector containing
  27. ## the powers of the denominators in the partial fraction terms.
  28. ## Assuming b and a represent polynomials P(s) and Q(s) we have:
  29. ##
  30. ##  P(s)    M       r(m)         N
  31. ##  ---- =  #  -------------  +  # k(n)*s^(N-n)
  32. ##  Q(s)   m=1 (s-p(m))^e(m)    n=1
  33. ##
  34. ## (# represents summation) where M is the number of poles (the length of
  35. ## the r, p, and e vectors) and N is the length of the k vector.
  36. ##
  37. ## [r p k e] = residue(b,a,tol)
  38. ##
  39. ## This form of the function call may be used to set a tolerance value.
  40. ## The default value is 0.001. The tolerance value is used to determine
  41. ## whether poles with small imaginary components are declared real. It is
  42. ## also used to determine if two poles are distinct. If the ratio of the
  43. ## imaginary part of a pole to the real part is less than tol, the
  44. ## imaginary part is discarded. If two poles are farther apart than tol
  45. ## they are distinct.
  46. ##
  47. ## Example:
  48. ##  b = [1,  1, 1];
  49. ##  a = [1, -5, 8, -4];
  50. ##
  51. ##  [r, p, k, e] = residue (b, a)
  52. ##
  53. ##  returns
  54. ##
  55. ##  r = [-2, 7, 3]; p = [2, 2, 1]; k = []; e = [1, 2, 1];
  56. ##
  57. ##  which implies the following partial fraction expansion
  58. ##
  59. ##        s^2 + s + 1         -2       7        3
  60. ##    ------------------- = ----- + ------- + -----
  61. ##    s^3 - 5s^2 + 8s - 4   (s-2)   (s-2)^2   (s-1)
  62. ##
  63. ## SEE ALSO: poly, roots, conv, deconv, polyval, polyderiv, polyinteg
  64.  
  65. ## Author: Tony Richardson <arichard@stark.cc.oh.us>
  66. ## Created: June 1994
  67. ## Adapted-By: jwe
  68.  
  69. function [r, p, k, e] = residue (b, a, toler)
  70.  
  71.   ## Here's the method used to find the residues.
  72.   ## The partial fraction expansion can be written as:
  73.   ##
  74.   ##
  75.   ##   P(s)    D   M(k)      A(k,m)
  76.   ##   ---- =  #    #    -------------
  77.   ##   Q(s)   k=1  m=1   (s - pr(k))^m
  78.   ##
  79.   ## (# is used to represent a summation) where D is the number of
  80.   ## distinct roots, pr(k) is the kth distinct root, M(k) is the
  81.   ## multiplicity of the root, and A(k,m) is the residue cooresponding
  82.   ## to the kth distinct root with multiplicity m.  For example,
  83.   ##
  84.   ##        s^2            A(1,1)   A(2,1)    A(2,2)
  85.   ## ------------------- = ------ + ------ + -------
  86.   ## s^3 + 4s^2 + 5s + 2    (s+2)    (s+1)   (s+1)^2
  87.   ##
  88.   ## In this case there are two distinct roots (D=2 and pr = [-2 -1]),
  89.   ## the first root has multiplicity one and the second multiplicity
  90.   ## two (M = [1 2]) The residues are actually stored in vector format as
  91.   ## r = [ A(1,1) A(2,1) A(2,2) ].
  92.   ##
  93.   ## We then multiply both sides by Q(s).  Continuing the example:
  94.   ##
  95.   ## s^2 = r(1)*(s+1)^2 + r(2)*(s+1)*(s+2) + r(3)*(s+2)
  96.   ##
  97.   ## or
  98.   ##
  99.   ## s^2 = r(1)*(s^2+2s+1) + r(2)*(s^2+3s+2) +r(3)*(s+2)
  100.   ##
  101.   ## The coefficients of the polynomials on the right are stored in a row
  102.   ## vector called rhs, while the coefficients of the polynomial on the
  103.   ## left is stored in a row vector called lhs.  If the multiplicity of
  104.   ## any root is greater than one we'll also need derivatives of this
  105.   ## equation of order up to the maximum multiplicity minus one.  The
  106.   ## derivative coefficients are stored in successive rows of lhs and
  107.   ## rhs.
  108.   ##
  109.   ## For our example lhs and rhs would be:
  110.   ##
  111.   ##       | 1 0 0 |
  112.   ## lhs = |       |
  113.   ##       | 0 2 0 |
  114.   ##
  115.   ##       | 1 2 1 1 3 2 0 1 2 |
  116.   ## rhs = |                   |
  117.   ##       | 0 2 2 0 2 3 0 0 1 |
  118.   ##
  119.   ## We then form a vector B and a matrix A obtained by evaluating the
  120.   ## polynomials in lhs and rhs at the pole values. If a pole has a
  121.   ## multiplicity greater than one we also evaluate the derivative
  122.   ## polynomials (successive rows) at the pole value.
  123.   ##
  124.   ## For our example we would have
  125.   ##
  126.   ## | 4|   | 1 0 0 |   | r(1) |
  127.   ## | 1| = | 0 0 1 | * | r(2) |
  128.   ## |-2|   | 0 1 1 |   | r(3) |
  129.   ##
  130.   ## We then solve for the residues using matrix division.
  131.  
  132.   if (nargin < 2 || nargin > 3)
  133.     usage ("residue (b, a [, toler])");
  134.   endif
  135.  
  136.   if (nargin == 2)
  137.     toler = .001;
  138.   endif
  139.  
  140.   ## Make sure both polynomials are in reduced form.
  141.  
  142.   a = polyreduce (a);
  143.   b = polyreduce (b);
  144.  
  145.   b = b / a(1);
  146.   a = a / a(1);
  147.  
  148.   la = length (a);
  149.   lb = length (b);
  150.  
  151.   ## Handle special cases here.
  152.  
  153.   if (la == 0 || lb == 0)
  154.     k = r = p = e = [];
  155.     return;
  156.   elseif (la == 1)
  157.     k = b / a;
  158.     r = p = e = [];
  159.     return;
  160.   endif
  161.  
  162.   ## Find the poles.
  163.  
  164.   p = roots (a);
  165.   lp = length (p);
  166.  
  167.   ## Determine if the poles are (effectively) real.
  168.  
  169.   index = find (abs (imag (p) ./ real (p)) < toler);
  170.   if (length (index) != 0)
  171.     p (index) = real (p (index));
  172.   endif
  173.  
  174.   ## Find the direct term if there is one.
  175.  
  176.   if (lb >= la)
  177.     ## Also returns the reduced numerator.
  178.     [k, b] = deconv (b, a);
  179.     lb = length (b);
  180.   else
  181.     k = [];
  182.   endif
  183.  
  184.   if (lp == 1)
  185.     r = polyval (b, p);
  186.     e = 1;
  187.     return;
  188.   endif
  189.  
  190.  
  191.   ## We need to determine the number and multiplicity of the roots.
  192.   ##
  193.   ## D  is the number of distinct roots.
  194.   ## M  is a vector of length D containing the multiplicity of each root.
  195.   ## pr is a vector of length D containing only the distinct roots.
  196.   ## e  is a vector of length lp which indicates the power in the partial
  197.   ##    fraction expansion of each term in p.
  198.  
  199.   ## Set initial values.  We'll remove elements from pr as we find
  200.   ## multiplicities.  We'll shorten M afterwards.
  201.  
  202.   e = ones (lp, 1);
  203.   M = zeros (lp, 1);
  204.   pr = p;
  205.   D = 1;
  206.   M(1) = 1;
  207.  
  208.   old_p_index = 1;
  209.   new_p_index = 2;
  210.   M_index = 1;
  211.   pr_index = 2;
  212.  
  213.   while (new_p_index <= lp)
  214.     if (abs (p (new_p_index) - p (old_p_index)) < toler)
  215.       ## We've found a multiple pole.
  216.       M (M_index) = M (M_index) + 1;
  217.       e (new_p_index) = e (new_p_index-1) + 1;
  218.       ## Remove the pole from pr.
  219.       pr (pr_index) = [];
  220.     else
  221.       ## It's a different pole.
  222.       D++;
  223.       M_index++;
  224.       M (M_index) = 1;
  225.       old_p_index = new_p_index;
  226.       pr_index++;
  227.     endif
  228.     new_p_index++;
  229.   endwhile
  230.  
  231.   ## Shorten M to it's proper length
  232.  
  233.   M = M (1:D);
  234.  
  235.   ## Now set up the polynomial matrices.
  236.  
  237.   MM = max(M);
  238.  
  239.   ## Left hand side polynomial
  240.  
  241.   lhs = zeros (MM, lb);
  242.   rhs = zeros (MM, lp*lp);
  243.   lhs (1, :) = b;
  244.   rhi = 1;
  245.   dpi = 1;
  246.   mpi = 1;
  247.   while (dpi <= D)
  248.     for ind = 1:M(dpi)
  249.       if (mpi > 1 && (mpi+ind) <= lp)
  250.         cp = [p(1:mpi-1); p(mpi+ind:lp)];
  251.       elseif (mpi == 1)
  252.         cp = p (mpi+ind:lp);
  253.       else
  254.         cp = p (1:mpi-1);
  255.       endif
  256.       rhs (1, rhi:rhi+lp-1) = prepad (poly (cp), lp);
  257.       rhi = rhi + lp;
  258.     endfor
  259.     mpi = mpi + M (dpi);
  260.     dpi++;
  261.   endwhile
  262.   if (MM > 1)
  263.     for index = 2:MM
  264.       lhs (index, :) = prepad (polyderiv (lhs (index-1, :)), lb);
  265.       ind = 1;
  266.       for rhi = 1:lp
  267.         cp = rhs (index-1, ind:ind+lp-1);
  268.         rhs (index, ind:ind+lp-1) = prepad (polyderiv (cp), lp);
  269.         ind = ind + lp;
  270.       endfor
  271.     endfor
  272.   endif
  273.  
  274.   ## Now lhs contains the numerator polynomial and as many derivatives as
  275.   ## are required.  rhs is a matrix of polynomials, the first row
  276.   ## contains the corresponding polynomial for each residue and
  277.   ## successive rows are derivatives.
  278.  
  279.   ## Now we need to evaluate the first row of lhs and rhs at each
  280.   ## distinct pole value.  If there are multiple poles we will also need
  281.   ## to evaluate the derivatives at the pole value also.
  282.  
  283.   B = zeros (lp, 1);
  284.   A = zeros (lp, lp);
  285.  
  286.   dpi = 1;
  287.   row = 1;
  288.   while (dpi <= D)
  289.     for mi = 1:M(dpi)
  290.       B (row) = polyval (lhs (mi, :), pr (dpi));
  291.       ci = 1;
  292.       for col = 1:lp
  293.         cp = rhs (mi, ci:ci+lp-1);
  294.         A (row, col) = polyval (cp, pr(dpi));
  295.         ci = ci + lp;
  296.       endfor
  297.       row++;
  298.     endfor
  299.     dpi++;
  300.   endwhile
  301.  
  302.   ## Solve for the residues.
  303.  
  304.   r = A \ B;
  305.  
  306. endfunction
  307.