matlab 劳斯稳定判据routh函数，同时实现对特殊情况的处理

matlab 劳斯稳定判据routh函数，同时实现对特殊情况的处理

首先，先明确两种特殊情况如何处理：

1.某行第一列元素为0，而该行元素不全为零

这时，计算劳斯表下一行第一个元素时，会出现分母是0的情况，将0用一个很小的正数代替即可，matlab里eps可以表示为一个很小的正数

2.出现全零行

此时就需要列辅助方程了，至于如何找辅助方程的方法就很简单了，不详细说了。这里具体看如何用代码实现

 m=r(i-1,:);    % m为全零行上面一行元素构成的数组
 s=n-i+1;       % s为全零行上面一行的最高次数
for k=1:len1    %len1为劳斯表列数 
     r(i,k)=s*m(k); 
     if(s-2>0)
     	s=s-2;
	 else
     	s=0;
     end
end

举例下图

则变量s=2，m=[5 25 0]

routh函数代码如下：
r为返回的矩阵形式的劳斯表，info为简要的系统稳定性分析结果

function [r,info]=routh(den)
info='';
n=length(den);
a1=den(1:2:end);
a2=den(2:2:end);
len1=length(a1);
len2=length(a2);
if(len1>len2)
    a2=[a2 0];
end
r=zeros(n,len1);
r(1,:)=a1;
r(2,:)=a2;

for i=3:n
    for j=1:ceil((n-i+1)/2)
        r(i,j)=(r(i-1,1)*r(i-2,j+1)-r(i-1,j+1)*r(i-2,1))/r(i-1,1);
    end
    if(r(i,1)==0&&sum(r(i,2:len1))~=0)
        r(i,1)=eps;
        info='某行第一列为0，而该行不全为零，将0用eps代替';
    end
    if(sum(r(i,:)==0)==len1)
        info='出现全零行';
        m=r(i-1,:);
        s=n-i+1;
        for k=1:len1
            r(i,k)=s*m(k);
            if(s-2>0)
                s=s-2;
            else
                s=0;
            end
        end
    end
end
if(sum(r(:,1)>0)==n)
    info=[info,'系统稳定'];
else
    info=[info,'系统不稳定'];
end