%
%  risoluzione di un sistema lineare fortemente mal-condizionato
%
stde = 1*1e-6;  % deviazione standard del rumore additivo aggiunto
n=64;
[U, R] = qr(rand(n)); 
[V, R] = qr(rand(n)); 
if 0
  S = diag(exp(1e1*[0:-1:-n+1]));
else
  S = diag(exp(0.6*[0:-1:-n+1]));
end
figure(1), semilogy(log10(diag(S)));
A = U*S*V';
K2_A = cond(A);
disp(['numero di condizionamento di A = ' num2str(K2_A,'%g')]);
x_vero = ones(n,1);  % soluzione vera (imposta da noi)
b = A*x_vero + stde*randn(n,1);
x_LU = A \ b;
figure(3); clf, plot(x_vero,'b-'); hold on, pause, plot(x_LU,'r'); pause,
%
n2res_hist = zeros(n,1);
for k=1:n
  disp(['k=' num2str(k)])
  x_TSVD = V(:,1:k)*diag(1./diag(S(1:k,1:k)))*U(:,1:k)' * b;
  n2xTSVD = norm(x_TSVD)
  A_TSVD = U(:,1:k)*S(1:k,1:k)*V(:,1:k)';
  n2res = norm(A_TSVD*x_TSVD - b)
  n2res_hist(k) = n2res;
  figure(3); plot(x_TSVD,'m'); %pause,
  figure(4); clf, hold on, plot(x_vero,'b-'); plot(x_TSVD,'m'); pause,
end

figure(10); clf, hold on,
            semilogy([1 n],log10([stde stde]),'r--');
            semilogy(log10(n2res_hist),'b.');
            axis([1 n log10(1e-18) log10(1e1)]);

            