%esercizio 1 21.01.20

clear all
close all
clc

x1=[-2 -1 0 1 2 3];
y1=[2.18 1.49 0.5 -0.49 -1.18 -0.07];
n=3;

[Q,R] = QR(x1,n);
[a] = Solve(Q,R,y1);

xx=linspace(-2,3);
P_1=@(xx)  P1(xx,a);
P_2=@(xx)  P2(xx,x1,y1);

fplot(P_1,[-2,3],'r')
hold on
fplot(P_2,[-2,3],'b')
legend('Approssimazione ai minimi quadrati','Metodo di Lagrange')
plot(x1,y1,'*r')
t=linspace(1,2);
f=@(t) P_2(t)+1;
figure
fplot(f,[1,2],'g')

Nmax=100;
toll=sqrt(eps);
[x0] = NewtonAlleDiff(f,t,Nmax,toll);
hold on
plot(x0,f(x0),'*r')


function [Q,R] = QR(x1,n)
x1=x1';
m=length(x1);

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

for i=1:m
    for j=3:n+1
        C(i,j)=C(i,j-1)*x1(i);
    end
end
n=n+1;

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)
%risoluzione sistema triangolare superiore

[m,n]=size(R);
R=R(1:n,1:n);
b=Q*y1';
b=b(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 [p1] = P1(xx,a)

n=length(a)-1;
p1 = xx-xx;

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

end

function [p2] = P2(xx,x1,y1)

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

end

function [x0] = NewtonAlleDiff(f,t,Nmax,toll)
x0=t(1);
x1=t(end);
r=(f(x0)-f(x1))/(x0-x1);
x1=x0;
x0= x0 -f(x0)/r;
k=1;
e(k)=abs(f(x0));

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

end