function [com] = stabanalysis(F,G)
%%%
%
% This function analysis the dimension and the noise structur
% of the linear system of SDEs
% dX(t) = F*X(t) dt + \sum_{r=1}^m G_r *X(t)*dW_r(t)
% where F and G_r are dxd-dimensional constant matrices.
%
% The mean-square stability matrix of that linear system is computed,
% cf. Equation (9) and the asympotical mean-square stability of its
% zero solution is checked in terms of the spectral abscissa, cf. Lemma 2.
%
% Input: - F is the drift matrix,
% - G is an array of diffusion matrices,
% Output: - com states if the systems has commutative noises terms (com=1)
% or non-commutative noise terms (com=0).
%
% For more details see the main function
% main_stability.m
%
% Copyright (C) 2010 T Sickenberger
% (Sent comments to t.sickenberger@hw.ac.uk)
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by the
% Free Software Foundation; either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%
% See the GNU General Public License for more details.
% http://www.gnu.org/licenses/
%
%%%
fprintf('\n Start analysing your system ...\n')
% Get the number of noise sources ...
m = length(G);
% First, check the dimension of the drift matrix ...
dimF = size(F);
if ( dimF(1) ~= dimF(2) )
error(' ERROR - Your drift matrix is not quadratic');
end
% Get the dimension ...
d = dimF(1);
% Secondly, check the dimension of all diffusion matrices ...
for r = 1:m
dimG{r} = size(G{r});
if ( dimG{r}(1) ~= dimG{r}(2) )
error(' ERROR - Diffusion matrix is not quadratic');
end
if ( dimG{r}(1) ~= d )
error(' ERROR - Diffusion matrix has not the same dimension as your drift matrix');
end
end
fprintf(' You have entered a %i dimensional system with %i multiplicative noise sources.\n',d,m);
% Check if the noise terms commute
com = true;
if (m > 1)
for r1 = 1:m
for r2 = 1:m
comTest = G{r1}*G{r2} - G{r2}*G{r1};
if (norm(comTest) > 0)
com = false;
fprintf(' ... the diffusion matrices G(%i) and G(%i) do not commute \n',r1,r2);
end
end
end
if (com == true)
fprintf(' Your system has pairwise commutative noise terms. \n');
else
fprintf(' Your system has non-commutative noise terms! \n',r1,r2);
end
end
fprintf('\n Compute asymptotic mean-square stability of the zero-solution of your system ...');
% Compute stabilty matrix, cf. Equation (9)
S = kron(F,eye(d)) + kron(eye(d),F);
for r = 1:m
S = S + kron(G{r},G{r});
end
% Compute spectral abscissa of that stability matrix, cf. Lemma 2.
scoef = max( real( eig(S)));
fprintf('\n The zero-solution of your system is');
if (scoef < 0)
fprintf(' asymptotically mean-square stable');
else
fprintf(' NOT asymptotically mean-square stable');
end
fprintf('\n (with spectral abscissa %2.5f).\n',scoef);