%es 1 21.01.20
clear all
close all
clc

x=[-2 -1 0 1 2 3];
y=[1.22 0.62 0.5 0.62 1.22 3.02];
n=2;
[Q,R] = QR2(x,n);
[a0] = SolveQR(Q,R,y);
xx=linspace(x(1),x(end));
PQ=@(xx) P0(a0,xx);
L=@(xx) Lagrange(xx,x,y);
f1=figure;
plot(x,y,'*r')
hold on 
fplot(PQ,[xx(1),xx(end)],'b')
fplot(L,[xx(1),xx(end)],'g')
legend('(x,f(x))','Approssimazione ai mininimi quadrati','Lagrange','Location','northwest')
title('Approssimazione funzionale', 'FontAngle','italic')

h=@(xx) L(xx) - 1;
a=-2;
b=-1;
Nmax=100;
toll=sqrt(eps);
[x0,e0,k0] = DIFFERENZEFINITE(h,a,b,Nmax,toll);
[x1,e1,k1] = SECANTI(h,a,b,Nmax,toll);

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

subplot(2,1,1)
hold on
fplot(h,'r')
title('Ricerca 0')
subplot(2,1,2)
fplot(h,'r',[a,b])
hold on
plot(x0,h(x0),'*k')
plot(x1,h(x1),'*b') %coincidono 
title('Zoom su intervallo di ricerca')

K0=linspace(1,k0,k0);
K1=linspace(1,k1,k1);
k=max(k0,k1);

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

subplot(2,1,1)
plot(K0,e0(K0),'g')
hold on
plot(K1,e1(K1),'m')
title('Convergenza Errore','Andamento canonico')
legend('Errore Newton','Errore Secanti')
subplot(2,1,2)
plot(K0,log10(e0(K0)),'g')
hold on
plot(K1,log10(e1(K1)),'m')
title(' ','Andamento logartimico')

%La velocità di convergenza di Newton è maggiore di quella delle secanti;

function [Q,R] = QR2(x,n)
m=length(x);
n=n+1;
x=x';

C(:,[1,2])=[x.^0 x];
for i=1:m
    for j=3:n
        C(i,j)=C(i,j-1)*x(i);
    end
end

Q=eye(m);

for i=1:n
    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 [a0] = SolveQR(Q,R,y)
[~,n]=size(R);

b=Q*y';
b=b(1:n);
R=R(1:n,1:n);
a0=eye(n,1);
a0(n)=b(n)/R(n,n);

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

function [y0] = P0(a0,xx)

y0=xx-xx;

for i=1:length(a0)
    y0=y0 + a0(i)*xx.^(i-1);
end
end

function [y1] = Lagrange(xx,x,y)
n=length(x)-1;
y1=xx-xx;

for j=0:n
    num=1;
    den=1;
    for i=[0:j-1,j+1:n]
        num=num.*(xx-x(i+1));
        den=den.*(x(j+1)-x(i+1));
    end
    y1=y1+y(j+1).*num./den;
end

end

function [x0,e0,k0] = DIFFERENZEFINITE(h,a,b,Nmax,toll)
x0=a;
k0=1;
e0(k0)=abs(h(x0));

while k0<Nmax && e0(k0)>toll
    if ~x0
        d=toll;
    else
        d=toll*x0;
    end
    r= (h(x0+d)-h(x0))/d;
    x0= x0 - h(x0)/r;
    k0=k0+1;
    e0(k0)=abs(h(x0));

end

end


function [x1,e1,k1] = SECANTI(h,a,b,Nmax,toll)
x1=a;
x=b;
r=(h(x1)-h(x))/(x1-x);
x=x1;
x1=x1-h(x1)/r;
k1=1;
e1(k1)=abs(h(x1));


while k1<Nmax && e1(k1)>toll
    r= (h(x1)-h(x))/(x1-x);
    x=x1;
    x1= x1 - h(x1)/r;
    k1=k1+1;
    e1(k1)=abs(h(x1));
end
end