%esercizio 1 22.06.22

clear all
close all
clc

x=[-2,-1,0,1,2,3];
y=[5,2,0.5,0,-1,-4];

plot(x,y,'*b')

n=3; %il grado minimo per averlo interpolante è 4
[Q,R]=QRC(x,n);
[a]=SOLVE(Q,R,y);

c=x(1);
d=x(end);
xx=linspace(c,d);

f=@(xx) PolVal(a,xx);

fplot(f,'g',[c,d])
hold on
plot(x,y,'*r')

legend('Approssimazione','nodi-immagine','Location','southwest')

g=@(xx) 1+5*xx;

h=@(xx) f(xx) -g(xx);

I=linspace(-1,1);
Nmax=100;
toll=sqrt(eps);
[x0] = SEC(h,I,Nmax,toll);

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

figure(f1)
subplot(2,1,1)
hold on
fplot(h,'r')
legend('h(x)')
title('Ricerca 0','FontAngle','italic','FontName','Times','FontSize',13)

subplot(2,1,2)
hold on
fplot(h,'r',[I(1),I(end)])
plot(x0,h(x0),'*g')
yline(0,'--k')

legend('h(x) in I','(x0,y0)','y=0')


xlabel('-1 \leq x \leq 1')
ylabel('h(x)')

function [Q,R]=QRC(x,n)
x=x';
m=length(x);

C=eye(m,n+1);
C(:,[1,2])=[x.^0, x];

for i=1:m
    for j=2:n
        C(i,j+1)=C(i,j)*x(i);
    end
end

Q=eye(m);

for i=1:n+1
    v=C(i:m,i);
    u=zeros(m,1);

    if v(1)==0
        u(i:m)=v + norm(v)*eye(m-i+1,1);
    else
        u(i:m)=v + sign(v(1))*norm(v)*eye(m-i+1,1);
    end

    Qi=eye(m) -2*(u*u')/(u'*u);

    Q=Qi*Q;
    C=Qi*C;
end
R=C;

end

function [a]=SOLVE(Q,R,y)
[~,n]=size(R);

y=Q*y';
y=y(1:n);
R=R(1:n,1:n);

a=eye(n,1);
a(n)=y(n)/R(n,n);

for i=n-1:-1:1
    a(i)=(y(i)-R(i,i+1:n)*a(i+1:n))/R(i,i);
end
end

function [p] = PolVal(a,xx)
p=xx-xx;
n=length(a);

for i=1:n
    p= p +a(i)*xx.^(i-1);
end
end

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

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