%esercizio 1 08.06.21

close all
clear all
clc

x1=[-2 -1 0 1 2 3];
y1=[7 5 1 -2 0 1];
n=4;

x2=x1;
y2=[1/3 -1/3 -1/3 1/3 5/3 11/3];

xx=linspace(x1(1),x1(end));
[Q,R] = QRC(x1,n);
[a] = SOLVE(Q,R,y1);
P1=@(xx) p1(a,xx);
P2=@(xx) Lagrange(xx,x2,y2);

figure
fplot(P1,[xx(1),xx(end)],'g')
hold on
fplot(P2,[xx(1),xx(end)],'r')
plot(x1,y1,'*b')
plot(x2,y2,'*c')

legend('App.min.quadrati per n=4','Metodo di Lagrange','(x1,y1)','(x2,y2)')
xlabel('-2 \leq x \leq 3')
title('Approssimazioni funzionali','FontAngle','italic')

%punto 3;
x3=linspace(0,1);
h=@(x3) P1(x3)+P2(x3) - 5*x3;

x0=0;
Nmax=100;
toll=sqrt(eps);
[x,e] = Secanti(h,x3,x0,Nmax,toll);

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

subplot(2,1,1)
hold on
plot(e,'b')
ylabel('e(k)')

title('Convergenza errore','Forma canonica','FontAngle','italic')
subplot(2,1,2);
plot(log10(e),'b');
title('','Forma logaritmica','FontAngle','italic')
xlabel('Iterazioni')
ylabel('log_1_0(e(k))')
function [Q,R] = QRC(x1,n)
m=length(x1);
n=n+1;
x1=x1';

C(:,[1,2])=[x1.^0 x1];

for i=1:m
    for j=3:n
        C(i,j)=C(i,j-1)*x1(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 [a] = SOLVE(Q,R,y1)
[~,n]=size(R);
b=Q*y1';
b=b(1:n);
R=R(1:n,1:n);
a=eye(n,1);
a(n)=b(n)/R(n,n);

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

end 

function [yy] = p1(a,xx)
yy=xx-xx;

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

end

function [zz] = Lagrange(xx,x2,y2)
n=length(x2)-1;
zz=xx-xx;

for j=0:n
    num=1;
    den=1;
    for i=[0:j-1,j+1:n]
        num=num.*(xx-x2(i+1));
        den=den.*(x2(j+1)-x2(i+1));
    end
    zz=zz + y2(j+1).*num./den;
end
end

function [x0,e] = Secanti(h,x3,x0,Nmax,toll)
x=x3(end);
k=1;
e(k)=abs(h(x0));

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