%   chaleur07.m
%   fd, 5 decembre 2010, 07 mars 2025  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   clear ; format long

fprintf('chaleur07 \n') ; 
fprintf('equation de la chaleur avec des differences finies \n') 
fprintf('schema implicite en temps \n')
fprintf('mesure de l ordre de convergence \n')


hold off 

sigma = 10
nmail  = 4 %%%  nombre de maillages 

nnn  = zeros (nmail,1) ; erreurs = zeros(nmail,1) ; 


ndf  = 32  ; nitdf  = 20 ; 

for kmail = 1:nmail ; %%%  boucle sur les maillages

ndf, nitdf   
dxdf = 1/ndf ; 

xxdf = zeros(1,ndf+1) ; yydf = xxdf ;  y1df = xxdf ;  yex  = xxdf ; 
xa    = 0 ; 
for j = 1:ndf+1 ; xxdf(j) = xa ; y1df(j) =  sin(pi*xa)  ; 
xa = xa + dxdf ; end 

%%%%  if kjob==1 ; 
%%%%  plot (xxdf,y1df,'k;;' ) 
%%%%  axis ([ -.02 , 1.02 , -.05 , 1.05 ]) ;
%%%%  pause (.5) 
%%%%  end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     boucle en temps 
 
dtdf  = sigma * dxdf**2                %  pas de temps         


uu = zeros (ndf-1,1) ; vv = uu ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%       matrice 
A = zeros (ndf-1,ndf-1) ;
for j=1:ndf-1 ; A(j,j) = 1+2*sigma ; end 
for j=1:ndf-2 ; A(j,j+1) = - sigma ; end 
for j=2:ndf-1 ; A(j,j-1) = - sigma ; end 
A ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for iter = 1:nitdf

for j=1:ndf-1 ;  yydf(j+1) = uu(j) ;   end 

%%%%%   jmaj = 1  ;  its10 = floor(iter/jmaj); 
%%%%%   if (jmaj*its10==iter)
%%%%%   plot ( xxdf,y1df,'r;inf;' ,   xxdf,yydf,'b*;t;' ) 
%%%%%   pause  (.5)
%%%%%   end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%       second membre  
sm = zeros (ndf-1,1) ;

for j=1:ndf-1 ; 
sm(j) =  uu(j) +  sigma * y1df(j+1) * (pi*dxdf)**2 ;
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%       resolution du systeme lineaire   
vv = A \ sm ; uu=vv ;

end    %%%%%%%%       fin de la boucle en temps


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%        solution exacte avec Fourier  

tt = nitdf * dtdf ; 

for j=1:ndf+1
yex(j)  =  sin(pi*xxdf(j)) * (1 - exp(-tt*pi**2)) ;    
end 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%       graphe   



yyt = zeros (ndf+1,1) ; 
for j=2:ndf ; yyt(j) = vv(j-1) ; end
yyt ; 

plot ( xxdf,y1df,'k;second membre;' , xxdf,yydf,'*-r;differences finies;' , ...
       xxdf,yex,'b;solution exacte;')
axis ([ -.02 , 1.02 , -.05 , 1.05 ]) ;
if kmail==1 ; 
title ('32 mailles, 20 pas de temps, sigma = 10') ;
print('07a-schema-explicite' , '-djpg','-S891,630'); end 
if kmail==2 ; 
title ('64 mailles, 80 pas de temps, sigma = 10') ;
print('07b-schema-explicite' , '-djpg','-S891,630'); end 
if kmail==3 ; 
title ('128 mailles, 320 pas de temps, sigma = 10') ;
print('07c-schema-explicite' , '-djpg','-S891,630'); end 
pause(1) ; hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%        calcul de l erreur

erlinf = max(abs(yydf-yex))

  nnn(kmail,1)       = log(ndf) / log(2)  ; 
erreur (kmail,1) = log(erlinf) / log(2) ;

ndf = 2*ndf ;   nitdf  = 4 * nitdf ;   
end  %%%%%%%%       fin de la boucle sur les maillages  

hold off 

nnn ;  erreur ; 
  plot (nnn(:,1) , erreur(:,1),'*k;erreur L infini;'  ) 
axis ([ 4.8 , 8.2, -11.5 , -4.8 ]) ;
title ('erreurs en norme L infini') ;
print('07d-convergence-schema-implicite' , '-djpg','-S891,630');


fprintf('ordre de convergence mesure  \n')
alpha = -codre (nmail , nnn , erreur) 
