%es 2 08.06.21

close all
clear all
clc

B=[8 1 0 0 1;1 5 2 1 0;0 2 10 1 1;0 1 1 6 2;1 0 1 2 -12];

%visualizzazione di Ghershgorin, svolgo solo per righe data la simmetria

RR=polyshape([0],[0]);
t=linspace(0,2*pi);
t=t(1:end-1);

n=length(B);

for i=1:n
    r(i) = sum(abs(B(i,[1:i-1,i+1:n])));
    R(i) = polyshape(r(i)*cos(t) + B(i,i), r(i)*sin(t));
    RR=union(R(i),RR);
end

f1=figure;
plot(RR)
title('Dominio di appartenenza','Theorem of Ghershgorin','FontAngle','italic')
axis equal

% 1. Svolgo verificando il # dei cambiamenti di segno nella sequenza di sturm
% Prima porto in forma tridiagonale
[T] = Hess(B);
n1=3;
n2=6;
[q] = sturm(T,n1,n2);
fprintf('Ho %d autovalori compresi tra %d e %d\n',q,n1,n2)

% 2. Attraverso la visualizzazione di Ghershgorin noto che il mio
% autovalore è compreso tra -16 e -8, prendo come stima iniziale -12

p=-12;
Nmax=100;
toll=sqrt(eps);
[l] = PotInv(T,p,Nmax,toll);
fprintf('Il mio unico autovalore negativo è %d\n',l) %-12.3067

function [T] = Hess(B)
n=length(B);

Q=eye(n);

for i=1:n-1
    v=B(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);
    B=Qi*B*Qi;
end
T=B;
end

function [q] = sturm(T,n1,n2)

D=diag(T);
C=diag(T,1);
n=length(D);
p1(1)=1;
p1(2)=n1-D(1);
P1=[p1(1)];
for i=2:n
    p1(i+1)=(n1-D(i))*p1(i) - C(i-1)^2*p1(i-1);
    if p1(i)~=0
        P1=[P1 p1(i)];
    end
end
if p1(n+1)~=0
    P1=[P1 p1(n+1)];
end
p2(1)=1;
p2(2)=n2-D(1);
P2=[p2(1)];
for i=2:length(D)
    p2(i+1)=(n2-D(i))*p2(i) - C(i-1)^2*p2(i-1);
    if p2(i)~=0
        P2=[P2 p2(i)];
    end
end
if p2(n+1)~=0
    P2=[P2 p2(n+1)];
end
l=0;
for i=1:length(P1)-1
    if P1(i+1)*P1(i)<0
        l=l+1;
    end
end
t=0;
for i=1:length(P2)-1
    if P2(i+1)*P2(i)<0
        t=t+1;
    end
end
q=l-t;
end

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

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

l = l +p;

end

function [L,U] = LU(H)
n=length(H);

L=eye(n);
U(1,:)=H(1,:);

for k=1:n
    U(k,k:n)= H(k,k:n) - L(k,1:k-1)*U(1:k-1,k:n);
    L(k+1:n,k)= (H(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(y);
t=eye(size(y));
t(1)=y(1)/L(1,1);

for i=2:n
    t(i)=(y(i) - L(i,1:i-1)*t(1:i-1))/L(i,i); %L(i,i) =1
end

w=eye(size(y));
w(n)=t(n)/U(n,n);

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