%esercizio 2 18.02.22
clear all
close all
clc

x=linspace(-3,3, 8);
f=@(x) 1./(x.^2 +3*x +4);
y=f(x);
n=4;

[C] = SISTNORM(x,n);
%Valutazione con fattorizzazione QR
[Q,R] = QRC(C);
R1= R(1:n+1,1:n+1);
d=Q*y';
d1=d(1:n+1);
aqr=R1\d1;

%Valutazione con Cholesky
[S,b] =CHOL_C(C,y);
[l]=SolveTrilCH(S,b);
[a]=SolveTriuCH(S',l);


g=@(x) p(a,x);

fplot(g,'g')
hold on
fplot(f,'b')
plot(x,y,'*r')

legend('Approssimazione polinomiale', 'Funzione f','nodi-immagine')
title ('Analisi Approssimazione')
xlabel('Valori della x')
ylabel('f(x) & p(x)')

a=-3;
b=3;
NP=100;
xx=linspace(a,b,NP);

h= @(xx) (f(xx) - p(a,xx)).^2;
Ie=integral(h,a,b); %integrazione esatta

[It]=Trapez(h,xx);

function [C] = SISTNORM(x,n)
%definisco sistema normale per grado 4;

m=length(x);
C=ones(m,n);
C(:,2)=x';

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

function [Q,R] = QRC(C)
%decomposizione QR, Q matrix ortogonale e simm, R trapezoidale superiore

[m,r]=size(C);
Q=eye(m);
v=eye(m,1);

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

    if round(v(1),8)==0
        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);
        Qi= eye(m) -2*(u*u')/(u'*u);
        Q=Qi*Q;
        C=Qi*C;
    end
end
R=C;
end

function [S,b] = CHOL_C(C,y)
% devo risolvere C*C'*a=C'*y
[m,n]=size(C);

b=C'*y';
A=C'*C;

S=zeros(n);
S(1,1)=sqrt(A(1,1));

for i=2:n
    for j=1:i-1
        S(i,j)=(A(i,j) - S(i,1:j-1)*S(j,1:j-1)')/S(j,j);
        S(i,i)=sqrt(A(i,i) - S(i,1:i-1)*S(i,1:i-1)');
    end
end

end


function [l]=SolveTrilCH(S,b)
% S*l=b

n=length(S);
l=eye(n,1);
l(1)=b(1)/S(1,1);
for i=2:n
    l(i)=(b(i) - S(i,1:i-1)*l(1:i-1))/S(i,i);
end
end

function [a]=SolveTriuCH(S,l)
%S'a=l

n=length(S);
a=eye(n,1);
a(n)=l(n)/S(n,n);

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


function [px] = p(a,x)
% creazione pol

n=length(a);
px=0;

for i=0:n-1
    px=px +a(i+1)*x.^i;
end

end

function [It] = Trapez(h,xx)
% Formula dei trapezi composita
w=ones(size(xx))*(xx(2)-xx(1));
w(1)=w(1)/2;
w(end)=w(end)/2;

H=h(xx);
It=w*H';
end