home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / dlyap.m < prev    next >
Text File  |  1999-04-29  |  2KB  |  102 lines
  1. # Copyright (C) 1993, 1994, 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 x = dlyap (a, b)
  20.  
  21. # Usage: x = dlyap (a, b)
  22. #
  23. # Solve a x a' - x + b = 0 (discrete Lyapunov equation) for square
  24. # matrices a and b.  If b is not square, then the function attempts 
  25. # to solve either
  26. #
  27. #  a x a' - x + b b' = 0
  28. #
  29. # or
  30. #
  31. #  a' x a - x + b' b = 0
  32. #
  33. # whichever is appropriate.  Uses Schur decomposition as in Kitagawa
  34. # (1977).
  35.  
  36. # Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993.
  37.  
  38.   if ((n = is_sqr (a)) == 0)
  39.     warning ("dlyap: a must be square");
  40.   endif
  41.  
  42.   if ((m = is_sqr (b)) == 0)
  43.     [n1, m] = size (b);
  44.     if (n1 == n)
  45.       b = b*b';
  46.       m = n1;
  47.     else
  48.       b = b'*b;
  49.       a = a';
  50.     endif
  51.   endif
  52.  
  53.   if (n != m)
  54.     warning ("dlyap: a,b not conformably dimensioned");
  55.   endif
  56.  
  57.   # Solve the equation column by column.
  58.  
  59.   [u, s] = schur (a);
  60.   b = u'*b*u;
  61.  
  62.   j = n;
  63.   while (j > 0)
  64.     j1 = j;
  65.  
  66. # Check for Schur block.
  67.  
  68.     if (j == 1)
  69.       blksiz = 1;
  70.     elseif (s (j, j-1) != 0)
  71.       blksiz = 2;
  72.       j = j - 1;
  73.     else
  74.       blksiz = 1;
  75.     endif
  76.  
  77.     Ajj = kron (s (j:j1, j:j1), s) - eye (blksiz*n);
  78.  
  79.     rhs = reshape (b (:, j:j1), blksiz*n, 1);
  80.  
  81.     if (j1 < n)
  82.       rhs2 = s*(x (:, (j1+1):n) * s (j:j1, (j1+1):n)');
  83.       rhs = rhs + reshape (rhs2, blksiz*n, 1);
  84.     endif
  85.  
  86.     v = - Ajj\rhs;
  87.     x (:, j) = v (1:n);
  88.  
  89.     if(blksiz == 2)
  90.       x (:, j1) = v ((n+1):blksiz*n);
  91.     endif
  92.  
  93.     j = j - 1;
  94.  
  95.   endwhile
  96.  
  97. # Back-transform to original coordinates.
  98.  
  99.   x = u*x*u';
  100.  
  101. endfunction
  102.