home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / lyap.m < prev    next >
Text File  |  1999-12-24  |  3KB  |  143 lines
  1. ## Copyright (C) 1996, 1997 Auburn University.  All Rights Reserved
  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} {} lyap (@var{a}, @var{b}, @var{c})
  22. ## @deftypefnx {Function File} {} lyap (@var{a}, @var{b})
  23. ##   Solve the Lyapunov (or Sylvester) equation via the Bartels-Stewart
  24. ##   algorithm (Communications of the ACM, 1972).
  25. ## 
  26. ##   If @var{a}, @var{b}, and @var{c} are specified, then @code{lyap} returns
  27. ##   the solution of the  Sylvester equation
  28. ##   @iftex
  29. ##     @tex
  30. ##       $$ A X + X B + C = 0 $$
  31. ##     @end tex
  32. ##   @end iftex
  33. ##   @ifinfo
  34. ##     @example
  35. ##       a x + x b + c = 0
  36. ##     @end example
  37. ##   @end ifinfo
  38. ##   If only @code{(a, b)} are specified, then @code{lyap} returns the
  39. ##   solution of the Lyapunov equation
  40. ##   @iftex
  41. ##     @tex
  42. ##       $$ A^T X + X A + B = 0 $$
  43. ##     @end tex
  44. ##   @end iftex
  45. ##   @ifinfo
  46. ##     @example
  47. ##       a' x + x a + b = 0
  48. ##     @end example
  49. ##   @end ifinfo
  50. ##   If @var{b} is not square, then @code{lyap} returns the solution of either
  51. ##   @iftex
  52. ##     @tex
  53. ##       $$ A^T X + X A + B^T B = 0 $$
  54. ##     @end tex
  55. ##   @end iftex
  56. ##   @ifinfo
  57. ##     @example
  58. ##       a' x + x a + b' b = 0
  59. ##     @end example
  60. ##   @end ifinfo
  61. ##   @noindent
  62. ##   or
  63. ##   @iftex
  64. ##     @tex
  65. ##       $$ A X + X A^T + B B^T = 0 $$
  66. ##     @end tex
  67. ##   @end iftex
  68. ##   @ifinfo
  69. ##     @example
  70. ##       a x + x a' + b b' = 0
  71. ##     @end example
  72. ##   @end ifinfo
  73. ##   @noindent
  74. ##   whichever is appropriate.
  75. ##
  76. ## Solves by using the Bartels-Stewart algorithm (1972).
  77. ## @end deftypefn
  78.  
  79. ## Author: A. S. Hodel <scotte@eng.auburn.edu>
  80. ## Created: August 1993
  81. ## Adapted-By: jwe
  82.  
  83. function x = lyap (a, b, c)
  84.  
  85.   if (nargin != 3 && nargin != 2)
  86.     usage ("lyap (a, b {,c})");
  87.   endif
  88.  
  89.   if ((n = is_sqr(a)) == 0)
  90.     error ("lyap: a is not square");
  91.   endif
  92.  
  93.   if (nargin == 2)
  94.  
  95.     ## Transform Lyapunov equation to Sylvester equation form.
  96.  
  97.     if ((m = is_sqr (b)) == 0)
  98.       if ((m = rows (b)) == n)
  99.  
  100.     ## solve a x + x a' + b b' = 0
  101.  
  102.     b = b * b';
  103.     a = a';
  104.       else
  105.  
  106.     ## Try to solve a'x + x a + b' b = 0.
  107.  
  108.     m = columns (b);
  109.     b = b' * b;
  110.       endif
  111.  
  112.       if (m != n)
  113.     error ("lyap: a, b not conformably dimensioned");
  114.       endif
  115.     endif
  116.  
  117.     ## Set up Sylvester equation.
  118.  
  119.     c = b;
  120.     b = a;
  121.     a = b';
  122.  
  123.   else
  124.  
  125.     ## Check dimensions.
  126.  
  127.     if ((m = is_sqr (b)) == 0)
  128.       error ("lyap: b must be square in a sylvester equation");
  129.     endif
  130.  
  131.     [n1, m1] = size(c);
  132.  
  133.     if (n != n1 || m != m1)
  134.       error("lyap: a,b,c not conformably dimensioned");
  135.     endif
  136.   endif
  137.  
  138.   ## Call octave built-in function.
  139.  
  140.   x = syl (a, b, c);
  141.  
  142. endfunction
  143.