%   chaleur05.m
%   fd, 5 decembre 2010, 06 mars 2025  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   clear ; format long

fprintf('chaleur05 \n') ; 
fprintf('la chaleur avec les differences finies \n') 
fprintf('cas ou  le schema explicite en temps est stable \n')
fprintf('mesure de l ordre de convergence \n')

nmail  = 4 %%%  nombre de maillages 

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


ndf  = 32  , nitdf  = 120  

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 

hold off
%%%  plot (xxdf,y1df,'k;second membre;' ) 
%%%  axis ([ -.02 , 1.02 , -.11 , 1.05 ]) ;
%%%  title ('second membre') ;
hold off

%%% fprintf('on suppose mu=1 \n')
  sigma = .49 ; %% .70 %% .49    

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


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

yyt = zeros (ndf+1,1) ; 

for iter = 1:nitdf

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

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

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


uu=vv ; 
iter ;

end %%%  fin de la boucle en temps 

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

yex  = xxdf ; tt = nitdf * dtdf ; 

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%        dernier graphe   

for j=1:ndf-1 ;  yydf(j+1) = uu(j) ;   end 
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, 120 pas de temps, sigma = 0.49') ;
print('05a-schema-explicite' , '-djpg','-S891,630'); end 
if kmail==2 ; 
title ('64 mailles, 480 pas de temps, sigma = 0.49') ;
print('05b-schema-explicite' , '-djpg','-S891,630'); end 
if kmail==3 ; 
title ('128 mailles, 1920 pas de temps, sigma = 0.49') ;
print('05c-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 

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

fprintf('le temps tres long d execution de ce programme peut etre reduit \n')
fprintf('avec une meilleure programmation ! \n\n')

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