%====================================================================
% Script per determinare la soluzione di un sistema lineare
% con i due metodi di rilassamento di Jacobi e di Gauss Seidel
% Necessita delle function jacobi e gausseidel
%====================================================================

clear all

% Definizione dati 
n = 50;
A = diag(4*ones(n,1)) + diag(-2*ones(n-1,1),1) + diag(ones(n-1,1),-1);
sol = ones(n,1);
b = A*sol;
toll = 1e-5;
x0 = zeros(n,1);
kmax = 50;

% Controlli generali
[mA,nA] = size(A);
if (mA~=nA)
    error('La matrice A non e'' quadrata')
end

if length(b) ~= nA
    error('Le dimensioni di A e b non sono compatibili')
end

% Risoluzione

% Metodo di Jacobi
[x,k,flag] = jacobi (A,b,x0,toll,kmax);

% Metodo di Gauss Seidel
[xg,kg,flagg] = gaussseidel (A,b,x0,toll,kmax);

% Controlli
if (flag ~= 0) || (flagg ~= 0)  
    error('Non si possono applicare i metodi di rilassamento!') 
end

disp('Numero di iterate per Jacobi')
disp(k)
disp('Numero di iterate per Gauss Seidel')
disp(kg)
if k == kmax 
    disp('Massimo numero di iterate raggiunto per Jacobi')
elseif kg == kmax
    disp('Massimo numero di iterate raggiunto per Gauss Seidel')
end

% Uscita risultati
disp('Dimensione del sistema')
disp(mA)
disp('Tolleranza')
disp(toll)

% Grafico

for j = 1:k+1
    err(j) = norm(x(:,j)-sol);
end
for j = 1:kg+1
    errg(j) = norm(xg(:,j)-sol);
end
xk = 0:k; xkg = 0:kg;
semilogy(xk,err,'r',xkg,errg,'b');
title('Metodi di Jacobi e di Gauss Seidel')
xlabel('iterate')
ylabel('norma dell''errore')
legend('Jacobi','Gauss Seidel','Location','best')

