Routh.m
上传用户:eighthdate
上传日期:2014-05-24
资源大小:270k
文件大小:2k
源码类别:

其他行业

开发平台:

Matlab

  1. % Function routh(a) constructs the Routh table for a polynomial of degree n
  2. % Hadi Saadat,  1998
  3. function  routh(a)
  4. head=['                        '
  5.       '    Routh-Hurwitz Array '
  6.       '                        '];
  7. n=length(a);
  8. jw=0;
  9. m=2;
  10. nc=round(n/2);
  11. b=zeros(n,nc);
  12. z=zeros(1,nc);
  13. if round(n/2) > n/2
  14. a(n+1)=0;
  15. else,end
  16. for i=1:2:n
  17. k=(i+1)/2;
  18. b(1,k)=a(i);
  19. b(2,k)=a(i+1);
  20. end
  21. if b(2 ,:)==z
  22.    fprintf('Elements of row %g',2)
  23.    fprintf(' are all zero.n')
  24.    fprintf('They are replaced by the auxiliary Eq. coefficients nn')
  25.    jw=1;
  26.    for k=1:nc
  27.    j=n-1;  d=j+2-2*k;
  28.    b(2,k)=d*b(1,k);
  29.    end
  30. else,end
  31. for i=1:n-2
  32.   for j=1:nc-1
  33.      if b(i+1,1)==0
  34.      b(i+1,1)=0.00001;
  35.      fprintf('Zero in the first column is replaced by 0.00001 nn')
  36.      else,end
  37.      b(m+i,j)=(b(i+1,1)*b(i, j+1)-b(i+1,j+1)*b(i,1)) /b(i+1,1);
  38.   end
  39.     if b(m+i,:) == z
  40.      if m+i <n
  41.       fprintf('Elements of row %g',m+i)
  42.       fprintf(' are all zero.n')
  43.       fprintf('They are replaced by the auxiliary Eq. coefficients nn')
  44.       jw=1;
  45.         for k=1:nc
  46.            j=n-m-i+1;
  47.            d=j+2-2*k;
  48.            if d< 0
  49.            b(m+i,k)= b(m+i-1,k);
  50.            else,
  51.            b(m+i,k)=d*b(m+i-1,k);
  52.            end
  53.         end
  54.       else, jw=3; end
  55.     else,end
  56. end
  57. disp(head)
  58. format short e
  59. disp(b)
  60. format short
  61. j=0 ;nr=0;
  62. for i= 1:n
  63.    if j ==0
  64.        if b(i,1) <0 ,nr=nr+1 ;j =1;
  65.        else,end
  66.      else,if b(i,1) > 0 ,nr=nr+1; j=0 ; else, end
  67.    end
  68. end
  69. if jw==3, fprintf('Characterisitc Equation contain root at originn'),end
  70. if jw==1, fprintf('Characteristic Equation include roots on jw-axis or n')
  71.      fprintf('pairs of real or complex roots symmetrical about jw-axis nn'),end
  72. if nr==0
  73.     if jw==0, fprintf('System is stable nn'),end
  74.     else,fprintf('There are %g',nr)
  75. fprintf(' roots in the right half s-plane nn'),end