%es 2 18.02.20

clear all
close all
clc
f1=figure;
subplot(2,1,2)
x=linspace(-3,3,8);
f=@(x) 1./(x.^2 + 3*x +4);

n=4;
y=f(x);
[M,b] = NORMSIST(x,y,n);
[S] = CHOL(M);
[c] = SolveTril(S,b);
[a] = SolveTriu(S',c);

P=@(x) PA(a,x);
subplot(2,1,1)
hold on
fplot(f,'g')
fplot(P,'b')
plot(x,f(x),'*r')
title('Approssimazione funzionale')

subplot(2,1,2)
hold on
plot(x,f(x),'*r')
fplot(f,[x(1),x(end)],'g')
fplot(P,[x(1),x(end)],'b')
title('Zoom su intervallo d''interesse','FontAngle','italic','FontSize',10)
xlabel ('-3 \leq x \leq 3')
ylabel('f(x) , P(x)')
legend('(x,f(x))','funzione data','Approssimazione','Location','southwest')

x=linspace(-3,3,2000);
g=@(x) (f(x)-P(x)).^2;

figure
fplot(g,[-3,3])
colors=[0.3010 0.7512 0.9330];
area(g(x),'LineStyle',':')
colororder(colors)
legend('Area Integrale')
title('Stima integrale sotteso')

[I] = TRAPEZCOMP(g,x);

function [M,b] = NORMSIST(x,y,n)
%ho un sitema C'*C*a=C'*y --> M*a=b --> M simm, def positiva per struttura.
n=n+1;
m=length(x);
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

b=C'*y';
M=C'*C;

end

function [S] = CHOL(M)

n=length(M);
S=zeros(n,n);
S(1,1)=sqrt(M(1,1));

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

function [c] = SolveTril(S,b)
n=length(b);
c=eye(n,1);
c(1)=b(1)/S(1,1);

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

end


function [a] = SolveTriu(S1,c)
n=length(c);
a=eye(n,1);
a(n)=c(n)/S1(n,n);

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

end

function [p] = PA(a,x)
n=length(a);
p=x-x;

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

function [I] = TRAPEZCOMP(g,x)

w=ones(size(x))*(x(2)-x(1));
w(1)=w(1)/2;
w(end)=w(end)/2;

F=g(x);
I=w*F';

end