% interpol06.m
% interpolation de Lagrange sur une grille reguliere
% fd, agreg, 08 octobre 2006, 04 mars 2025 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   clear , format long
   
printf("interpol06 : phenomene de Runge\n");

%  n  nombre de donnees brutes ; c est un nombre d intervalles...

xmin=-5 ; deltax = 2*abs(xmin) ; xmax = xmin + deltax ; 
alfa = 1 ; 
% alfa = .4 ; 
 
nex = 300 ; dxex = deltax  / nex  ;
xex =zeros(nex+1 , 1) ; yex =zeros(nex+1 , 1) ; 
for j=1:nex+1 ; xx = xmin  +  (j-1)*dxex ;   
xex (j:j , 1) = xx ; 
yex (j:j , 1) = 1 / (alfa**2 +xx**2 ) ; end 
hold off ;
%%%  if (xmin==-1) ; axis ( [ xmin-.2 , xmax+.2 , -.1 , 1/(alfa**2) + .3  ]) ; end ; 
%%%  if (xmin==-5) ; axis ( [ xmin-.2 , xmax+.2 , -1.1 , 1.3  ]) ; end ; 
plot(xex,yex,"k;fonction de Lorentz;") ;
axis ( [ -5.1 , 5.1 , -1.1 , 1.3  ] )
hold on ; 


kmax = 5
for k=1:kmax ; 
n = 2 * k ; 
x = zeros(n+1 , 1) ;  u0 = zeros(n+1 , 1) ; 
dx = deltax / n ;
for j=1:n+1 ; xx = xmin + (j-1)*dx ;  
x (j:j , 1) = xx ;
u0(j:j , 1) = 1 / (alfa**2 + xx**2 ) ;
end ; 
%%%%   if (n==2) 
%%%%   plot(x,u0(:,1),"*;donnees;") ; 
%%%%   else ;  plot(x,u0(:,1),"*;;") ;  end ; 
%   pause ; 
%             utilisation des fonctions usuelles de matlab  
%   [P, S] = polyfit (x, u0, n) ; 

%                        interpolation de Lagrange 'home made' 
newt=diffdiv(n,x,u0) ; 

%
nn=nex ;  dy=deltax/nn ; 
xx =zeros(nn+1, 1) ;  yy =zeros(nn+1, 1) ; 
for j=1:nn+1 ; xxx = xmin + (j-1)*dy ;
xx (j:j , 1) = xxx ; 
yy(j:j,1) = evalpol(n,newt,x,xxx) ;
end ;  
% yy=polyval(P,xx) ;


% calcul de l erreur d interpolation en norme l infini
err = 0; 
for j=1:nex+1 ; 
err = max ( err , abs ( yex (j:j , 1) - yy (j:j , 1) ) ) ; 
end ;  n ; err ;
xexp (k:k,1) = n ;  
zexp (k:k,1) = log ( err ) / log(10) ; 
 
hold on ; 
if (n==2)  ; plot(xx,yy ,"-b;n=2;") ; pause (.5)  ; end ; 
if (n==4)  ; plot(xx,yy ,"-r;n=4;") ;  pause (.5) ; end ; 
if (n==6)  ; plot(xx,yy ,"-m;n=6;") ;  pause (.5) ; end ; 
if (n==8)  ; plot(xx,yy ,"-k;n=8;") ;  pause (.5) ; end ; 
if (n==10)  ; plot(xx,yy ,"-g;n=10;") ; pause (.5) ; end ; 
if (n==12)  ; plot(xx,yy ,"-k;n=12;") ; end ; 
end ; 

title ('interpolations de la Lorentzienne sur grilles regulieres') ;
print('06a-phenomene-de-Runge', '-djpg','-S891,630');
hold off 


pause (2) 
%%%   if (xmin==-1) ; axis ( [ 1.8 , 12.2 ,  -4.2 , .2  ]) ;  end ; 
%%%   if (xmin==-5) ; axis ( [ 1.8 , 12.2 ,  -1.1 , 1.5  ]) ;  end ; 
%% xlabel (' interpol06. Phenomene de Runge. Nombre de noeuds equidistants. ') ;
%% ylabel (' logarithme en base 10 de l erreur en norme L infini ') ;
plot( xexp  ,  zexp  ,"*k;erreur mesuree;" )
axis ( [ 1.8 , 10.2, -.41 , .38  ] )
title ('erreurs d interpolations de la Lorentzienne sur grilles regulieres') ;
print('06b-phenomene-de-Runge', '-djpg','-S891,630');
