home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-bin / lib / octave / 1.1.1 / m / control / lyap.m < prev    next >
Encoding:
Text File  |  1996-10-12  |  2.2 KB  |  106 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 = lyap (a, b, c)
  20.  
  21. # Usage: x = lyap (a, b {,c})
  22. #
  23. # If (a, b, c) are specified, then lyap returns the solution of the
  24. # Sylvester equation
  25. #
  26. #   a x + x b + c = 0
  27. #
  28. # If only (a, b) are specified, then lyap returns the solution of the 
  29. # Lyapunov equation
  30. #
  31. #   a' x + x a + b = 0
  32. #
  33. # If b is not square, then lyap returns the solution of either
  34. #
  35. #   a' x + x a + b' b = 0     
  36. #
  37. # or
  38. #
  39. #   a x + x a' + b b' = 0
  40. #
  41. # whichever is appropriate.
  42. #
  43. # Solves by using the Bartels-Stewart algorithm (1972).
  44.  
  45. # Written by A. S. Hodel (scotte@eng.auburn.edu) August 1993.
  46.   
  47.  
  48.   if (nargin != 3 && nargin != 2)
  49.     usage ("lyap (a, b {,c})");
  50.   endif
  51.  
  52.   if ((n = is_square(a)) == 0)
  53.     error ("lyap: a is not square");
  54.   endif
  55.  
  56.   if (nargin == 2)
  57.  
  58. # Transform Lyapunov equation to Sylvester equation form.
  59.  
  60.     if ((m = is_square (b)) == 0)
  61.       if ((m = rows (b)) == n)
  62.  
  63. # solve a x + x a' + b b' = 0
  64.  
  65.     b = b * b';
  66.     a = a';
  67.       else 
  68.  
  69. # Try to solve a'x + x a + b' b = 0.
  70.  
  71.     m = columns (b);
  72.     b = b' * b;
  73.       endif
  74.  
  75.       if (m != n)
  76.     error ("lyap: a, b not conformably dimensioned");
  77.       endif
  78.     endif
  79.  
  80. # Set up Sylvester equation.
  81.  
  82.     c = b;
  83.     b = a;
  84.     a = b'
  85.  
  86.   else 
  87.  
  88. # Check dimensions.
  89.  
  90.     if ((m = is_square (b)) == 0)
  91.       error ("lyap: b must be square in a sylvester equation");
  92.     endif
  93.  
  94.     [n1, m1] = size(c);
  95.  
  96.     if (n != n1 || m != m1)
  97.       error("lyap: a,b,c not conformably dimensioned");
  98.     endif;
  99.   endif
  100.  
  101. # Call octave built-in function.
  102.  
  103.   x = syl (a, b, c);
  104.  
  105. endfunction
  106.