%esercizio 2 22.06.22

clear all
close all
clc

A=[10 0 1 0; 0 5 2 -1; 1 2 8 2; 0 -1 2 7];
b=[11; 6;13;8];
[H] = QRH(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);
% 
% n=length(H);
% 
% for i=1:n
%     r(i) = sum(abs(H(i,[1:i-1,i+1:n])));
%     R(i) = polyshape(r(i)*cos(t)+ H(i,i), r(i)*sin(t));
%     figure(f1)
%     subplot(2,1,1)
%     hold on
%     plot(R(i))
%     text(H(i,i),0, strcat('R_',num2str(i)'))
% 
%     c(i) = sum(abs(H([1:i-1,i+1:n],i)));
%     C(i) = polyshape(c(i)*cos(t)+ H(i,i), c(i)*sin(t));
%     figure(f2)
%     subplot(2,1,1)
%     hold on
%     plot(C(i))
%     text(H(i,i),0, strcat('C_',num2str(i)'))
%     
%     RR=union(RR,R(i));
%     CC=union(CC,C(i));
% end
% 
% figure(f1)
% subplot(2,1,2)
% hold on
% plot(RR)
% 
% figure(f2)
% subplot(2,1,2)
% hold on
% plot(CC)
% 
% U=intersect(RR,CC)
% 
% f3=figure
% 
% figure(f3)
% hold on
% plot(U)
% plot(complex(eig(H)),'*r')

xx=linspace(1,12);
f=@(xx) Psturm(xx,H)
figure
fplot(f,[1,12])
I1=linspace(3,4);
I2=linspace(7,8);
I3=linspace(9,10);
I4=linspace(10,11);


Nmax=200;
toll=sqrt(eps);

[l1]=SEC(f,I1,Nmax,toll);
[l2]=SEC(f,I2,Nmax,toll);
[l3]=SEC(f,I3,Nmax,toll);
[l4]=SEC(f,I4,Nmax,toll);

l= [l1, l2, l3, l4]
hold on
plot(l,f(l),'*g')

[y1,~] = POT(H,l1,toll,Nmax);
[y2,~] = POT(H,l2,toll,Nmax);
[y3,~] = POT(H,l3,toll,Nmax);
[y4,~] = POT(H,l4,toll,Nmax);

[a,b]=eigs(H)

Y=[y1,y2,y3,y4];
a1=a(:,4)
a2=a(:,3)
a3=a(:,2)
a4=a(:,1)

rank([y1,a1])
rank([y2,a2])
rank([y3,a3])
rank([y4,a4])

[xc]=CHOL(A,b);

B=A'*A;
[~,t] = POT(H,0,toll,Nmax);
norm2=sqrt(t)

function [H] = 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
H=A;
end

function [P] = Psturm(xx,H)

D=diag(H);
C=diag(H,1);
n=length(D);
j=1;
for x=xx
p(1)=1;
p(2)=x-D(1)*ones(1,length(x));

for i=2:n
    p(i+1)= (x-D(i)*ones(1,length(x))).*p(:,i) - C(i-1)^2*p(:,i-1);
end
P(j)=p(n+1);
j=j+1;
end

end 

function [x0] = SEC(f,I,Nmax,toll)

x0=I(1);
x=I(end);
r=(f(x0)-f(x))/(x0-x);
x=x0;
x0= x0 - f(x0)/r;
k=1;
e(k)=abs(f(x0));

while k<Nmax && e(k)<toll

    r=(f(x0)-f(x))/(x0-x);
    x=x0;
    x0= x0 - f(x0)/r;
    k=k+1;
    e(k)=abs(f(x0));
end
end

function [w,l] = POT(H,p,toll,Nmax)

y=ones(length(H),1);
k=1;
l=1;
A=H-p*eye(length(H));
e=1;

if p==0

while k<Nmax && e>toll
    l0=l;
    w=A*y;
    [~,m]=max(w);

    l=y(m)/w(m);
    y=w/norm(w,inf);
    k=k+1;
    e=abs(l-l0);
end
else

while k<Nmax && e>toll
    l0=l;
    w=A\y;
    [~,m]=max(w);

    l=y(m)/w(m);
    y=w/norm(w,inf);
    k=k+1;
    e=abs(l-l0);
end

end

end

function [xc]=CHOL(A,b)

n=length(A);
S=zeros(n);

S(1,1)=sqrt(A(1,1));

for i=2:n
    for j=1:i-1
    S(i,j) = (A(i,j) - S(i,1:j-1)*S(j,1:j-1)')/S(j,j);
    S(i,i)= sqrt(A(i,i)-S(i,1:i-1)*S(i,1:i-1)');
    end
end

y=S\b;
xc=S'\y;

end