%es 1 18.02.2020

close all
clear all
clc

A=[4 1 2; 0 7 2; 2 1 -5];
B=A'*A;

% Porto in forma triadiagonale (prima verifico la simmetria);

q=0;
n=length(B);

for i=2:n
    for j=1:i-1
        if B(i,j)==B(j,i)
            q=q+1;
        end
    end
end

if q~=(n^2/2 - n/2)
    disp('Matrice non simmetrica')
    return
end
    
[T]=HESS(B);

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);

n=length(T);

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

    c(i)=sum(abs(T([1:i-1,i+1:n],i)));
    C(i)=polyshape(r(i)*cos(t) + T(i,i),r(i)*sin(t));
    figure(f2)
    hold on
    plot(C(i))

    RR=union(RR,R(i));
    CC=union(CC,C(i));
end

close(f1)
close(f2)
f3=figure;
U=intersect(RR,CC);
plot(U)
title('Dominio di appartenenza dei miei autovalori','Theorem of Ghershgorin')

%dunque so che ho tre autovalori, il massimo in modulo è quello che mi
%interessa per conoscere la norma;
%
% noto anche che sono tutti positivi (vista la definita positività di una
%matrice ottenuta attraverso il prodotto matriciale tra una matrice e la
%sua trasposta.
%
%Come metodo alternativo dunque potrei anche plottare il mio termine
%n-esimo della mia espansione in serie di Schurm e cercare il suo zero in
%un intervallo scelto propriamente;
Nmax=100;
toll=sqrt(eps);

[l0] = Pot(T,Nmax,toll)
xx=linspace(30,60);
f=@(xx) SCHURM(T,xx);

f4=figure;
subplot(2,1,2)

%scelgo come intervallo I=[50,60]; 
I=linspace(50,60);
[l1] = SECANTI(f,I,Nmax,toll)

subplot(2,1,1)
hold on
fplot(f,[10,70],'g')
subplot(2,1,2)
hold on
fplot(f,[50,60],'g')
plot(l1,f(l1),'*r')
yline(0,'--k')

norm2_0=sqrt(l0);
norm2_1=sqrt(l1);

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);
    Q=Qi*Q;
    B=Qi*B*Qi;
end
T=B;
end

function [l0] = Pot(T,Nmax,toll)
n=length(T);
y=ones(n);
l=1;
k=1;
e=1;

while k<Nmax && e>toll
    l0=l;
    w=T*y;
    [M,m]=max(abs(w));
    l=w(m)/y(m);
    y=w/M;
    e=abs(l-l0);
end
l0=l;
end

function [yy] = SCHURM(T,xx)
D=diag(T);
C=diag(T,1);

n=length(D);
yy=[];
p(1)=1;
for x=xx
    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
    yy=[yy p(n+1)];
end
end

function [x0] = SECANTI(f,I,Nmax,toll)
x0=I(1);
x=I(end);
e=1;
k=0;
while k<Nmax && e>toll
    k=k+1;
    r=(f(x0)-f(x))/(x0-x);
    x=x0;
    x0=x0-f(x0)/r;
    e=abs(f(x0));
end
end