clc; close all; 
format long
% Laboratorio 11

%
% definisco la matrice
% A = [3,0,0;
%     1,2,0;
%     0,1,4];
% disp(A)   

% caso 1
A = load('matrice.txt');
b = [4;11;29;20];

% caso 2
% G = numgrid('S',12);
% A = delsq(G);

% caso 3
% n = 50;
% A = diag(4*ones(n,1)) + diag(-2*ones(n-1,1),1) + diag(ones(n-1,1),-1);

% dimensione del sistema
ndim = size(A,1);

% fisso il termine noto, sapendo che la soluzione esatta e' [1,1,1]
%x_exact = ones(ndim,1);
%b = A*x_exact;


%
% calcolo con il solutore lineare di Matlab
x = A\b;
disp('Soluzione Matlab')
disp(x)


%
% calcolo con il metodo di Jacobi
x0=zeros(ndim,1);
toll=1e-6;
imax=500;

[x_j,scarto_j,iter_j] = jacobi(A,b,x0,toll,imax);
disp('Soluzione Jacobi')
disp(x_j)

disp('Iterazioni')
disp(iter_j)

disp('residuo')
disp(scarto_j(iter_j))

[x_s,scarto_s,iter_s] = seidel(A,b,x0,toll,imax);
disp('Soluzione Gauss-Seidel')
disp(x_s)
disp('Iterazioni')
disp(iter_s)
disp('residuo')
disp(scarto_s(iter_s))


% grafico convergenza
figure(1)
semilogy(1:iter_j,scarto_j,'b'); hold on
semilogy(1:iter_s,scarto_s,'r')
title('Convergenza metodi iterativi')
legend('Jacobi','Gauss-Seidel')