%Esercizio 1 8/06/21

clear all
close all
clc

x1=[-2 -1 0 1 2 3];
y1=[7 5 1 -2 0 1];
x2=x1;
y2=[1/3 -1/3 -1/3 1/3 5/3 11/3];

n1=4;
n2=2;

% 1. Devo applicare l'approssimazione ai minimi quadrati alla prima serie e
%   newton o lagrange alla seconda;
%
% 2. Graficare le due funzioni coi dati iniziali;
%
% 3. Verificare quanti 0 ho per la funzione P1(x)+P2(x)-5x in I=[0,1];
%   Poi cercarli;
%
% 4. Graficare errore;
[Q,R] = QR(x1,y1,n1);
[a] = SOLVE(Q,R,y1);
xx=linspace(-2,3);

f1=@(xx)  P1(xx,a);
f2=@(xx)  P2(xx,x2,y2);
%fN=@(xx)  PN(xx,x2,y2,F);

fplot(f1,'r')
hold on
fplot(f2,'g')
%fplot(fN, 'b')

I=linspace(0,1);
h=@(I) f1(I) + f2(I) -5*I;
Nmax=100;
toll=sqrt(eps);

[x0,e,k] = SEARCH(h,I,Nmax,toll);
K=linspace(0,k, k+1);
fplot(h)
xline(x0)

figure
plot(K,e)
hold on
plot(K,e,'*')

function [Q,R] = QR(x1,y1,n1)

% C   : matrice del sistema normale;
% Q,R : matrici a seguito della decomposizione, per come è costruito
% l'algoritmo in realtà ho poi Q'*R=C;
x1=x1';
m=length(x1);

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

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

Q=eye(m);
n=n1+1;

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
a=a';

end

function [p1] = P1(xx,a)
l=length(xx);
n=length(a);
P=zeros(1,l);

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

function [p2] = P2(xx,x2,y2)
% metodo di Lagrange
n=length(x2)-1;
P=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
    P=P + y2(j+1).*num./den;
end
p2=P;
end
% 
% function [F] = DIFFDIV(x2,y2)
% 
% outputArg1 = inputArg1;
% outputArg2 = inputArg2;
% end
% 
% function [pn] = PN(xx,x2,y2,F)
% 
% outputArg1 = inputArg1;
% outputArg2 = inputArg2;
% end

function [x0,e,k] = SEARCH(h,I,Nmax,toll)
x0=I(1);
x1=I(end);
r=(h(x1)-h(x0))/(x1 -x0);
x1=x0;
x0=x0 - h(x0)/r;
k=1;
e(k)=abs(h(x0));

while k<Nmax && e(k)>toll
    r=(h(x1)-h(x0))/(x1 -x0);
    x1=x0;
    x0=x0 - h(x0)/r;

    k=k+1;
    e(k)=abs(h(x0));
end
k=k-1;
end