Hurw.m

% Hurw.m
% Hurwitz stability criterion
% from characteristic polynomial
%
% As an input, enter the vector of arbitrary members,
% for example [1 5 4 7 2 1], which is 5th order characteristic polynomial

% Author: Karel Perutka

clear all

C=input('Enter the coefficients of characteristic polynomial, for instance [1 5 4 7 2 1]:');
if length(C)<3
    disp('You have entered the polynomial with its order less than two. Stopping.')
elseif length(C)>=3
    znamenka=sign(C);
    rozhodni=(znamenka==1);
    delka=length(C);
    vektorJednicek=ones(1,delka);
    vektorNul=zeros(1,delka);
    if (isequal(rozhodni,vektorJednicek) || isequal(rozhodni,vektorNul))
        if isequal(rozhodni,vektorNul)
            C=-1*C;
        end
        disp('Testing using Hurwitz criterion is starting, please wait ...')
        prvniRadek=zeros(1,delka-1); 
        druhyRadek=prvniRadek;% faster than direct function calling
        HurwMatice=zeros(delka-1);
        for i=2:2:delka
            prvniRadek(i/2)=C(i);
            druhyRadek(i/2)=C(i-1);
        end
        if rem(delka,2)==1 % for odd number of coefficients
            druhyRadek(ceil(delka/2))=C(delka);
        end
        %naplnn Hurwitzovy matice
        for i=2:2:delka-1
            HurwMatice(i-1,:)=prvniRadek;
            HurwMatice(i,:)=druhyRadek;
            %rotation:
            prvniRadek=[0 prvniRadek];
            druhyRadek=[0 druhyRadek];         
            prvniRadek(end)=[];
            druhyRadek(end)=[];
        end 
        if rem(delka,2)==0 % for even number of coefficients
            HurwMatice(end,:)=prvniRadek;
        end
        HurwMatrix=HurwMatice %displays Hurwitz matrix
        %computation of main minors:
        delkaHMatice=length(HurwMatice);
        for i=1:delkaHMatice
            hlavniMinory(i)=det(HurwMatice(1:i,1:i));
        end
        mainMinors=hlavniMinory %displays main minors of Hurwitz matrix
        pomocnyHM=ones(1,delkaHMatice);
        if isequal(sign(hlavniMinory),pomocnyHM)
            disp('Stability test is finished, system is stable.')
        else
            disp('Stability test is finished, system is unstable.')
        end
    else
        disp('Necessary but insufficient condition does not hold on, program is terminated.')
    end
end