
///  Exercice 1.
//disp('Effacement des graphiques')
xdel(winsid());

//On appelle s la variable polynomiale.
s=poly(0,"s");

disp('Exercice 1');


function[P]=lagrange(x,f,v)
  // x est un vecteur contenant les points d'interpolations
  // f est un vecteur contenant les valeur de la fonction aux points
  // d'interpolation
  // v est un parametre sur lequel on veut evaluer le polynome d'interpolation.
  n=length(x);
  l=ones(1,n);
  for i=1:n
    J=[(1:i-1),(i+1:n)];
    l(J)=l(J).*(v-x(i))./(x(J)-x(i));
  end
  P=sum(f.*l);
endfunction

function[diff_div]=differences_divisees(x,f)
  diff_div=f;
  diff_aux=f;
  n=length(x);
  for k=2:1:n
    diff_aux=(diff_aux(2:$)-diff_aux(1:$-1))./(x(k:$)-x(1:$-k+1));
    diff_div(k)=diff_aux(1);
  end
endfunction

function[y]=serie_newton(x,diff_div,v)
  diff_aux=x;
  y=diff_div($);
  for i=length(diff_div)-1:-1:1
    y=diff_div(i)+(v-x(i))*y;
  end
endfunction

function[y]=newton(x,f,v)
  diff_div=differences_divisees(x,f);
  y=serie_newton(x,diff_div,v);
endfunction

function[res]=neville(x,f,v)
  n=length(x);
  x=ones(v)*x;res=ones(v)*f;
  v=v*ones(f);
  for k=1:n-1
    I=1:n-k;
    v=v(:,I);
//    disp([size(v),size(x),size(res)])
    res=((v-x(:,I)).*res(:,I+1)-res(:,I).*(v-x(:,I+k)))./(x(:,I+k)-x(:,I));
  end
endfunction

N=4;
disp('N=4, f(x)=exp()');
x=linspace(-1,1,N);
disp(lagrange(x,exp(x),s),'Lagrange: ');
disp(newton(x,exp(x),s),'Newton: ');
disp(neville(x,exp(x),s),'Neville: ');


N=80;u=log(2);
disp('N='+string(N)+', f(x)=exp(x), valeur du polynome d interpolation en '+string(u));
x=linspace(-1,1,N);

disp(exp(u),'exp('+string(u)+')');

disp(lagrange(x,exp(x),u),'Lagrange: ');
disp(newton(x,exp(x),u),'Newton: ');
disp(neville(x,exp(x),u),'Neville: ');

disp('Apres etre repasse dans la base standard')
disp(horner(lagrange(x,exp(x),s),u),'Lagrange: ');
disp(horner(newton(x,exp(x),s),u),'Newton: ');
disp(horner(neville(x,exp(x),s),u),'Neville: ');


disp('');
disp('');

disp('Observation: les series de Newton sont instables. En revanche, Lagrange et Neville se comportent mieux. Neville resiste au passage dans la base standard mais pas Lagrange. On utilisera donc Neville dans la suite.')

//Nouvelle fonction renvoyant le polynome d'intepolation
function[P]=lagrange_neville(x,f)
  P=neville(x,f,s);
endfunction


figure();
N=4;
x=linspace(-2,2,N);
P=lagrange_neville(x,exp(x));
x=linspace(-2,2,100);
plot2d(x,[horner(P,x);exp(x)]');
legend('$P_{'+string(N+1)+'}$','$\exp(x)$');
xlabel('$x$');
xtitle('$\exp(x) \mbox{ et son interpolee}$');



// Exercice 2
disp('Exercice 2');
function[y]=cheb(N)
  //Calcul des points de Chebyshev.
  I=0:1:N
  y=cos(%pi*(2*I+1)/(2*N+2))
endfunction

figure();
subplot(2,1,1)
N=8
x=linspace(-1,1,N+1);
Pequi=lagrange_neville(x,exp(x));
x=cheb(N);
Pcheb=lagrange_neville(x,exp(x));
x=linspace(-1,1,200);
plot2d(x,[abs(horner(Pequi,x)-exp(x));abs(horner(Pcheb,x)-exp(x))]');
legend(['Points equidistants','Points de Chebyshev']);
title('$\mbox{Erreur d interpolation pour } \exp(x)$');
xlabel('x');




subplot(2,1,2)
Nt=2:80; err=[];
for N=Nt
  v=linspace(-1,1,200)';
  x=linspace(-1,1,N);
  xcheb=cheb(N-1);
  err=[err,[norm(neville(x,exp(x),v)-exp(v),'inf');...
            norm(neville(xcheb,exp(xcheb),v)-exp(v),'inf')]];
end
plot2d(Nt,err',logflag='nl')
xtitle('Erreur d interpolation','$\mbox{Nombre de points de collocation:} N$','$\|f-P_N\|_{\infty}$')
legend(['Points equidistants','Points de Chebyshev']);
//On remarque que l'interpolee converge tres vite vers exp(x) mais
//que l instabilite prend ensuite le dessus pour les points equidistants.

///Exercice 3
disp('Exercice 3');

function[]=runge(alpha,pts,str)
Nt=[6,10,16];
Pval=[];
for N=Nt
  x=eval(pts);
  P=lagrange_neville(x,(x.^2+alpha^2).^(-1))
  x=linspace(-1,1,200);
  Pval=[Pval;horner(P,x)];
end
Pval=[(x.^2+alpha^2).^(-1);Pval];
plot2d(x,Pval');
xtitle('$\mbox{Interpolee de } \frac{1}{x^2+\alpha^2} \mbox{ lorsque } \alpha='+string(alpha)+'\mbox{ ('+str+')}$');
xlabel('x');
legend(['$\frac{1}{x^2+\alpha^2}$','$P_6(x)$','$P_{10}(x)$','$P_{16}(x)$'])
endfunction

figure();
subplot(2,2,1)
runge(1/sqrt(8),'linspace(-1,1,N+1)','Pts equidistants');
subplot(2,2,2)
runge(1/sqrt(8),'cheb(N)','Pts de Chebyshev');

subplot(2,2,3)
runge(sqrt(8),'linspace(-1,1,N+1)','Pts equidistants');
subplot(2,2,4)
runge(sqrt(8),'cheb(N)','Pts de Chebyshev');



function[u]=potentiel_equi(x,y)
  z=x+%i*y;
  u=-1+real((z+1)/2*log(z+1)-(z-1)/2*log(z-1));
endfunction


function[u]=potentiel_cheb(x,y)
  z=x+%i*y;
  u=abs(log((abs(z-sqrt(z^2-1)))))-log(2);
endfunction


function[]=runge_complexe(alpha,pts,str,potentiel)
N=20;
  x=eval(pts);
  P=lagrange_neville(x,(x.^2+alpha^2).^(-1))
  function[u]=ecart(x,y)
    z=x+y*%i;
    u=min(abs(horner(P,z)-1/(z^2+alpha^2)),1);
  endfunction  
  x=linspace(-1.5,1.5,100);
  y=linspace(-0.5-alpha,0.5+alpha,200);
  fgrayplot(x,y,ecart);
  colorbar(0,1)
  contour(x,y,potentiel,(-2:.1:2));
  plot(0,alpha,'*');
  plot(0,-alpha,'*');
  legend('$\mbox{Poles de } f_{\alpha}$')
xtitle('$\mbox{Erreur sur } \frac{1}{x^2+\alpha^2} \mbox{ lorsque } \alpha='+string(alpha)+'\mbox{ ('+str+')}$');
xlabel('$\Re z$');
ylabel('$\Im z$');
endfunction


figure();
subplot(2,2,1)
runge_complexe(1/sqrt(8),'linspace(-1,1,N+1)','Pts equidistants',potentiel_equi);
subplot(2,2,2)
runge_complexe(1/sqrt(8),'cheb(N)','Pts de Chebyshev',potentiel_cheb);

subplot(2,2,3)
runge_complexe(sqrt(8),'linspace(-1,1,N+1)','Pts equidistants',potentiel_equi);
subplot(2,2,4)
runge_complexe(sqrt(8),'cheb(N)','Pts de Chebyshev',potentiel_cheb);



/// Exercice 4
disp('Exercice 4')

function[y]=fourier_interp(f,N)
  //f: valeurs de f aux points de la grille.
  //N: nombre de points sur lequel on souhaite calculer l'interpolee de f.
  y=fftw(f,-1);
  k=ceil(length(f)/2);
  y=[y(1:k),zeros(1,N-2*k+1),y($-k+2:$)];
  y=fftw(y,1)*N/length(f);
endfunction

function[y]=ff(x)
  x=modulo(2*%pi+modulo(x,2*%pi),2*%pi);
  y=x-%pi;
  I=find(x<%pi);
  y(I)=%pi-x(I);
endfunction


N=11;
h=1/N;
fe=ff(2*%pi*(0:h:1-h));
finterp=fourier_interp(fe,1000);
N=1000;
h=1/N;
fe=ff(2*%pi*(0:h:1-h));
figure();
plot2d(2*(0:h:1-h)*%pi,real([fe;finterp]'))
xlabel('x')
legend('f','Interpolee de Fourier de f');
