% interpol07.m
% interpolation de Lagrange 
% fd, agreg, 08 octobre 2006, 04 mars 2025 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   clear , format long
   
printf("interpol07 : interpolation aux points de Tchebycheff\n");

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

xmin=-5 ; deltax = 2*abs(xmin) ; xmax = xmin + deltax ; 
alfa = 1 ; 
% cv facile :    xmin=-1 , alfa = 1 ;  
% cv difficile : xmin=-2 , alfa = .5 ;  
% demo : xmin=-5 ; alfa = 1 ;  
 
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 ;
%%% axis ( [ xmin-.2 , xmax+.2 , -.1 , 1/(alfa**2) + .3  ]) ; 
%%% axis ( [ xmin-.2 , xmax+.2 , -6 , 1/(alfa**2) + 3  ]) ; % cas qui diverge
plot(xex,yex,"k;fonction de Lorentz;") ; 
axis ( [ -5.1 , 5.1 , -.21 , 1.11  ] )
hold on ; 
%
kmax=8 ; nmax=2*kmax ;
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 = .5* deltax * cos( ((2*j-1)*pi) / (2*n+2))  ;  
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 ; 
%                        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 ;  

% 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 ,"c;n=2;") ;  pause (.5) ; end ; 
if (n==4)  ; plot(xx,yy ,"b;n=4;") ;  pause (.5)  ; end ; 
if (n==6)  ; plot(xx,yy ,"r;n=6;") ;  pause (.5)  ; end ; 
if (n==8)  ; plot(xx,yy ,"m;n=8;") ;  pause (.5)  ; end ; 
if (n==10)  ; plot(xx,yy ,"c;n=10;") ; pause (.5)  ; end ; 
if (n==12)  ; plot(xx,yy ,"k;n=12;") ; pause (.5)  ; end ; 
if (n==14)  ; plot(xx,yy ,"b;n=14;") ;  pause (.5)  ; end ; 
if (n==16)  ; plot(xx,yy ,"r;n=16;") ;  pause (.5)  ; end ; 
if (n==18)  ; plot(xx,yy ,"m;n=18;") ;  pause (.5)  ; end ; 
if (n==20)  ; plot(xx,yy ,"c;n=20;") ;  pause (.5)  ; end ; 
if (n==22)  ; plot(xx,yy ,"w;n=22;") ;  pause (.5)  ; end ; 
if (n==24)  ; plot(xx,yy ,"b;n=24;") ;  end ; 
end ; 

title ('interpolations de la Lorentzienne aux points de Tchbycheff') ;
print('07a-points-de-Tchebycheff', '-djpg','-S891,630');
hold off 


pause (2) ; 
axis ( [ 1.8 , nmax+.2 ,  -3 , 1  ]) ; % convergence difficile 
% axis ( [ 1.8 ,  nmax+.2 ,  -9 , 0  ]) ; % convergence facile 
xlabel (' interpol07. Nombre de noeuds d interpolation.') ;
ylabel (' logarithme en base 10 de l erreur en norme L infini ') ;
plot( xexp  ,  zexp  ,"*k;erreur mesuree;" )
axis ( [ 1.8 , nmax+.2 ,  -1.6 , -.17 ]) ; % convergence difficile 
title ('erreurs d interpolations de la Lorentzienne sur grilles de Tchebycheff') ;
print('07b-points-de-Tchebycheff', '-djpg','-S891,630');


%  calcul du coefficient avec une droite de regression
printf("premiere estimation de l ordre de convergence\n");
alpha=-codre(kmax,xexp,zexp)
