clear;
%%%
%
% Consider the following 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 real-valued matrices. 
%
% This function computes the mean-square stability matrix
% of that linear system as well as the mean-square stability matrix
% of Theta-Maruyama and Theta-Milstein approximations. 
% If the spectral abscissa of the stability matrices is less than 0, 
% the zero solution of the system is asympotically mean-square stable.
% If the spectral radius of the stability matrices is less than 1, 
% the Theta-Maruyama and Theta-Milstein approximation
% is asympotically mean-square stable.
%
% Calls: kron, eig, max, inv 
%
% For more details see: 
% E. Buckwar, T. Sickenberger: A structural analysis of asymptotic mean-square  
% stability for multi-dimensional linear stochastic differential systems.
% Heriot-Watt Mathematics Report HWM10-21, 2010.
% Link: http://www.ma.hw.ac.uk/maths/deptreps/HWM10-21.html
%
% 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/
%
%%%

% Select a test equation or define your own linear system
% 1: two-dim lin system with one noise term
% 2: two-dim lin system with two commutative noise terms
% 3: two-dim lin system with two non-commutative noise terms
% 0: your own lin system (see below)
test_equation = 1;

% Define a step-size sequence, which you like to test
h = [ 1, 0.5, 0.2, 0.1, 0.05, 0.02, 0.01, 0.005, 0.002, 0.001 ];

% Define a sequence for the method parameter Theta
% theta = 0:   drift-explicit discretization
% theta = 0.5: trapazoidal rule
% theta = 1:   drift-implizit discretization
theta = [ 0, 0.5, 1];

% Use a test-equation or implement your own test-equation (see else part)
% === test-equation (28) ===
if (test_equation == 1)
  % Two-dimensional system with one noise term
  % Define some parameters for the drift and diffusion matrices
  lambda = -3; sigma = 0.5; epsilon = sqrt(3);

  % Define drift matrix F
  F = [ lambda, 0;
        0, lambda ];
    
  % Define diffusion matrices G_r (r=1,...,m)
  G{1} = [ sigma, epsilon;
           epsilon, sigma ];

% === test-equation (29) ===
elseif (test_equation == 2)
  % Two-dimensional system with two commutative noise terms
  % Define some parameters for the drift and diffusion matrices
  lambda = -2; sigma = 0.5; epsilon = sqrt(3);

  % Define drift matrix F
  F = [ lambda, 0 ;
        0, lambda ];
    
  % Define diffusion matrices G_r (r=1,...,m)
  G{1} = [ sigma, 0;
           0, sigma ];
  G{2} = [ 0, epsilon;
           -epsilon, 0 ];

% === test-equation (30) ===
elseif (test_equation == 3)
  % Two-dimensional system with two non-commutative noise terms
  % Define some parameters for the drift and diffusion matrices
  lambda = -2; sigma = 0.5; epsilon = sqrt(3);

  % Define drift matrix F
  F = [ lambda, 0 ;
        0, lambda ];
    
  % Define diffusion matrices G_r (r=1,...,m)
  G{1} = [ - sigma, 0;
           0, sigma ];    
  G{2} = [ 0, epsilon;
           epsilon, 0 ];
         
% === define your own linear system ===       
else     
  % Implement your own system of linear SDEs
  % F and G_r must be d x d matrices of the same dimension
  % Define drift matrix F
  F = [ -3,  0,  0 ;
         0, -3,  0 ;
         0,  0, -3 ];
    
  % Define diffusion matrices G_r (r=1,...,m)
  G{1} = [ 0, 1, 0 ;
           0, 0, 0 ;
           0, 0, 0 ];

  G{2} = [ 0, 0, 0 ;
           0, 0, 1 ; 
           0, 0, 0 ];

  G{3} = [ 0, 0, 0 ;
           0, 0, 0 ; 
           1, 0, 0 ];
end

% Analysing the system 
com = stabanalysis(F,G);

% Compute the stability coefficients
% for different step-sizes and Theta values
for it=1:length(theta)
  for ih=1:length(h)
      [sMar(ih,it),sMil(ih,it)] = stabcalculation(F,G,h(ih),theta(it),com);
  end
end

% Print results ...
fprintf('\n And here are the results...');
fprintf('\n The Theta-Maruyama approximation is stable (o) or NOT stable (x)');
fprintf('\n step-size ');
for it = 1:length(theta)  fprintf('| theta = %2.4f    ',theta(it)); end
fprintf('\n =====================================================================');
for ih = 1:length(h)
  fprintf('\n  h= %2.3f ',h(ih));
  for it = 1:length(theta)
    fprintf('| spRad = %2.4f ',sMar(ih,it));
    if (sMar(ih,it) <1)  fprintf(' o ')
    else fprintf(' x ');
    end
  end
end

fprintf('\n\n The Theta-Milstein approximation is stable (o) or NOT stable (x)');
fprintf('\n step-size ');
for it = 1:length(theta)  fprintf('| theta = %2.4f    ',theta(it)); end
fprintf('\n =====================================================================');
for ih = 1:length(h)
  fprintf('\n  h= %2.3f ',h(ih));
  for it = 1:length(theta)
    fprintf('| spRad = %2.4f ',sMil(ih,it));
    if (sMil(ih,it) <1)  fprintf(' o ')
    else fprintf(' x ');
    end
  end
end

fprintf('\n');  

