function [inital,routh_list,conclusion] = routhtable(chara_equ)
% ---------------------------------------------------------
% 输入:
% chara_equ = 特征方程向量
% 输出:
% inital=直接计算时的劳斯表
% routh_list = 劳斯表
% conclusion = 给出系统是否稳定或存在多少个不稳定的根的结论
% ----------------------------------------------------------
n=length(chara_equ);
chara_equ=reshape(chara_equ,1,n);
if mod(n,2)==0
n1=n/2;
else
n1=(n+1)/2;
chara_equ=[chara_equ,0];
end
routh=reshape(chara_equ,2,n1);
routh_list=zeros(n,n1);
routh_list(1:2,:)=routh;i=3;s=0;
disp('直接求解时的劳斯表如下面的第一个表所示。');
D=chara_equ;L=length(D);
if rem(L,2)==1
k=L+1;
else
k=L;
end
Q=zeros(L,k/2);
if rem(L,2)==0
for j=1:k/2
Q(1,j)=D(2*j-1);Q(2,j)=D(2*j);
end
else
for j=1:(k/2)-1
Q(1,j)=D(2*j-1);Q(1,k/2)=D(k-1);Q(2,j)=D(2*j);
end
end
for x=3:L
for j=1:(k/2)-1
Q(x,j)=(Q(x-1,1)*Q(x-2,j+1)-Q(x-2,1)*Q(x-1,j+1))/Q(x-1,1);
end
end
inital=Q;
while 1;
% 特殊情况1(第一列中有0元素,但其所在行不为0)
if routh_list(i-1,1)==0 && sum(routh_list(i-1,2:n1))~=0
chara_equ = conv(chara_equ,[1 3]);
s=1;n=length(chara_equ);
if mod(n,2)==0
n1=n/2;
else
n1=(n+1)/2; chara_equ=[chara_equ,0];
end
routh=reshape(chara_equ,2,n1);
routh_list=zeros(n,n1);
routh_list(1:2,:)=routh;i=3;
end
% 计算劳斯表中元素值
ai=routh_list(i-2,1)/routh_list(i-1,1);
for j=1:n1-1
routh_list(i,j)=routh_list(i-2,j+1)-ai*routh_list(i-1,j+1);
end
% 特殊情况2(有全0行)
if sum(routh_list(i,:))==0
s=2; k=0;l=1; F=zeros(1,n1);
while n-i-k>=0
F(l)=n-i+1-k;
k=k+2;l=l+1;
end
routh_list(i,:)=routh_list(i-1,:).*F(1,:);
end
i=i+1;
if i>n
break;
end
end
% 结果输出
if s==1
disp('在特征方程基础上乘上因式(s+3),重新开始计算劳斯表中的元素值。');
end
if s==2
disp('对全0行上面一行的因子F(s)求导,重新计算劳斯表中的元素值。');
end
r=find(routh_list(:,1)<0);
if isempty(r)==1
conclusion='系统稳定!';
else
n2=length(r);m=n2;
for i=1:n2-1
if r(i+1)-r(i)==1
m=m-1;
end
end
str1='系统不稳定!存在';
if r(n2)==n
str2=num2str(m*2-1);
else
str2=num2str(m*2);
end
str3=' 个具有正实部的极点!';
conclusion = [str1,str2,str3];
end