1## Copyright (C) 2009-2016   Lukas F. Reichlin
2##
3## This file is part of LTI Syncope.
4##
5## LTI Syncope is free software: you can redistribute it and/or modify
6## it under the terms of the GNU General Public License as published by
7## the Free Software Foundation, either version 3 of the License, or
8## (at your option) any later version.
9##
10## LTI Syncope is distributed in the hope that it will be useful,
11## but WITHOUT ANY WARRANTY; without even the implied warranty of
12## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13## GNU General Public License for more details.
14##
15## You should have received a copy of the GNU General Public License
16## along with LTI Syncope.  If not, see <http://www.gnu.org/licenses/>.
17
18## -*- texinfo -*-
19## @deftypefn{Function File} {@var{x} =} lyap (@var{a}, @var{b})
20## @deftypefnx{Function File} {@var{x} =} lyap (@var{a}, @var{b}, @var{c})
21## @deftypefnx{Function File} {@var{x} =} lyap (@var{a}, @var{b}, @var{[]}, @var{e})
22## Solve continuous-time Lyapunov or Sylvester equations.
23##
24## @strong{Equations}
25## @example
26## @group
27## AX + XA' + B = 0      (Lyapunov Equation)
28##
29## AX + XB + C = 0       (Sylvester Equation)
30##
31## AXE' + EXA' + B = 0   (Generalized Lyapunov Equation)
32## @end group
33## @end example
34##
35## @strong{Algorithm}@*
36## Uses SLICOT SB03MD, SB04MD and SG03AD by courtesy of
37## @uref{http://www.slicot.org, NICONET e.V.}
38##
39## @seealso{lyapchol, dlyap, dlyapchol}
40## @end deftypefn
41
42## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
43## Created: January 2010
44## Version: 0.2.1
45
46function [x, scale] = lyap (a, b, c, e)
47
48  scale = 1;
49
50  switch (nargin)
51    case 2                                      # Lyapunov equation
52
53      if (! is_real_square_matrix (a, b))
54        ## error ("lyap: a, b must be real and square");
55        error ("lyap: %s, %s must be real and square", ...
56                inputname (1), inputname (2));
57      endif
58
59      if (rows (a) != rows (b))
60        ## error ("lyap: a, b must have the same number of rows");
61        error ("lyap: %s, %s must have the same number of rows", ...
62                inputname (1), inputname (2));
63
64      endif
65
66      [x, scale] = __sl_sb03md__ (a, -b, false);     # AX + XA' = -B
67
68      ## x /= scale;                            # 0 < scale <= 1
69
70    case 3                                      # Sylvester equation
71
72      if (! is_real_square_matrix (a, b))
73        ## error ("lyap: a, b must be real and square");
74        error ("lyap: %s, %s must be real and square", ...
75                inputname (1), inputname (2));
76      endif
77
78      if (! is_real_matrix (c) || rows (c) != rows (a) || columns (c) != columns (b))
79        ## error ("lyap: c must be a real (%dx%d) matrix", rows (a), columns (b));
80        error ("lyap: %s must be a real (%dx%d) matrix", ...
81                rows (a), columns (b), inputname (3));
82      endif
83
84      x = __sl_sb04md__ (a, b, -c);  # AX + XB = -C
85
86    case 4                                      # generalized Lyapunov equation
87
88      if (! isempty (c))
89        print_usage ();
90      endif
91
92      if (! is_real_square_matrix (a, b, e))
93        ## error ("lyap: a, b, e must be real and square");
94        error ("lyap: %s, %s, %s must be real and square", ...
95                inputname (1), inputname (2), inputname (4));
96      endif
97
98      if (rows (b) != rows (a) || rows (e) != rows (a))
99        ## error ("lyap: a, b, e must have the same number of rows");
100        error ("lyap: %s, %s, %s must have the same number of rows", ...
101                inputname (1), inputname (2), inputname (4));
102      endif
103
104      if (! issymmetric (b))
105        ## error ("lyap: b must be symmetric");
106        error ("lyap: %s must be symmetric", ...
107                inputname (2));
108      endif
109
110      [x, scale] = __sl_sg03ad__ (a, e, -b, false);  # AXE' + EXA' = -B
111
112      ## x /= scale;                            # 0 < scale <= 1
113
114    otherwise
115      print_usage ();
116
117  endswitch
118
119  if (scale < 1)
120    warning ("lyap: solution scaled by %g to prevent overflow", scale);
121  endif
122
123endfunction
124
125
126## Lyapunov
127%!shared X, X_exp
128%! A = [1, 2; -3, -4];
129%! Q = [3, 1; 1, 1];
130%! X = lyap (A, Q);
131%! X_exp = [ 6.1667, -3.8333;
132%!          -3.8333,  3.0000];
133%!assert (X, X_exp, 1e-4);
134
135## Sylvester
136%!shared X, X_exp
137%! A = [2.0   1.0   3.0
138%!      0.0   2.0   1.0
139%!      6.0   1.0   2.0];
140%!
141%! B = [2.0   1.0
142%!      1.0   6.0];
143%!
144%! C = [2.0   1.0
145%!      1.0   4.0
146%!      0.0   5.0];
147%!
148%! X = lyap (A, B, -C);
149%!
150%! X_exp = [-2.7685   0.5498
151%!          -1.0531   0.6865
152%!           4.5257  -0.4389];
153%!
154%!assert (X, X_exp, 1e-4);
155
156## Generalized Lyapunov
157%!shared X, X_exp
158%! A = [  3.0     1.0     1.0
159%!        1.0     3.0     0.0
160%!        1.0     0.0     2.0];
161%!
162%! E = [  1.0     3.0     0.0
163%!        3.0     2.0     1.0
164%!        1.0     0.0     1.0];
165%!
166%! B = [-64.0   -73.0   -28.0
167%!      -73.0   -70.0   -25.0
168%!      -28.0   -25.0   -18.0];
169%!
170%! X = lyap (A.', -B, [], E.');
171%!
172%! X_exp = [-2.0000  -1.0000   0.0000
173%!          -1.0000  -3.0000  -1.0000
174%!           0.0000  -1.0000  -3.0000];
175%!
176%!assert (X, X_exp, 1e-4);
177