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);

