% -----------------------------------%
% Esercitazione del 23 gennaio 2006  %
% -----------------------------------%

% Primo esercizio

alfa=0.5:0.001:.999;
for k=1:length(alfa),
  A=[1 alfa(k) alfa(k); alfa(k) 1 alfa(k); alfa(k) alfa(k) 1];
 % matrice di Jacobi
 D=diag(diag(A));
 DI=inv(D);
 J=-DI*(A-D);
 
 % matrice di GS
L=tril(A); 
U=triu(A);
L=L-D;
U=U-D;
DI=inv(D+L);
GS=-DI*U;

rhoJ(k)=max(abs(eig(J)));
rhoGS(k)=max(abs(eig(GS)));


end;
plot(alfa,rhoJ,'.-r',alfa,rhoGS,'.-b');
legend('Raggio spett. di Jacobi','Raggio spett. di Gauss-Seidel');

alpha=2/3;
b=[1 -1 3];
A=[1 alpha alpha;  alpha 1 alpha; alpha alpha 1];
[sol,iter]=GaussSeidel(A,b,zeros(3,1),1.e-6,100)

disp('un tasto qualsiasi')
pause

% -------------- VARIANTE ---------------------------------
% Al variare di alfa facciamo vedere come varia il numero
% delle iterazioni che dipenderanno dal valore del raggio
% spettrale, ovvero ci aspettiamo 
% valori minori per piccoli valori del raggio spettrale
% della matrice d'iterazione.
%----------------------------------------------------------
alfa=0.5:0.01:.999;
for k=1:length(alfa),
  A=[1 alfa(k) alfa(k); alfa(k) 1 alfa(k); alfa(k) alfa(k) 1];
  [sol,iter(k)]=GaussSeidel(A,b,zeros(3,1),1.e-6,200);
end
bar(0.5:0.01:.999,iter);
title('Iterazioni con Gauss Seidel');


%------------ Le due funzioni ------------------

function [sol,iter]=Jacobi(A,b,x0,tol,kmax)
% ----------------------------------------------
% Funzione che implementa il metodo iterativo
% di Jacobi
%-----------------------------------------------
% Inputs
% A, b: matrice e termine noto, rispettivamente
% x0  : guess iniziale
% tol : tolleranza calcoli
% 
% Outputs
% sol : vettore soluzione
% iter: numero delle iterazioni
%-----------------------------------------------

D=diag(diag(A));
DI=inv(D);
J=-DI*(A-D);
b1=DI*b';
x1=J*x0+b1;
k=1;

while(abs(x1-x0)>tol*abs(x1) & k<=kmax)
    x0=x1;
    x1=J*x0+b1;
    k=k+1;
end

sol=x1;
iter=k-1;

%--------------------------------------------------------


function [sol,iter]=GaussSeidel(A,b,x0,tol,kmax)
% ----------------------------------------------
% Funzione che implementa il metodo iterativo
% di Gauss-Seidel
%-----------------------------------------------
% Inputs
% A, b: matrice e termine noto, rispettivamente
% x0  : guess iniziale
% tol : tolleranza calcoli
% 
% Outputs
% sol : vettore soluzione
% iter: numero delle iterazioni
%-----------------------------------------------

L=tril(A); U=triu(A);
D=diag(diag(A));
L=L-D;
U=U-D;

DI=inv(D+L);
GS=-DI*U;
b1=DI*b';
x1=GS*x0+b1;
k=1;

while(abs(x1-x0)>tol*abs(x1) & k<=kmax)
    x0=x1;
    x1=GS*x0+b1;
    k=k+1;
end

sol=x1;
iter=k-1;