% Soluzione di un sistema lineare mal condizionato
% con metodi vari e confronto degli errori,
% si veda l'esercizio 4 del 09/12/2008 
% (http://www.dm.unibo.it/~achilles/calc/calc08-12-09.pdf)
% e Quarteroni-Saleri, esempio 5.9, pag. 147.
% I programmi da utilizzare sono:
% http://www.dm.unibo.it/~achilles/calc/Programmi_CS4a/itermeth.m
% http://www.dm.unibo.it/~achilles/calc/Programmi_CS4a/sol.m


A = hilb(13)
% numero di condzionamento rispetto alla norma 2
cond(A)
b=A*ones(13,1)
[L, U, P] = lu(A)
bperm = P*b;
% sostituzione in avanti
y(1) = bperm(1);
for k = 2:13 y(k)=bperm(k)-L(k,1:k-1)*y(1:k-1)';end;
% sostituzione all'indietro
xlu(13)=y(13)/U(13,13);
for k=12:-1:1 xlu(k)=(y(k)-U(k,k+1:13)*xlu(k+1:13)')/U(k,k);end;
xlu'
errlu = norm(ones(13,1) - xlu')

% Il calcolo precedente puo' esssere eseguito piu' comodamente 
% con lo script sol.m.
xsol = sol(A,b)
errsol = norm(ones(13,1) - xsol)

% Soluzione del sistema lineare usando \.
x = A\b
err = norm (ones(13,1) - x)

% Soluzione del sistema lineare usando il metodo iterativo di Gauss-Seidel,
% implementato nello script itermeth.m.
[xgs, iter]  = itermeth(A,b,zeros(13,1),20000,1.e-10,'G')
errgs = norm (ones(13,1) - xgs)