% interpol04.m
% interpolation de Lagrange 
% fd, agreg, 08 octobre 2006, 04 mars 2025 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   clear , format long
   
printf("interpol04 : introduction des differences divisees\n");

%  n  nombre de donnees brutes ; c est un nombre d intervalles...
nex = 100 ; dxex = 3*pi  / nex  ;
xex =zeros(nex+1 , 1) ; yex =zeros(nex+1 , 1) ; 
for j=1:nex+1 ; xx =  (j-1)*dxex ;   
xex (j:j , 1) = xx ; 
yex (j:j , 1) = sin(xx) ; end 
hold off ;
plot(xex,yex,"k;fonction sinus;") ;  
axis ( [ -.1 , 3*pi+.1 , -1.1 , 1.4  ]) ;

n=16; 
n=22; 
n=30; 
x =zeros(n+1 , 1) ; 
u0=zeros(n+1 , 1) ; 
dx =  3*pi  / n ;
for j=1:n+1 ; xx = (j-1)*dx ;  
x (j:j , 1) = xx ;
u0(j:j , 1) =  sin(xx) ; end 
hold on ;
plot(x,u0(:,1),"*-b;interpolation affine par morceaux;") ;  
%   pause ; 

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


nn=300 ;  dy=3*pi/nn ; poly=zeros(nn+1,1) ; 
xx =zeros(nn+1, 1) ; 
for j=1:nn+1 ; xxx =  (j-1)*dy ;
xx (j:j , 1) = xxx ; end ;  
for k=1:nn+1 ; 
poly(k:k,1)=evalpol(n,newt,x,xx(k)) ;
end ; 

%             utilisation des fonctions usuelles de matlab  
[P, S] = polyfit (x, u0, n) ; 
nn=300 ;  dy=3*pi/nn ; 
xx =zeros(nn+1, 1) ; 
for j=1:nn+1 ; xxx =  (j-1)*dy ;
xx (j:j , 1) = xxx ; end ;  

yy=polyval(P,xx) ; 


hold on ;   
%%%  plot(xx,yy ,"r;polyfit et polyval : interpolation de Lagrange;") ; 
plot(xx,poly ,"r;interpolation de Lagrange avec Horner;") ; 
title ('interpolation du sinus de degre 30') ;
print('04-algorithme-Horner', '-djpg','-S891,630');
hold off 


printf("warning: matrix singular to machine precision, rcond = 1.75424e-39\n");
printf("warning: called from polyfit at line 121 column 5\n");
