%es 1 19/09/2022
clear all
close all
clc

A=[0 0 1 0 0; 0 10 2 1 1; 1 2 6 0 0;0 1 0 4 0;0 1 0 0 2];
eig(A)
f1=figure;
hold on
f2=figure;
hold on

RR=polyshape([0],[0]);
CC=polyshape([0],[0]);

t=linspace(0,2*pi);
t=t(1:end-1);

for i=1:length(A)
    r(i)=sum(abs(A(i,[1:i-1,i+1:end])));
    R(i)=polyshape(r(i)*cos(t)+A(i,i),r(i)*sin(t));
    figure(f1)
    hold on
    plot(R(i))

    c(i)=sum(abs(A([1:i-1,i+1:end],i)));
    C(i)=polyshape(c(i)*cos(t)+A(i,i),c(i)*sin(t));
    figure(f2)
    hold on
    plot(C(i))
    RR=union(RR,R(i));
    CC=union(CC,R(i));
end
%coincidono perchè è una matrice simmetrica;

figure(f1)
axis equal
title ('Cerchi per righe')

figure(f2)
axis equal
title ('Cerchi per colonne')

f3=figure;
U=intersect(RR,CC);
plot(U)
axis equal

title('Dominio di appartenenza','Theorem of Ghershgorin','FontAngle','italic')
%si vede che ho al massimo un solo autovalore di A negativo

%Applico fattorizzazione per portare matrice in forma tridiagonale
[T] = QRH(A);
[r] = STURM(T);

fprintf('ho %d autovalori positivi\n',r)

p=4;
Nmax=100;
toll=sqrt(eps);
[l] = MetPot(T,p,Nmax,toll)

function [T] = QRH(A)
n=length(A);

Q=eye(n);

for i=1:n-1
    v=A(i+1:n,i);
    u=zeros(n,1);

    if v(1)==0
        u(i+1:n)=v + norm(v)*eye(n-i,1);
    else
        u(i+1:n)=v + sign(v(1))*norm(v)*eye(n-i,1);
    end

    Qi=eye(n) -2*(u*u')/(u'*u);
    Q=Qi*Q;
    A=Qi*A*Qi;
end
T=A;
end

function [r] = STURM(T)
D=diag(T);
C=diag(T,1);
n=length(D);

p(1)=1;
p(2)=-D(1);
P=[p(1)];

for i=2:n
    p(i+1)=-D(i)*p(i) -C(i-1)^2*p(i-1);
    if (p(i)~=0)
        P=[P p(i)];
    end
end

if (p(n+1)~=0)
        P=[P p(n+1)];
end
r=0;
for i=1:length(P)-1
    if P(i+1)*P(i)<0
        r=r+1;
    end
end
end

function [l] = MetPot(T,p,Nmax,toll)
n=length(T);
y=ones(n,1);
l=1;
k=1;
F=T-p*eye(n);
[L,U]=luM(F);
e=1;

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

l=1/l +p;

end

function [L,U] =luM(F)

n=length(F);
L=eye(n);
U(1,:)=F(1,:);

for k=1:n
    U(k,k:n)=F(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);
    L(k+1:n,k)=(F(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,k))/U(k,k);
end

end

function [w] = solve(L,U,y)

n=length(L);
a=eye(n,1);
a(1)=y(1)/L(1,1);

for j=2:n
    a(j)=(y(j) - L(j,1:j-1)*a(1:j-1))/L(j,j);
end

w=eye(n,1);
w(n)=a(n)/U(n,n);

for i=n-1:-1:1
    w(i)=(a(i) - U(i,i+1:n)*w(i+1:n))/U(i,i);
end

end