%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         Řešení lin. soustavy citlivé na chyby v datech
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Vytvoříme sekvenci Hilbertových matic A rostoucí dimenze. Budeme řešit
% soustavy Ax=b pro x=(1, ...,1) a sledovat chyby řešení.

help hilb                      % Hankelovská, velmi špatně podmíněná

kappa=[];                      % podmíněnost
res = [];                      % norma residua pro aproximace řešení
bcerror=[];                    % odhad relativní zpětné chyby

error=[];                      % relativní norma skutečné chyby řešení
est  =[];                      % odhad relativní chyby řešení z teorie citlivosti

for n=1:20 
    A = hilb(n);               % vytvoř matici
    kappa = [kappa, cond(A)];  
    
    xex=ones(n,1);             % přesné řešení ze samých jedniček
    b=A*xex;
    x=A\b;                     % spočtená aproximace řešení

    res=[res,norm(b-A*x)];     
    bcerror=[bcerror, norm(b-A*x)/(norm(A)*norm(x)+norm(b)) ];

    error=[error,(norm(xex-x)/norm(xex))];
    est=[est, kappa(n)*eps];
end

subplot(1,4,1)                
semilogy(kappa,'*')             % podmíněnost rychle roste
title('Podminenost A')
subplot(1,4,2) 
semilogy(res,'-*')              % norma residua je malá
title('||b-A*x||')
subplot(1,4,3)
semilogy(bcerror,'-*')          % relativní zpětná chyba je malá - vyřešili 
title('Odhad zpetne chyby')                    % jsme úlohu blízkou původní

subplot(1,4,4)
semilogy(error,'-*')            % přesto je skutečná chyba řešení velká -> 
hold on                         % citlivost úlohy na perturbace dat
semilogy(est,'g-o')             % odpovídá odhadu z teorie citlivosti 
title('Skutecna relativni chyba a jeji odhad') 

% zopakujeme experiment pro single precision ... A = single(hilb(n));
% opravte eps
% -> skutečná chyba roste ještě rychleji