% Interpolazione polinomiale di una funzione trigonometrica
clear
figure(1)
clf
n = 4;
x = linspace(-pi,pi,n+1);
y = sin(x);
X = linspace(-pi,pi);
Y = sin(X);
p = polyfit(x,y,n);
Zp = polyval(p,X);
plot(x,y,'-o',X,Y,X,Zp,'--r')
pause

n1 = 8;
x1 = linspace(-pi,pi,n1+1);
y1 = sin(x1);
p1 = polyfit(x1,y1,n1);
Zp1 = polyval(p1,X);
plot(x1,y1,'-*',X,Y,X,Zp1,'--r')
pause
% Altra funzione
figure(2)
clf
n = 16;
x = linspace(-1,1,n+1);
y = 1./(1+25*x.^2);
X = linspace(-1,1,1000);
Y = 1./(1+25*X.^2);
plot(X,Y,x,y,'o')
legend('Cauchy','Nodi')
p = polyfit(x,y,n);
Yp = polyval(p,X);
hold on
plot(X,Yp,'r');
pause
% Calcolo del polinomio omega
omega = (X-x(1));
for k = 2:n+1
    omega = omega.*(X-x(k));
end

% Scegliamo i nodi di Chebyshev
I = 1:n+1;
phi = (2*I-1)/(2*(n+1))*pi;
xc = -cos(phi);
omegac = (X-xc(1));
for k = 2:n+1
    omegac = omegac.*(X-xc(k));
end
figure(3)
clf
plot(X,omega,X,omegac,'r')
pause
axis([-1 1 -2e-2 2e-2])
hold off
yc = 1./(1+25*xc.^2);
pc = polyfit(xc,yc,n);
Yc = polyval(pc,X);
figure(2)
hold on
plot(xc,yc,'*',X,Yc,'g')