%es 1 29.06.20

close all
clear all
clc

a=-2;
b=2;
NV=[8 10 12 14];
k=1;
x=linspace(a,b);
f=@(x) 1./(1+10*x.^2);

for n=NV
%come P(x) costruisco polinomio interpolante nella forma di
%Lagrange(grado=4), P1(x) sarà sempre in forma di Lagrange,ma come nodi
%utilizzerò i nodi chebishev;

x1=linspace(a,b,n);
y1=f(x1);

if n==12
    xx=x1;
    yy=y1;
end

for i=0:n-1
    x2(i+1) = (a+b)/2 +(b-a)/2*cos((2*i+1)/(2*(n-1)+2)*pi);
end
y2=f(x2);

P=@(x) Lagrange(x,x1,y1);
P1=@(x) Lagrange(x,x2,y2);

h=@(x) abs(f(x)-P(x));
h1=@(x) abs(f(x)-P1(x));

[e(k)] = TCOMP(h,x);
[e1(k)] = TCOMP(h1,x);
k=k+1;

end

figure
plot(NV,e,'r')
hold on
plot(NV,e1,'g')

temp=xx;
xx=x;
x=temp;
y=yy;

g=@(xx) s(xx,x,y);
figure
fplot(g,[a,b])
hold on
plot(x,y,'*')

fplot(g,[a,b])

function [y] = Lagrange(x,xi,yi)
n=length(xi)-1;
y=x-x;
for j=0:n
    num=1;
    den=1;
    for i=[0:j-1,j+1:n]
        num=num.*(x-xi(i+1));
        den=den.*(xi(j+1)-xi(i+1));
    end
    y=y + yi(j+1).*num./den;
end
end

function [I] = TCOMP(h,x)
w=ones(1,length(x))*(x(2)-x(1));
w(1)=w(1)/2;
w(end)=w(end)/2;
H=h(x);

I=w*H';

end

function [yy]= s(xx,x,y)
x=x';
y=y';

h=diff(x);
v=diff(y)./h;

Dd= 2*(x(3:end)-x(1:end-2));
Cd=h(2:end-1);

A=diag(Dd) +diag(Cd,1)+diag(Cd,-1);
b=6*diff(v);

M=A\b;
M=[0;M;0];

yy=zeros(size(xx));

for i=1:length(x)-1
    I=(x(i)<=xx)&(xx<=x(i+1));

    C=y(i+1)./h(i) -h(i)/6*M(i+1);
    D=y(i)/h(i) - h(i)/6*M(i);

    yy=yy+(M(i+1)/6/h(i).*(xx-x(i)).^3 +M(i)/6/h(i).*(x(i+1)-xx).^3 +C.*(xx-x(i))+D.*(x(i+1)-xx)).*I;

end

end