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