home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / SZEROS.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.8 KB  |  66 lines
  1. % function out = szeros(sys,epp)
  2. %
  3. %   Finds the transmission zeros of a SYSTEM matrix. 
  4. %   Occasionally, large zeros are included which may 
  5. %   actually be at infinity. Solving for the transmission 
  6. %   zeros of a system involves two generalized eigenvalue 
  7. %   problems. EPP (optional) defines if the difference between 
  8. %   two generalize eigenvalues is small.
  9. %
  10. %   See also: EIG and SPOLES. 
  11.  
  12. %   Algorithm based on Laub & Moore 1978 paper, Automatica
  13.  
  14. function out = szeros(sys,epp)
  15.   if nargin < 1
  16.    disp('usage: szeros(sys)')
  17.    return
  18.   end
  19.   if nargin == 1
  20.    epp = eps;
  21.   end
  22.   [mtype,mrows,mcols,mnum] = minfo(sys);
  23.   if mtype == 'syst'
  24.     [a,b,c,d] = unpck(sys);
  25.     if mnum == 0
  26.       disp('SYSTEM has no states')
  27.     end
  28.     sysu = [a b; c d];
  29.  
  30. % find generalized eigenvalues of a square system matrix
  31.  
  32.     if mrows == mcols
  33.       x = zeros(mnum+mcols,mnum+mcols);
  34.       x(1:mnum,1:mnum) = eye(mnum);
  35.       z = eig(sysu,x);
  36.       out = z(~isnan(z) & finite(z));
  37.     else
  38.       nrm = norm(sysu,1);
  39.       if mrows > mcols
  40.         x1 = [ sysu  nrm*(rand(mnum+mrows,mrows-mcols)-.5)];
  41.         x2 = [ sysu  nrm*(rand(mnum+mrows,mrows-mcols)-.5)];
  42.       else
  43.         x1 = [ sysu; nrm*(rand(mcols-mrows,mnum+mcols)-.5)];
  44.         x2 = [ sysu; nrm*(rand(mcols-mrows,mnum+mcols)-.5)];
  45.       end
  46.  
  47.       x = zeros(x1);
  48.       x(1:mnum,1:mnum) = eye(mnum);
  49.       z1 = eig(x1,x);
  50.       z2 = eig(x2,x);
  51.       z1 = z1(~isnan(z1) & finite(z1));
  52.       z2 = z2(~isnan(z2) & finite(z2));
  53.       nz = length(z1);
  54.       for i=1:nz
  55.         if any(abs(z1(i)-z2) < nrm*sqrt(epp))
  56.           out = [out; z1(i)];
  57.         end
  58.       end
  59.     end
  60.   else
  61.      error('matrix is not a SYSTEM matrix')
  62.      return
  63.   end 
  64. %
  65. % Copyright MUSYN INC 1991,  All Rights Reserved
  66.