%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        Citlivost vl. cisel a vliv FPA na vypocty
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Priklad 1: citlivost vl. cisel a vypocet pres charakteristicky polynom

A = diag(20:-1:1);              % konstrukce diagonalni matice
p = charpoly(A);                % spocti char. polynom v presnosti double
q = single(charpoly(A));        % a single

semilogy(abs(p),'o')            % vykreslime, porovname
hold on
semilogy(abs(q),'*')
pause

norm(p-q)                       % norma chyby je velka! ->
pause                           % pozor na pocet platnych cifer 

figure(2)
subplot(2,1,1)
semilogy(abs(p-q),'*')          % abs. chyba je velka pro velke koeficienty
subplot(2,1,2)
semilogy(abs(p-q)./abs(p),'*')  % rel. chyba je radu presnosti pro single
pause

% vypocet vl. cisel pres koreny char. polynomu v double

r1=roots(p);

figure(3)                       
subplot(2,2,1)                  % vykreslime a porovname s presnymi vl.c.
plot(r1,'o')
hold on
plot(diag(A),'*')

pert1=abs((diag(A)-r1));        % vykreslime abs. chybu -> ruzna, radove mala
subplot(2,2,2)
semilogy(pert1,'o')
pause

% vypocet vl. cisel pres koreny char. polynomu v single

r2=roots(q);
                     
subplot(2,2,3)                  % vykreslime a porovname s presnymi vl.c.
plot(r2,'o')
hold on
plot(diag(A),'*')               

pert2=abs((diag(A)-r2));        % vykreslime chybu -> obrovska!
subplot(2,2,4)
semilogy(pert2,'o')
pause

% srovnani s funkci eig

lambda=eig(A);                  % presna aproximace
pause

%% Priklad 2: citlivost vl. cisel a vypocet pres funkci eig

A = gallery(3);                 % spatne podminena, 3x3 

figure(4)                       % vykresleni matice
surf(A)
pause

cond(A)                         % velke -> bude citliva k perturbacim

[X,lambda] = eig(A);            % spocti vl. cisla
lambda = diag(lambda);

cond(X)                         % podminenost vl. vekt. je velka 1e-3
                                % -> lze cekat citlivost na perturbace

% perturbujeme matici a sledujeme vliv na vl. cisla

format long
lambda

delta = 1.e-6;                  % velikost perturbace
lambda1 = eig(A + delta*randn(3,3))

figure(5)  

plot(lambda,'o')                % porovname vl. cisla A a perturbovane A
hold on
plot(lambda1,'*')

pert = lambda1 - lambda         % zmena vl. c. zavisi na podminenosti X  
norm(pert)                      % o 2 rady vetsi nez velikost perturbace
pause

% zvetsujte perturbaci, sledujte vl. cisla


%% Priklad 3: citlivost vl. cisel a vypocet pres funkci eig

A = gallery(5);                 % spatne podminena, 5x5 
                                % presna vl. c. = {0,0,0,0,0}!

figure(6)                       % vykresleni matice
surf(A)
pause

cond(A)                         % velke -> bude citliva k perturbacim

% overime, ze A ma vl. c. 0 s geom. nasobnosti 1 (ma jeden Jord. blok)

A^5                             % nulova matice


% vliv FPA na vypocet vl.c. pomoci funkce eig

[X,lambda] = eig(A);            % spocti vl. cisla
lambda = diag(lambda);

cond(X)                         % podminenost vl. vekt. je velka 1e+18
                                % -> velka citlivost na perturbace

figure(7)                       % vykreslime a porovname s presnymi vl.c.

plot(real(lambda),imag(lambda),'r*',0,0,'bo')
axis square
pause

% pridame navic perturbaci a sledujeme zmeny vl.c.

eps                             % strojova presnost double

lambda1 = eig(A + eps*randn(5,5).*A);  % simulujeme single precision

hold on
plot(real(lambda1),imag(lambda1),'m*') % sleduj posun vl. cisel od nuly

% zvetsujte perturbaci, sledujte vl. cisla