test_lyap2.m

function pass = test_lyap2()

% Test LYAP for some precalculated matrices.

tol = 1e-10;
n = 3;

A = gallery('frank', n);
Q = A;
I = eye(n);

% A = rand(n) + 1i*rand(n);
% Q = rand(n);

%%
U = lyap(A, Q);
V = lyap(A, A', Q);
% W = lyap3(A, A', Q);
W = [ ...
  -0.500000000000000   0.071428571428571  -0.214285714285714
  -0.071428571428572  -0.500000000000000   0.142857142857143
   0.214285714285715  -0.142857142857143  -0.500000000000000];
err(1) = norm(U - V) + norm(U - W);

%%
A = A + 1i*I;
U = lyap(A, A.', Q);
% V = lyap3(A, A.', Q);
V = [ ...
  -0.337804878048780  -0.163414634146341  -0.075609756097561
  -0.178048780487805  -0.200000000000000  -0.143902439024390
   0.104878048780488  -0.173170731707317  -0.207317073170732] + ...
    1i*[ ...
   0.159756097560976  -0.070731707317073   0.019512195121951
  -0.002439024390244   0.100000000000000   0.004878048780488
  -0.156097560975610   0.141463414634146   0.134146341463415];
err(2,1) = norm(U - V);

%%
A = real(A);
B = rot90(A);
U = lyap(A, B, Q);
% V = lyap3(A, rot90(A), Q);
V = [ ...
   0.469283276450512   0.315699658703072  -1.187713310580205
   0.085324232081911  -0.488054607508532  -0.034129692832765
  -2.116040955631399  -0.696245733788396   1.846416382252561];

err(3,1) = norm(U - V);

%%

[ZA, TA] = schur(A, 'complex');
[ZB, TB] = schur(B, 'complex');
U = ZA * lyap(TA, TB, ZA'*Q*ZB) * ZB';
% V = lyap3(A, B, Q);
V = [ ...
   0.469283276450512   0.315699658703072  -1.187713310580205
   0.085324232081911  -0.488054607508532  -0.034129692832765
  -2.116040955631399  -0.696245733788396   1.846416382252561];

err(4,1) = norm(U - V);

%%

pass = err < tol;

end

function varargout = lyap3(varargin)
currDir = cd;
cd(fileparts(which('sylvslv')))
[varargout{1:nargout}] = lyap(varargin{:});
cd(currDir)
end