home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts.fat / polynom / residue.m < prev    next >
Text File  |  1999-12-24  |  10KB  |  348 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. ## -*- texinfo -*-
  21. ## @deftypefn {Function File} {} residue (@var{b}, @var{a}, @var{tol})
  22. ## If @var{b} and @var{a} are vectors of polynom coefficients, then
  23. ## residue calculates the partial fraction expansion corresponding to the
  24. ## ratio of the two polynoms.
  25. ## @cindex partial fraction expansion
  26. ## 
  27. ## The function @code{residue} returns @var{r}, @var{p}, @var{k}, and
  28. ## @var{e}, where the vector @var{r} contains the residue terms, @var{p}
  29. ## contains the pole values, @var{k} contains the coefficients of a direct
  30. ## polynom term (if it exists) and @var{e} is a vector containing the
  31. ## powers of the denominators in the partial fraction terms.
  32. ## 
  33. ## Assuming @var{b} and @var{a} represent polynoms
  34. ## @iftex
  35. ## @tex
  36. ## $P(s)$ and $Q(s)$
  37. ## @end tex
  38. ## @end iftex
  39. ## @ifinfo
  40. ##  P (s) and Q(s)
  41. ## @end ifinfo
  42. ##  we have:
  43. ## @iftex
  44. ## @tex
  45. ## $$
  46. ## {P(s)\over Q(s)} = \sum_{m=1}^M {r_m\over (s-p_m)^e_m}
  47. ##   + \sum_{i=1}^N k_i s^{N-i}.
  48. ## $$
  49. ## @end tex
  50. ## @end iftex
  51. ## @ifinfo
  52. ## 
  53. ## @example
  54. ##  P(s)    M       r(m)         N
  55. ##  ---- = SUM -------------  + SUM k(i)*s^(N-i)
  56. ##  Q(s)   m=1 (s-p(m))^e(m)    i=1
  57. ## @end example
  58. ## @end ifinfo
  59. ## 
  60. ## @noindent
  61. ## where @var{M} is the number of poles (the length of the @var{r},
  62. ## @var{p}, and @var{e} vectors) and @var{N} is the length of the @var{k}
  63. ## vector.
  64. ## 
  65. ## The argument @var{tol} is optional, and if not specified, a default
  66. ## value of 0.001 is assumed.  The tolerance value is used to determine
  67. ## whether poles with small imaginary components are declared real.  It is
  68. ## also used to determine if two poles are distinct.  If the ratio of the
  69. ## imaginary part of a pole to the real part is less than @var{tol}, the
  70. ## imaginary part is discarded.  If two poles are farther apart than
  71. ## @var{tol} they are distinct.  For example,
  72. ## 
  73. ## @example
  74. ## @group
  75. ##  b = [1, 1, 1];
  76. ##  a = [1, -5, 8, -4];
  77. ##  [r, p, k, e] = residue (b, a);
  78. ##      @result{} r = [-2, 7, 3]
  79. ##      @result{} p = [2, 2, 1]
  80. ##      @result{} k = [](0x0)
  81. ##      @result{} e = [1, 2, 1]
  82. ## @end group
  83. ## @end example
  84. ## 
  85. ## @noindent
  86. ## which implies the following partial fraction expansion
  87. ## @iftex
  88. ## @tex
  89. ## $$
  90. ## {s^2+s+1\over s^3-5s^2+8s-4} = {-2\over s-2} + {7\over (s-2)^2} + {3\over s-1}
  91. ## $$
  92. ## @end tex
  93. ## @end iftex
  94. ## @ifinfo
  95. ## 
  96. ## @example
  97. ##         s^2 + s + 1       -2        7        3
  98. ##    ------------------- = ----- + ------- + -----
  99. ##    s^3 - 5s^2 + 8s - 4   (s-2)   (s-2)^2   (s-1)
  100. ## @end example
  101. ## @end ifinfo
  102. ## @end deftypefn
  103.  
  104. ## SEE ALSO: poly, roots, conv, deconv, polyval, polyderv, polyintg
  105.  
  106. ## Author: Tony Richardson <arichard@stark.cc.oh.us>
  107. ## Created: June 1994
  108. ## Adapted-By: jwe
  109.  
  110. function [r, p, k, e] = residue (b, a, toler)
  111.  
  112.   ## Here's the method used to find the residues.
  113.   ## The partial fraction expansion can be written as:
  114.   ##
  115.   ##
  116.   ##   P(s)    D   M(k)      A(k,m)
  117.   ##   ---- =  #    #    -------------
  118.   ##   Q(s)   k=1  m=1   (s - pr(k))^m
  119.   ##
  120.   ## (# is used to represent a summation) where D is the number of
  121.   ## distinct roots, pr(k) is the kth distinct root, M(k) is the
  122.   ## multiplicity of the root, and A(k,m) is the residue cooresponding
  123.   ## to the kth distinct root with multiplicity m.  For example,
  124.   ##
  125.   ##        s^2            A(1,1)   A(2,1)    A(2,2)
  126.   ## ------------------- = ------ + ------ + -------
  127.   ## s^3 + 4s^2 + 5s + 2    (s+2)    (s+1)   (s+1)^2
  128.   ##
  129.   ## In this case there are two distinct roots (D=2 and pr = [-2 -1]),
  130.   ## the first root has multiplicity one and the second multiplicity
  131.   ## two (M = [1 2]) The residues are actually stored in vector format as
  132.   ## r = [ A(1,1) A(2,1) A(2,2) ].
  133.   ##
  134.   ## We then multiply both sides by Q(s).  Continuing the example:
  135.   ##
  136.   ## s^2 = r(1)*(s+1)^2 + r(2)*(s+1)*(s+2) + r(3)*(s+2)
  137.   ##
  138.   ## or
  139.   ##
  140.   ## s^2 = r(1)*(s^2+2s+1) + r(2)*(s^2+3s+2) +r(3)*(s+2)
  141.   ##
  142.   ## The coefficients of the polynoms on the right are stored in a row
  143.   ## vector called rhs, while the coefficients of the polynom on the
  144.   ## left is stored in a row vector called lhs.  If the multiplicity of
  145.   ## any root is greater than one we'll also need derivatives of this
  146.   ## equation of order up to the maximum multiplicity minus one.  The
  147.   ## derivative coefficients are stored in successive rows of lhs and
  148.   ## rhs.
  149.   ##
  150.   ## For our example lhs and rhs would be:
  151.   ##
  152.   ##       | 1 0 0 |
  153.   ## lhs = |       |
  154.   ##       | 0 2 0 |
  155.   ##
  156.   ##       | 1 2 1 1 3 2 0 1 2 |
  157.   ## rhs = |                   |
  158.   ##       | 0 2 2 0 2 3 0 0 1 |
  159.   ##
  160.   ## We then form a vector B and a matrix A obtained by evaluating the
  161.   ## polynoms in lhs and rhs at the pole values. If a pole has a
  162.   ## multiplicity greater than one we also evaluate the derivative
  163.   ## polynoms (successive rows) at the pole value.
  164.   ##
  165.   ## For our example we would have
  166.   ##
  167.   ## | 4|   | 1 0 0 |   | r(1) |
  168.   ## | 1| = | 0 0 1 | * | r(2) |
  169.   ## |-2|   | 0 1 1 |   | r(3) |
  170.   ##
  171.   ## We then solve for the residues using matrix division.
  172.  
  173.   if (nargin < 2 || nargin > 3)
  174.     usage ("residue (b, a [, toler])");
  175.   endif
  176.  
  177.   if (nargin == 2)
  178.     toler = .001;
  179.   endif
  180.  
  181.   ## Make sure both polynoms are in reduced form.
  182.  
  183.   a = polyred (a);
  184.   b = polyred (b);
  185.  
  186.   b = b / a(1);
  187.   a = a / a(1);
  188.  
  189.   la = length (a);
  190.   lb = length (b);
  191.  
  192.   ## Handle special cases here.
  193.  
  194.   if (la == 0 || lb == 0)
  195.     k = r = p = e = [];
  196.     return;
  197.   elseif (la == 1)
  198.     k = b / a;
  199.     r = p = e = [];
  200.     return;
  201.   endif
  202.  
  203.   ## Find the poles.
  204.  
  205.   p = roots (a);
  206.   lp = length (p);
  207.  
  208.   ## Determine if the poles are (effectively) real.
  209.  
  210.   index = find (abs (imag (p) ./ real (p)) < toler);
  211.   if (length (index) != 0)
  212.     p (index) = real (p (index));
  213.   endif
  214.  
  215.   ## Find the direct term if there is one.
  216.  
  217.   if (lb >= la)
  218.     ## Also returns the reduced numerator.
  219.     [k, b] = deconv (b, a);
  220.     lb = length (b);
  221.   else
  222.     k = [];
  223.   endif
  224.  
  225.   if (lp == 1)
  226.     r = polyval (b, p);
  227.     e = 1;
  228.     return;
  229.   endif
  230.  
  231.  
  232.   ## We need to determine the number and multiplicity of the roots.
  233.   ##
  234.   ## D  is the number of distinct roots.
  235.   ## M  is a vector of length D containing the multiplicity of each root.
  236.   ## pr is a vector of length D containing only the distinct roots.
  237.   ## e  is a vector of length lp which indicates the power in the partial
  238.   ##    fraction expansion of each term in p.
  239.  
  240.   ## Set initial values.  We'll remove elements from pr as we find
  241.   ## multiplicities.  We'll shorten M afterwards.
  242.  
  243.   e = ones (lp, 1);
  244.   M = zeros (lp, 1);
  245.   pr = p;
  246.   D = 1;
  247.   M(1) = 1;
  248.  
  249.   old_p_index = 1;
  250.   new_p_index = 2;
  251.   M_index = 1;
  252.   pr_index = 2;
  253.  
  254.   while (new_p_index <= lp)
  255.     if (abs (p (new_p_index) - p (old_p_index)) < toler)
  256.       ## We've found a multiple pole.
  257.       M (M_index) = M (M_index) + 1;
  258.       e (new_p_index) = e (new_p_index-1) + 1;
  259.       ## Remove the pole from pr.
  260.       pr (pr_index) = [];
  261.     else
  262.       ## It's a different pole.
  263.       D++;
  264.       M_index++;
  265.       M (M_index) = 1;
  266.       old_p_index = new_p_index;
  267.       pr_index++;
  268.     endif
  269.     new_p_index++;
  270.   endwhile
  271.  
  272.   ## Shorten M to it's proper length
  273.  
  274.   M = M (1:D);
  275.  
  276.   ## Now set up the polynom matrices.
  277.  
  278.   MM = max(M);
  279.  
  280.   ## Left hand side polynom
  281.  
  282.   lhs = zeros (MM, lb);
  283.   rhs = zeros (MM, lp*lp);
  284.   lhs (1, :) = b;
  285.   rhi = 1;
  286.   dpi = 1;
  287.   mpi = 1;
  288.   while (dpi <= D)
  289.     for ind = 1:M(dpi)
  290.       if (mpi > 1 && (mpi+ind) <= lp)
  291.         cp = [p(1:mpi-1); p(mpi+ind:lp)];
  292.       elseif (mpi == 1)
  293.         cp = p (mpi+ind:lp);
  294.       else
  295.         cp = p (1:mpi-1);
  296.       endif
  297.       rhs (1, rhi:rhi+lp-1) = prepad (poly (cp), lp);
  298.       rhi = rhi + lp;
  299.     endfor
  300.     mpi = mpi + M (dpi);
  301.     dpi++;
  302.   endwhile
  303.   if (MM > 1)
  304.     for index = 2:MM
  305.       lhs (index, :) = prepad (polyderv (lhs (index-1, :)), lb);
  306.       ind = 1;
  307.       for rhi = 1:lp
  308.         cp = rhs (index-1, ind:ind+lp-1);
  309.         rhs (index, ind:ind+lp-1) = prepad (polyderv (cp), lp);
  310.         ind = ind + lp;
  311.       endfor
  312.     endfor
  313.   endif
  314.  
  315.   ## Now lhs contains the numerator polynom and as many derivatives as
  316.   ## are required.  rhs is a matrix of polynoms, the first row
  317.   ## contains the corresponding polynom for each residue and
  318.   ## successive rows are derivatives.
  319.  
  320.   ## Now we need to evaluate the first row of lhs and rhs at each
  321.   ## distinct pole value.  If there are multiple poles we will also need
  322.   ## to evaluate the derivatives at the pole value also.
  323.  
  324.   B = zeros (lp, 1);
  325.   A = zeros (lp, lp);
  326.  
  327.   dpi = 1;
  328.   row = 1;
  329.   while (dpi <= D)
  330.     for mi = 1:M(dpi)
  331.       B (row) = polyval (lhs (mi, :), pr (dpi));
  332.       ci = 1;
  333.       for col = 1:lp
  334.         cp = rhs (mi, ci:ci+lp-1);
  335.         A (row, col) = polyval (cp, pr(dpi));
  336.         ci = ci + lp;
  337.       endfor
  338.       row++;
  339.     endfor
  340.     dpi++;
  341.   endwhile
  342.  
  343.   ## Solve for the residues.
  344.  
  345.   r = A \ B;
  346.  
  347. endfunction
  348.