%es 4 29.06.20

close all
clear all
warning off
clc

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

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

f1=figure;
RR=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)
    subplot(2,1,1)
    plot(R(i))
    %text(A(i,i),0,strcat('R_',num2str(i)'))
    hold on
    RR=union(RR,R(i));
end
axis equal
title('Dominio di appartenenza','Theorem of Ghershgorin')
subplot(2,1,2)
plot(RR)
axis equal

fprintf('Avrò un autovalore compreso tra -17 e -8 \n --> non è definita positiva\n')

%porto in forma tridiagonale e cerco zeri attraverso sturm (per quelli
%positivi) per quello negativo uso il Metodo delle potenze e per max
%positivo

[T] = Hess(A);
%come approssimazione prendo p1=-12, centro della mia circonferenza e p2=20, estremo massimo;
Nmax=100;
toll=sqrt(eps);
p1=-12;
p2=20;
[l1]=MetPot(T,p1,Nmax,toll);
[l2]=MetPot(T,p2,Nmax,toll);

%ora procedo a calcolare gli altri 3
x=linspace(0,l2);

f=@(x) sturm(T,x);

figure
fplot(f,[0,l2])

I3=linspace(2.5,3);
I4=linspace(4,4.5);
I5=linspace(6,7);

[l3]=secanti(f,I3,Nmax,toll);
[l4]=secanti(f,I4,Nmax,toll);
[l5]=secanti(f,I5,Nmax,toll);
L=[l3 l4 l5];
hold on
plot(L,f(L),'*g')
title('Ricerca autovalori interni')

L=[l1; L'; l2]
figure(f1)
subplot(2,1,2)
hold on
plot(complex(L),'*r')
legend('cerchi per righe','autovalori')

function [T] = Hess(A)
n=length(A);
Q=eye(n);

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

    if abs(v(1))<sqrt(eps)
        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);
A=Qi*A*Qi;
end
T=A;
end

function [l]=MetPot(T,p,Nmax,toll)
n=length(T);
H=T-p*eye(n);
y=ones(n,1);
l=1;
k=1;
e=1;
while k<Nmax && e>toll
    l0=l;
    w=H\y;
    [~,m]=max(abs(w));
    l=y(m)/w(m);
    y=w/norm(w,inf);
    k=k+1;
    e=abs(l-l0);

end
l=l+p;
end

function [pp]=sturm(T,xx)
D=diag(T);
C=diag(T,1);
pp=[];
for x=xx
p=x-x;
n=length(D);
p(1)=1;
p(2)=x-D(1);
for i=2:n
    p(i+1)=(x-D(i)).*p(i) -C(i-1)^2*p(i-1);
end
pp=[pp p(n+1)];
end

end

function [x0]=secanti(f,I,Nmax,toll)
x0=I(1);
x=I(end);
k=1;
e=abs(f(x0));

while k<Nmax && e>toll
    r=(f(x0)-f(x))/(x0-x);
    x=x0;
    x0=x0-f(x0)/r;
    k=k+1;
    e=abs(f(x0));
end
end