// comparaison des interpolation polynomiales avec 
// points uniformement repartis
// points de chebyshev

function y=f(x)
   y = sin(x)
  // y = exp(x)
  // y = 1 ./ (1+x.*x)
endfunction;
  
function [x,pdet]=p(n,I,t,m) // degre n du polynome, intervalle I, t vecteur, 
                         //m=0 uniforme, m=1 chebyshev
  n1 = n+1 // nb de points
  a = I(1) // extremite gauche
  b = I(2) // extremite droite I=[a,b]
  
  pdet = zeros(t);
  if m==0
    x = a+(0:n).*((b-a)/n) // les points x_i uniformement repartis
  else
    x = (a+b)/2+(b-a)/2.*cos( (2 .* (0:n) +1).*(%pi/(2*n1))) // chebyshev
  end  
  
  for i = 1:n1
    lidet = ones(t)
    for j = 1:n1
      if (j<>i) 
        lidet = lidet .* ( (t-x(j))./(x(i)-x(j)))
      end
    end // on sort avec le vecteur li(t) pour les t points a afficher
    pdet = pdet + f(x(i)) .* lidet 
  end  
endfunction


t=-10:0.1:10; 


clf();
plot(t,f(t),'k-');

[xi,pt]=p(8,[-6,6],t,0);
plot(xi,f(xi),'bx',t,pt,'b-');
[xi,pt]=p(8,[-6,6],t,1);
plot(xi,f(xi),'gx',t,pt,'g-');


replot(rect=[-10 -5 10 5]);