%esempio di successione di sturm
clear all
close all

clc

A=[1/6 0 1/3 0 1/3; 0 1/12 1/3 1/6 0; 1/3 1/3 0 0 0; 0 1/6 0 1/3 0; 1/3 0 0 0 1/3];
%Costruisco una matrice tridiagonale essendo A simmetrica;
[H] = Hessemberg(A);

%gli autovalori sono i valori che mi annullano il termine n-esimo della
%successione di sturm. Posso dunque dire che tutte i lambda che mi annullano
%la 'funzione' sono tra -1 ed 1;
x=linspace (-1,1);
f=@(x) ST(H,x);

figure

fplot(f,'g',[-1,1])
yline(0)

legend('termine n-esimo successione di sturm','Location','southeast')
title ('Polinomio caratteristico')

figure
fplot(f,'g',[0.2,0.3])
xx=linspace(0.2,0.3);
Nmax=100;
toll=sqrt(eps);

[lambda0]=MetodoSecanti(f,xx,Nmax,toll);
hold on
plot(lambda0,f(lambda0),'*r')
yline(0,'--k')
title('Zoom del mio polinomio caratteristico nell''intervallo d''interesse', 'FontSize',10)
legend('Zoom pol caratt','( \lambda_0 , p(\lambda_0) )','Location','southeast')
xlabel('0.2 \leq \lambda \leq 0.3')
ylabel('p(\lambda)')

function [H] = Hessemberg(A)
%costruisco tridiag

[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);
    A=Qi*A*Qi;
    Q=Qi*Q;
end
H=A;
end

function [p] = ST(H,x)
%costruisce polinomio caratteristico
n=length(H);
d=diag(H);
c=diag(H,-1);
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
p=p(n+1);
end

function [x0]=MetodoSecanti(f,xx,Nmax,toll)
%cerco autovalore tra 0.2 e 0.3, ce n'è uno solo (cerco uno 0)

x0=xx(1);
x1=xx(end);
r=(f(x0)-f(x1))/(x0-x1);
x0=x0-f(x0)/r;
k=1;
e(k)=abs(f(x0));

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

end
end