function [sMar, sMil] = stabcalculation(F,G,h,theta,com)
%%%
%
% This function computes the mean-square stability matrix for 
% Theta-Maruyama and Theta-Milstein approximation 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.
%
% An approximation is asymptotically mean-square stable
% if and only if the spectral radius of that stability matrix 
% is less than 1, rho(S) < 1, cf. Lemma 3.
%
% Input:  - F is the drift matrix,
%         - G is an array of diffusion matrices,
%         - h is the stepsize,
%         - theta is the method parameter,
%           which controls the discretization of the drift part, and
%         - com is a boolean variable, stating if the system has
%           commutative noise terms)
% Output: - sMar is the stability coeffient of Maruyama-type approx.
%         - sMil is the stability coeffient of Milstein-type approx.
%       
% 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/
%
%%%

% First, get the number of noise sources ...
m = length(G);

% Get the dimension of the drift matrix ...
dimF = size(F);
d = dimF(1);

% Start constructing the stability matrix S,
% compute matrices A and B_r, cf. Equations (15) and (16)
Inv= inv( eye(d) - h*theta*F );
A = Inv*( eye(d) + h*(1-theta)*F );
S = kron(A,A);

for r = 1:m
  B{r} = sqrt(h) * Inv * G{r};
  S = S + kron(B{r},B{r});
end

% S is now the stability matrix of Theta-Maruyama approximation,
% cf. Theorem 1.
SMar = S;

% Compute additional C matrices for Theta-Milstein approximation, 
% cf. Equation (23), compute all C's and mix them up for (24) ...
Cmix = zeros(d);
for r1 = 1:m
  for r2 = 1:m
    C{r1,r2} = 0.5 * h * Inv * G{r1}*G{r2};
    if  (r1 ~= r2)
      Cmix = Cmix + C{r1,r2};
    end
  end
end

% Compute the stability matrix S for Theta-Milstein approximation 
% with commutative noise, cf. Theorem 3.
% add the Kronecker product of the equal terms ...
for r = 1:m
    S = S + 2*kron(C{r,r},C{r,r});
end

% add the Kronecker product of the mixed terms ...
S = S + kron(Cmix,Cmix);

if (com == true)
  % done, cf Theroem 3.
  SMil = S;
else
  % add an additional term ...
  % cf. Theorem 4.
  Cmix = zeros(d);
  for r1 = 1:m
    for r2 = 1:m
      if  (r1 < r2)
        Cmix = Cmix + C{r1,r2};
      elseif (r1 > r2)
        Cmix = Cmix - C{r1,r2};
      end
    end
  end
  % compute K = sum_{k=1}^p (1/k^2)
  % where p = ceil(1/h)
  p = ceil(1/h);
  K = 0;
  for k=1:p
    K = K + 1/(k^2);
  end
  % add additional term to the stability matrix, cf. Theorem 4. 
  SMil = S + (6/pi^2) *K* kron(Cmix,Cmix);
end

% compute spectral radii, cf. Lemma 3,
% (maximum of the eigenvalues in modules)
sMar = max( abs( eig(SMar)));
sMil = max( abs( eig(SMil)));

