%metodi iterativi per sistemi lineari con
%interpretazione grafica per il pb di minimizzazione associato (se dim(A)=2)
%Ax=b  <--> min 1/2 <Ax,x> - <b,x>

clear
clc

N=2; %dimensione matrice
toll=10e-15; %precisione approssimazione su || x^(k+1)-x^k ||
it_max=1000;  %numero massimo di iterazioni

%matrice A (diag dominante)
A=rand(N,N)*6; 
for i=1:N
  A(i,i)=sum(A(i,:))+1; %correzione elementi sulla diagonale
end

%vettore b
b=rand(1,N);

%punto iniziale x^0
x0=rand(1,N)*4;

%::::::::::::: CORE ::::::::::::::::::::::
%P e' una matrice che contiene il risultato di ogni iterazione: La prima
%riga e' x^0, la k-ma riga e' x^(k-1).
%it e' il numero di righe di P (numero di iterazioni effettuate +1)

[P,it]=GaussSeidel(A,b,x0,toll,it_max);
%[P,it]=Jacobi(A,b,x0,toll,it_max);
%:::::::::::::::::::::::::::::::::::::::::

soluz=P(it,:)'  %soluzione ( ultima delle x^k) )

norm(A*soluz - b',1)  %residuo finale ( norma 1 del vettore A*soluz-b )

%-------------------------------------
%plot (solo se N=2)
if N==2
 [x,y]=meshgrid(-5:0.1:5,-4:0.1:4);
 %creazione della matrice z contenente i valori della funzione 1/2 <Ax,x> - <b,x>, x in R^2
 for i=1:size(x,1)
    for j=1:size(y,2)
       v=[x(i,j);y(i,j)];
       z(i,j)=( 0.5*dot(A*v,v) ) - dot(b,v);
    end
 end
 contour(x,y,z,60)  %disegna le curve di livello di 1/2 <Ax,x> - <b,x>
 hold on
 plot(P(:,1),P(:,2),'ks-','markersize',2) %disegna gli x^k
 figure %apre una nuova finestra grafica
 mesh(x,y,z) %disegna il grafico di 1/2 <Ax,x> - <b,x>
end

