lyap

LYAP  Solve continuous-time Lyapunov equations.

   X = LYAP(A,C) solves the special form of the Lyapunov matrix 
   equation:

           A*X + X*A' = -C

   X = LYAP(A,B,C) solves the general form of the Lyapunov matrix
   equation (also called Sylvester equation):

           A*X + X*B = -C

   See also  DLYAP.

0001 function X = lyap(A, B, C)
0002 %LYAP  Solve continuous-time Lyapunov equations.
0003 %
0004 %   X = LYAP(A,C) solves the special form of the Lyapunov matrix
0005 %   equation:
0006 %
0007 %           A*X + X*A' = -C
0008 %
0009 %   X = LYAP(A,B,C) solves the general form of the Lyapunov matrix
0010 %   equation (also called Sylvester equation):
0011 %
0012 %           A*X + X*B = -C
0013 %
0014 %   See also  DLYAP.
0015 
0016 %    S.N. Bangert 1-10-86
0017 %    Copyright (c) 1986-1999 The Mathworks, Inc. All Rights Reserved.
0018 %    $Revision: 1.1.1.1 $  $Date: 2006/08/29 07:25:30 $
0019 %    Last revised JNL 3-24-88, AFP 9-3-95
0020 
0021 ni = nargin;
0022 
0023 if ni==2,
0024    C = B;
0025    B = A';
0026 end
0027 
0028 [ma,na] = size(A);
0029 [mb,nb] = size(B);
0030 [mc,nc] = size(C);
0031 
0032 % A and B must be square and C must have the rows of A and columns of B
0033 if (ma ~= na) | (mb ~= nb) | (mc ~= ma) | (nc ~= mb)
0034    error('Dimensions do not agree.')
0035 elseif ma==0,
0036    X = [];
0037    return
0038 end
0039 
0040 % Perform schur decomposition on A (and convert to complex form)
0041 [ua,ta] = schur(A);
0042 [ua,ta] = rsf2csf(ua,ta);
0043 if ni==2,
0044    % Schur decomposition of A' can be calculated from that of A.
0045    j = ma:-1:1;
0046    ub = ua(:,j);
0047    tb = ta(j,j)';
0048 else
0049    % Perform schur decomposition on B (and convert to complex form)
0050    [ub,tb] = schur(B);
0051    [ub,tb] = rsf2csf(ub,tb);
0052 end
0053 
0054 % Check all combinations of ta(i,i)+tb(j,j) for zero
0055 p1 = diag(ta).'; % Use .' instead of ' in case A and B are not real
0056 p1 = p1(ones(mb,1),:);
0057 p2 = diag(tb);
0058 p2 = p2(:,ones(ma,1));
0059 sum = abs(p1) + abs(p2);
0060 if any(any(sum == 0)) | any(any(abs(p1 + p2) < 1000*eps*sum))
0061    error('Solution does not exist or is not unique.');
0062 end
0063 
0064 % Transform C
0065 ucu = -ua'*C*ub;
0066 
0067 % Solve for first column of transformed solution
0068 y = zeros(ma,mb);
0069 ema = eye(ma);
0070 y(:,1) = (ta+ema*tb(1,1))\ucu(:,1);
0071 
0072 % Solve for remaining columns of transformed solution
0073 for k=2:mb,
0074    km1 = 1:(k-1);
0075    y(:,k) = (ta+ema*tb(k,k))\(ucu(:,k)-y(:,km1)*tb(km1,k));
0076 end
0077 
0078 % Find untransformed solution
0079 X = ua*y*ub';
0080 
0081 % Ignore complex part if real inputs (better be small)
0082 if isreal(A) & isreal(B) & isreal(C),
0083    X = real(X);
0084 end
0085 
0086 % Force X to be symmetric if ni==2 and C is symmetric
0087 if (ni==2) & isequal(C,C'),
0088    X = (X+X')/2;
0089 end
0090 
0091 % end lyap