home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / RIC_EIG.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.0 KB  |  69 lines
  1. % function [x1,x2,fail,reig_min] = ric_eig(ham,epp)
  2. %
  3. %  Solves the eigenvalue problem associated with the
  4. %  stabilizing solution (A+R*X stable) of the Riccati equation 
  5. %
  6. %               A'*X + X*A + X*R*X - Q = 0.
  7. %
  8. %  An eigenvalue decomposition is used to obtain the stable
  9. %  invariant subspace of the Hamiltonian matrix, HAM,
  10. %  which contains the above variables in the following format:
  11. %
  12. %       HAM = [A  R; Q  -A'].
  13. %
  14. %  If HAM has no jw axis eigenvalues, there exists n x n matrices
  15. %  X1 and X2 such that [ X1 ; X2 ] spans the n-dimensional
  16. %  stable-invariant subspace of HAM. IF X1 is invertible, then
  17. %  X = X2/X1 satisfies the Riccati equation, and results in A+RX being
  18. %  stable. FAIL is returned with a value of 1 if jw-axis eigenvalues 
  19. %  of HAM are found. An eigenvalue is considered to be purely 
  20. %  imaginary if the magnitude of the real part is less than EPP.
  21. %  The minimum real part of the eigenvalues is returned in REIG_MIN.
  22. %  The default value for EPP is 1e-10.
  23. %
  24. %  This can give incorrect results if HAM is not diagonalizable,
  25. %  in this case it is better to use RIC_SCHR.
  26. %
  27. %  See also: EIG and RIC_SCHR.
  28.  
  29. function [r1,r2,fail,reig_min] = ric_eig(ham,epp)
  30.   if nargin <= 0
  31.     disp('usage: [x1,x2] = ric_eig(ham,epp)')
  32.     return
  33.   end
  34.   [mtype,mrows,mcols,mnum] = minfo(ham);
  35.   if mtype ~= 'cons'
  36.     error(['RIC_EIG is valid only with CONSTANT matrices'])
  37.     return
  38.   end
  39.   if nargin == 1
  40.     epp = 1e-10;
  41.   end
  42.   np = mrows/2;
  43.   fail = 0;
  44.   [v,d] = eig(ham);
  45.   reig_min = min(abs(real(diag(d))));
  46. %
  47. % check there are no jw-axis eigenvalues in the Hamiltonian matrix
  48. %
  49.   if reig_min < epp
  50.     fail = 1;
  51.   else
  52.     v2 = zeros(2*np,np);
  53.     i2 = 1;
  54.     for i=1:2*np
  55.       if real(d(i,i))<0
  56.         v2(:,i2)=v(:,i);
  57.         i2=i2+1;
  58.       end 
  59.     end
  60.     if i2 == np+1
  61.       r1 = v2(1:np,:);
  62.       r2 = v2((np+1):2*np,:);
  63.     else
  64.       fail = 1;
  65.     end
  66.   end
  67. %
  68. % Copyright MUSYN INC 1991,  All Rights Reserved
  69.