%es 2 23.07.21
close all
clear all
clc

alphav=linspace(-0.5,3);
Nmax=100;
toll=sqrt(eps);
[A,b,J,G,a] = Aa(alphav,Nmax,toll);
[xj,ej,rj] = JACOBI(A,J,b,Nmax,toll);
[xg,eg,rg] = GaussSeidel(A,G,b,Nmax,toll);

f1=figure;

subplot(2,1,1)
hold on
plot (rj,'g')
plot(rg,'b')
title('CONVERGENZE')
subplot(2,1,2)
hold on
plot (log10(rj),'g')
plot(log10(rg),'b')

f2=figure;

subplot(2,1,1)
hold on
plot (ej,'g')
plot(eg,'b')
title('CONVERGENZE')
subplot(2,1,2)
hold on
plot (log10(ej),'g')
plot(log10(eg),'b')


function [A,b,J,G,a] = Aa(alphav,Nmax,toll)

for av=alphav
    Ai=[8 1 6 0 5; 1 5 2 1 av; 6 2 10 5 1; 0 1 1 6 10; av 0 1 2 -4];
    D=diag(diag(Ai));
    J=-D\(Ai-D);
    G=-tril(Ai)\(Ai-tril(Ai));
    [esitoj]=MetPotenze(J,Nmax,toll);
    [esitog]=MetPotenze(G,Nmax,toll);

    if esitoj==1 && esitog==1
        a=av;
        A=Ai;
        b=[a 1 1+a 0 1]';
        return
    end
end

end

function [esito] = MetPotenze(T,Nmax,toll)
k=1;
l=1;
y=ones(length(T),1);
k=1;
e=1;

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

function [xj,ej,rj] = JACOBI(A,J,b,Nmax,toll)
hj=diag(diag(A))\b;
x0=zeros(size(b));
xj=hj;
k=1;
rj(k)=norm(A*xj-b);
ej(k)=norm(xj-x0);

while k<Nmax && rj(k)>sqrt(toll)
    xj=J*xj +hj;
    k=k+1;
    rj(k)=norm(A*xj-b);
    ej(k)=norm(xj-x0);
end

end

function [xg,eg,rg] = GaussSeidel(A,G,b,Nmax,toll)
hg=tril(A)\b;
x0=zeros(size(b));
xg=hg;
k=1;
rg(k)=norm(A*xg-b);
eg(k)=norm(xg-x0);
while k<Nmax && rg(k)>sqrt(toll)
    xg=G*xg +hg;
    k=k+1;
    rg(k)=norm(A*xg-b);
    eg(k)=norm(xg-x0);
end

end