% test_lyap
%
% Berechnung des Cholesky-Faktors L der L"osung Q=L*L' der 
% Lyapunov-Gleichung

clear all;
figure(1); clf;
q = 200; 
tol = 1e-10;

% Systemmatrizen
[A, B, C] = baum(q);
ev = eig(A);
subplot(2,2,1)
plot(real(ev), imag(ev), 'b*');
title('Eigenvalues von A')

% Frequenzgang G(iw)
disp('Berechnung des Frequenzgangs');
wmin  = 0;
wmax  = 2;
numw  = 100;
wh    = (wmax-wmin)/numw;
w     = logspace(wmin,wmax,numw);
i     = sqrt(-1);
for j = 1:numw
   iwj = i*w(j);
   G   = sparse(C)*(sparse(iwj*speye(2*q)-A)\sparse(B));
   nrmG(j) = norm(full(G));
end 
subplot(2,2,2);
semilogx(w, nrmG(:), 'LineWidth', 1.5);
xlabel('Frequenz w');
ylabel('Magnitude');
title('Frequenzgang G(iw)');

% L"osung der Lyapunov-Gleichung A`*Q+Q*A=-C`*C
disp('  ');
disp('Loesung der Lyapunov-Gleichung A`*Q+Q*A=-C`*C');
[L, convC] = sfm_lyap(A, C, tol);
nresC = norm(A'*L*L'+L*L'*A+C'*C,'fro')/norm(C,'fro')^2
subplot(2,2,3);
semilogy(1:length(convC), convC, 'b*-');
xlim([1 length(convC)]);
xlabel('Iteration');
ylabel('||Ai+I||_F');
title('Konvergenz der Signumfunktionsmethode');

% L"osung der Lyapunov-Gleichung A*P+P*A`=-B*B`
disp('  ');
disp('Loesung der Lyapunov-Gleichung A*P+P*A`=-B*B`');
% ACHTUNG!!! 
%Auskommentieren Sie diese Zeile und setzen die richtigen Matrizen ein
%[R, convB] = sfm_lyap(??, ??, tol);
nresB = norm(A*R*R'+R*R'*A'+B*B','fro')/norm(B,'fro')^2
subplot(2,2,4);
semilogy(1:length(convB), convB, 'b*-');
xlim([1 length(convB)]);
xlabel('Iteration');
ylabel('||Ai+I||_F');
title('Konvergenz der Signumfunktionsmethode');