routh.m

基于ＭＡＴＬＡＢ的辅助设计 简述了矩阵分析的重要性

	
function [rtab,np,imagr] = routh(den)
%
% ROUTH  ROUTH is the function which is used to configurate
%        the routh table RTAB of a given system whose charac-
%        teristic polynomial is DEN.  It is also used to
%        calculate the number of positive eigenvalues NP.
%
%         [RTAB, NP, IMAGR] = ROUTH(DEN)
%
%   RTAB - returns the Routh table of the system
%   NP - returns the number of positive poles
%   IMAGR - returns the imaginary poles, if no, return []
%   DEN -  the denominator of the system
%

%    Author:  Ole Barup Sorensen, Rapid Data Ltd 

%    Copyright (c) 1989-94 by Rapid Data Ltd
%    Revision 14:57  26/01/94
	
imagr = []; key = 1;
nc = size(den); den1 = den;
den1(nc(2)+1) = 0;
for i=1:2:nc(2)
   i0 = (i+1)/2;
   vec1(i0) = den1(i);
   vec2(i0) = den1(i+1);
end
iend = nc(2)-1;
irow = 2;
n1 = size(vec1);
rtab = zeros(nc(2),n1(2));
rtab(1,:) = vec1(1,:);
while iend ~= 0
   if sum(abs(vec2)) < eps
      if key == 1,
%         vv = zeros(1,2*iend-1);
%         for ii = 1:iend,ii, vv(2*ii-1) = vec1(ii); end
%         aa = roots(vv);
         aa = roots(den);
         iii = 0; n0 = size(aa); key = 0;
         for ii = 1:n0(1)
             if (abs(real(aa(ii))) <= 1e-10)
                 iii = iii+1; imagr(iii) = aa(ii);
             end
         end
      end
      n0 = iend;
      for i =1:n1(2)
         vec2(i) = n0*vec1(i);
         n0 =n0-2;
      end
   elseif abs(vec2(1)) < eps
      vec2(1) = eps;
   end
   for i =1:n1(2)
      rtab(irow,i) = vec2(i);
   end
   irow = irow +1;
   for i =2:n1(2)
      vec3(i-1) = vec1(i)-vec2(i)*vec1(1)/vec2(1);
   end
   vec3(n1(2)) = 0;
   iend = iend-1;
   vec1 = vec2; vec2 = vec3;
end
np = 0;
for i=1:nc(2)-1
   if rtab(i,1)*rtab(i+1,1) < 0
      np = np+1;
   end
end
for i=1:nc(2)
   for j=1:n1(2)
      if abs(rtab(i,j)) > 1d10
          rtab(i,j) = inf * sign(rtab(i,j));
      end
   end
end