home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / LYAP.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.7 KB  |  75 lines
  1. function X = lyap(A, B, C)
  2. %LYAP    Lyapunov equation.
  3. %    X = LYAP(A,C) solves the special form of the Lyapunov matrix 
  4. %    equation:
  5. %
  6. %        A*X + X*A' = -C
  7. %
  8. %    X = LYAP(A,B,C) solves the general form of the Lyapunov matrix
  9. %    equation:
  10. %
  11. %        A*X + X*B = -C
  12. %
  13. %    See also DLYAP.
  14.  
  15. %    S.N. Bangert 1-10-86
  16. %    Copyright (c) 1986-93 by The MathWorks, Inc.
  17. %    Last revised JNL 3-24-88
  18.  
  19.     if nargin == 2
  20.         C = B;
  21.         B = A';
  22.     end
  23.     [ma,na] = size(A);
  24.     [mb,nb] = size(B);
  25.     [mc,nc] = size(C);
  26.  
  27.     if (ma ~= na) | (mb ~= nb) | (mc ~= ma) | (nc ~= mb)
  28.         error('Dimensions do not agree.')
  29.     elseif ma==0,
  30.         X = []; return
  31.     end
  32.  
  33. % check if problem has any complex inputs
  34.     real_flg = 1;
  35.     if any(any(imag(A))) | any(any(imag(B))) | any(any(imag(C)))
  36.         real_flg = 0;
  37.     end
  38.  
  39. % perform schur decomposition on A and B (note: complex schur form forced by
  40. % adding small complex part so ua and ub are complex upper triangular)
  41.     [ua,ta] = schur(A);
  42.     [ua,ta] = rsf2csf(ua,ta);
  43.     [ub,tb] = schur(B);
  44.     [ub,tb] = rsf2csf(ub,tb);
  45.  
  46. % check all combinations of ua(i,i)+ub(j,j) for zero
  47.     p1 = ones(mb,1)*diag(ta)';
  48.     p2 = diag(tb)*ones(1,ma);
  49.     sum = abs(p1) + abs(p2);
  50.     if any(any(sum == 0)) | any(any(abs(p1 + p2) < 1000*eps*sum))
  51.         error('Solution is not unique.');
  52.     end
  53.  
  54. % transform C
  55.     ucu = -ua'*C*ub;
  56.  
  57. % solve for first column of transformed solution
  58.     y(ma,mb) = 0;
  59.     ema = eye(ma);
  60.     y(:,1) = (ta+ema*tb(1,1))\ucu(:,1);
  61.  
  62. % solve for remaining columns of transformed solution
  63.     for k=2:mb
  64.         km1 = 1:(k-1);
  65.         y(:,k) = (ta+ema*tb(k,k))\(ucu(:,k)-y(:,km1)*tb(km1,k));
  66.     end
  67.  
  68. % find untransformed solution 
  69.     X = ua*y*ub';
  70.  
  71. % ignore complex part if real inputs (better be small)
  72.     if real_flg
  73.         X = real(X);
  74.     end
  75.