%Gauss-Seidel
clear all
clc

A=[4 1;1 1];
b=sum(A,2);

[x,e,k]=GaussSeidel(A,b,100,eps);
iter=linspace(1,k,k);
plot(e,'g');
hold on
plot(iter,e(iter),'*r');
legend('Andamento asintotico del mio errore','Corrispondenza iter-errore')
xlabel('numero di iterazioni')

function [x,e,k] = GaussSeidel(A,b,Nmax,tol)

%metodo di iterazione di Gauss-Seidel
%A=L+U, dove L=parte triangolare di A
%B=-inv(L)*U;
%Gauss Seidel converge se il raggio spettrale di B è <1

n=length(A);
x0=zeros(n,1);
L=tril(A);
B=-L\(A-L);
h=L\b;
x=B*x0+h;
e(1)=norm(x-x0)/norm(x);
RHO=max(abs(eig(B)));

if RHO<1
    disp('Il metodo converge!');
    k=1;
    while k<Nmax && e(k)>eps
        x0=x;
        x=B*x + h;
        k=k+1;
        e(k)=norm(x-x0)/norm(x);
    end
else
    disp('Il metodo di Gauss-Seidel non converge!');
end
end