【matlab】劳斯判据实现

matlab实现劳斯判据表的绘制

考虑全零行、首列有0情况下的劳斯表绘制

% 劳斯判据的matlab实现
% 功能：考虑全零行、首列有0情况下的劳斯表绘制
 
%% 初始化
clear;close all; clc
coeffVector = input('输入系数矩阵：\n [an an-1 an-2 ... a0] = ');

%coeffVector =[1 3 12 20 35 25]; %全0行 示例系数向量
%coeffVector =[1 0 -3 2 ];%首列有0示例系数向量

coeffLength = length(coeffVector); %rh表的行数
rhTableColumn = round(coeffLength/2);  %rh表的列数
rhTable = zeros(coeffLength,rhTableColumn); %不赋值的都是0
rhTable(1,:) = coeffVector(1,1:2:coeffLength);

% 根据系数向量长度的奇偶 给rh表前两行赋值
if (rem(coeffLength,2) ~= 0)
    rhTable(2,1:rhTableColumn-1) = coeffVector(1,2:2:coeffLength);
else
    rhTable(2,:)=coeffVector(1,2:2:coeffLength);
end
%首列有0
if rhTable(2,1) ==0
    rhTable(2,1) = 1e-3;
end
%% 计算rh表
for i = 3:coeffLength
    %全0行
    if abs(rhTable(i-1,:)) <=1e-10 %首项为0，小数计算有时为小量不等于0
        power_first = coeffLength - i+2 ;% 全0行的上一行首项的幂数
        power_decrease= 0;%幂数减少量
        for j = 1:rhTableColumn - 1
            rhTable(i-1,j) = (power_first - power_decrease) * rhTable(i-2,j);
            power_decrease = power_decrease + 2;
        end
    fprintf('\n第%1.0f行为全零行\n',i-1);
    end

    %首列有0
    if rhTable(i,1) ==0
        rhTable(i,1) = 1e-3;
    end

    %非全0行
    for x = 1:rhTableColumn-1
        rhTable(i,x) =(rhTable(i-1,1)* rhTable(i-2,x+1)-...
            (rhTable(i-2,1)*rhTable(i-1,x+1)))/rhTable(i-1,1);
    end
end
%% 判断稳定性及不稳定极点个数
unstablePoles = 0;
for i = 1:coeffLength-1;
    if sign(rhTable(i,1))*sign(rhTable(i+1,1))==-1
        unstablePoles = unstablePoles+1;
    end
end

fprintf('\n劳斯表为：')
rhTable

if unstablePoles == 0
    fprintf('该系统稳定\n')
else
    fprintf('该系统不稳定\n')
end
fprintf('\n虚轴右侧的极点数为 %2.0f \n',unstablePoles)

程序运行截图