% risolvere un sistema lineare mal condzionato
A = hilb(13);
% A e' simmetrica e definita positiva
% partiamo dalla soluzione precisa e ci calcoliamo i termini
% noti b
x = ones(13,1);
b = A * x;
cond(A)
cond(A,1)
cond(A,'fro')
[L, U, P] = lu(A);
function y = avanti(L,b)
y = b(1);
for k = 2 : length(b)
y = [y; b(k) - L(k, 1: k-1) * y];
end
end
function x = indietro(U,y)
n = length(y);
x = y(n)/U(n,n);
for k = n-1 : -1 : 1
x = [(y(k)-U(k,k+1 : n)*x)/U(k,k); x];
end
end
y = avanti(L, P*b);
x_lu = indietro(U, y)
eps_x_lu = norm(x - x_lu)/norm(x)
norm(x - x_lu)
norm(x)
x_div = A\b
eps_x_div = norm(x - x_div)/norm(x)
help itermeth
[x_gs, iter] = itermeth(A,b,zeros(13,1),100000,1.e-16,'G')
[x_gs, iter] = itermeth(A,b,zeros(13,1),200000,1.e-16,'G')
eps_x_gs = norm(x - x_gs)/norm(x)
quit
# Octave 3.6.2, Wed Nov 21 16:57:20 2012 CET <achilles@achilles-TOSH>
