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