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