%esercizio 2 21.01.20
clear all
close all
clc

alphav=linspace(5,10);
[A,a]=Aa(alphav);
b=[26;14;7];

Nmax=100;
toll=sqrt(eps);

[xj] = Jacobi(A,b,Nmax,toll);
[xg] = GaussSeidel(A,b,Nmax,toll);


function [A,a]=Aa(alphav)

%norm(A,2) = sqrt(RhoS(A'*A)) dove rhos(A'*A) = raggio spettr
n=length(alphav);
Nmax=100;
toll=sqrt(eps);
for av=alphav
    Ai=[7 2 1;2 av -2;1 -2 8];
    B=Ai'*Ai;
    [l]=MetPot(B,Nmax,toll);

    if av==alphav(1)
        l0=l;
        a=av;
        A=Ai;
    else 
        if l<l0
        l0=l;
        a=av;
        A=Ai;
        end
    end 
end

end

function [l]=MetPot(B,Nmax,toll)
y=ones(length(B),1);
k=1;
l=1;
e=1;

while k<Nmax && e>toll
    l0=l;
    w=B*y;
    [~,m]=max(w);
    l=w(m)/y(m);
    y=w/norm(w,inf);
    e=abs(l-l0);
    k=k+1;
end
   
end

function [xj] = Jacobi(A,b,Nmax,toll)

J=-diag(1./diag(A))*(A-diag(diag(A)));
[rhoj]=MetPot(J,Nmax,toll);

if rhoj>=1
    disp('Metodo di Jacobi non convergente')
    return
end

k=1;
hj=diag(1./diag(A))*b;
x0=zeros(size(b));
xj=hj;
e=1;

while k<Nmax && e>toll
    xj= J*xj +hj;
    e=norm(xj-x0)/norm(xj);
    k=k+1;
end
   
end

function [xg] = GaussSeidel(A,b,Nmax,toll)

G=-tril(A)\(A-tril(A));
[rhog]=MetPot(G,Nmax,toll);

if rhog>=1
    disp('Metodo di Gauss-Seidel non convergente')
    return
end

k=1;
hg=tril(A)\b;
x0=zeros(size(b));
xg=hg;
e=1;

while k<Nmax && e>toll
    xg= G*xg +hg;
    e=norm(xg-x0)/norm(xg);
    k=k+1;
end

end