% lu2.m
% fd, iacs, 04 septembre  2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% title ('Cnam Paris. Travaux pratiques Iacs.') ;
xlabel (' Interpolation par splines cubiques  ') ;
'Factorisation L*U d une matrice tridiagonale' , 
'Puis resolution des systemes lineaires tridiagonaux' ,
'Puis affichage de la spline cubique ' ,  pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=5
%   reservation de place memoire pour stocker les vecteurs
a=zeros(n+1,1) ; b=zeros(n,1) ; c=zeros(n,1);
al=zeros(n+1,1); be=zeros(n,1); ga=zeros(n,1);
%   donnees pour les splines cubiques
a(1:1,1)=2 ; 
for i=2:n; a(i:i,1)=4 ; end ; 
a(n+1:n+1)=2 ; 
for i=1:n; b(i:i,1)=1 ;  c(i:i,1)=1 ; end ;      
%   algorithme L*U
al(1:1,1)=a(1:1,1) ; 
for i=1:n ;
be(i:i,1)=b(i:i,1) ; 
ga(i:i,1)=c(i:i,1)/ al(i:i,1) ; 
al(i+1:i+1,1)=a(i+1:i+1,1) - ga(i:i,1)*be(i:i,1) ; 
end ;     
%   donnees a interpoler
% pdebut = [-1 1 0] ; 
pdebut = [-1 1 -1 1  0] ; 
x =zeros(n+1 , 1) ;  u0=zeros(n+1 , 1) ;  dx = 1 / n ;
for j=1:n+1 ; xx = (j-1)*dx ;   x (j:j , 1) = xx ;
u0(j:j , 1) = polyval(pdebut,xx) ; end ;
u0
%   formation du second membre
'Attention. On se trompe toujours dans le second membre' , pause
sm =zeros(n+1,1) ;  
sm(1:1,1)=3*(u0(2:2,1)-u0(1:1,1)) ; 
for j=2:n ;  sm(j:j,1)=3*(u0(j+1:j+1,1)-u0(j-1:j-1,1)) ; end
sm(n+1:n+1,1)=3*(u0(n+1:n+1,1)-u0(n:n,1)) ; 
%   descente
yy =zeros(n+1,1) ; 
yy(1:1,1)=sm(1:1,1) ; 
for j=1:n ;  
yy(j+1:j+1,1)=sm(j+1:j+1,1) - ga(j:j,1)*yy(j:j,1) ; 
end ; 
%   remontee
xx =zeros(n+1,1) ;
xx(n+1:n+1,1)=yy(n+1:n+1,1) / al(n+1:n+1,1) ; 
for i=1:n ; j=n+1-i ; 
xx(j:j,1)=(yy(j:j,1) - be(j:j,1)*xx(j+1:j+1,1)) /  al(j:j,1) ;
end ;
xx
%   affichage sur une seconde grille plus fine
nn=22
dy=1/nn
xxx =zeros(nn+1,1);
uu1 =zeros(nn+1,1); uu2 =zeros(nn+1,1); uu3 =zeros(nn+1,1); uu4 =zeros(nn+1,1);
for i=1:nn+1;
xxx(i) = (i-1) * dy ;  
uu1(i) = (1+2*xxx(i))*(xxx(i)-1)*(xxx(i)-1) ; 
uu2(i) = (3-2*xxx(i))*xxx(i)*xxx(i)  ;  
uu3(i) =  xxx(i)*(xxx(i)-1)*(xxx(i)-1)  ; 
uu4(i) =  xxx(i)*xxx(i)*(xxx(i)-1)  ;  end ; 
m=n*nn ,  xxxx=zeros(m+1,1) ;  uuuu=zeros(m+1,1) ; 
dxxxx=1/m ; k=1 ;  xxxx(k:k,1)=0 ; uuuu(k:k,1)=u0(1:1,1) ; 
for j=1:n;
for i=2:nn+1;
k=k+1 ; xxxx(k:k,1)=(k-1)*dxxxx ; 
uuuu(k:k,1)=u0(j:j,1)*uu1(i)+u0(j+1:j+1,1)*uu2(i)+ \
            xx(j:j,1)*uu3(i)+xx(j+1:j+1,1)*uu4(i);
end;end;
k ,   xxxx(k:k,1) ,  uuuu(k:k,1)
% gset yrange[-0.1:0.4] ; gset xrange[-0.1:1.1] ; 
hold off ; 
axis ( [ -.1 , 1.1 , -.1 , .4 ]) ;
plot(x,u0(:,1),"*;donnees;",x,u0(:,1),"b;interpolation affine;",\
% plot(x,u0(:,1),"*;donnees;",\
xxxx,uuuu(:,1),"g;spline cubique;") ;  





