% dd3.m
% fd, 9, 10, 11, 13  novembre 2010, 25, 26 fevrier, 02, 08 mars 2025 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   clear, format long 


fprintf ('dd3 \n') ;
fprintf ('probleme de Dirichlet pour le Laplacien \n') ;
fprintf ('differences finies et volumes finis  a trois points \n') 
fprintf ('solution exacte = x (1-x) (1/2-x) \n') 
fprintf ('octave 6.4 \n')    

nbnn = 6 
nnn = zeros (nbnn,1) ; erreur = zeros(nbnn,4) ; 

n = 8  ; 
for lll = 1:nbnn  
n
dx = 1/n ; 

adf = zeros(n-1,n-1) ;          %  la matrice 'df' ou 'ef'
for j=1:n-1; 
adf(j,j) = 2 ; end 
for j=1:n-2 ;   adf(j,j+1) = -1 ;  end 
for j=2:n-1 ;   adf(j,j-1) = -1 ;  end
adf ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   scaling 
for i = 1:n-1 , for j=1:n-1 
adf(i,j) =  adf(i,j) / (dx**2) ; end , end 



avf = zeros(n,n) ;          %  la matrice 'vf' 
for j=1:n; 
avf(j,j) = 2 ; end 
for j=1:n-1 ;   avf(j,j+1) = -1 ;  end 
for j=2:n  ;   avf(j,j-1) = -1 ;  end
avf(1,1) = 3 ; avf (1,2) = -1 ;       %%%%%%%%%  premiere version 
avf(n,n) = 3 ; avf (n,n-1) = -1 ;     %%%%%%%%%  premiere version 
% avf(1,1) = 4 ; avf (1,2) = -4/3 ;       %%%%%%%%%  seconde version 
% avf(n,n) = 4 ; avf (n,n-1) = -4/3 ;     %%%%%%%%%  seconde version 
avf  ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   scaling 
for i = 1:n , for j=1:n 
avf(i,j) =  avf(i,j) / (dx**2) ; end , end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  second membre 
bdf = zeros(n-1,1) ;

xa    = dx ;  
for j = 1:n-1 ; 
bdf(j) = 3-6*xa ; xa = xa + dx ;    end 
bdf'   ; 


%%%%%%%%%       formule de Simpson pour les volumes finis. 

bvf = zeros(n,1) ; 
xm = 0 ; xa = dx/2 ;  xp = dx ;    
for  j = 1:n ; 
bvf(j) = 3 - 6*((xm/6) + (2*xa/3) + (xp/6))  ;   %%%%%%%%%%%%%%    Simpson 
% bvf(j) = 3 - 6*xa  ;       %%%%%%%%%%%%%%    point milieu
xm = xm + dx ; xa = xa + dx ; xp = xp + dx ;    end 
bvf'   ; 



udf = adf \ bdf ; 
uvf = avf \ bvf ; 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%   graphique   des champs  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 

xx = zeros(1,n+1) ; yy = xx ;  ye = xx ; 
for j=1:n+1 ; xx(j) = (j-1)*dx ; end  

xvf = zeros(1,n) ; yevf = xvf ;  yvf = xvf ; 
xa  = dx/2 ; 
for j=1:n ; 
xvf(j) = xa ;  yevf(j) = xa * (1-xa) * (1/2-xa) ; 
yvf(j) = uvf(j) ;     xa = xa + dx  ;  end  
 
for j = 1:n+1 ; 
ye(j) = xx(j)*(1-xx(j))*(1/2-xx(j)) ; end 
for j = 2:n 
yy(j) = udf(j-1) ; end 
yy(n+1)  = 0 ; 


%%%%%%%%%%%%%%%%%%%%%%  evaluation des erreurs 
%%%%%%%%%%%%%%%%%%%%%%  differences finies  

errli = 0 ;  errl2 = 0 ; 
for j = 1:n
tt    = abs(yy(j) - ye(j)) ; 
errli = max (errli , tt ) ; 
errl2 = errl2 + tt**2 ; end 
errl2 = sqrt ( dx * errl2 )  ; 
errli ; 

%%%%%%%%%%%%%%%%%%%%%%  volumes finis  


evfli = 0 ; evfl2 = 0 ; 
for j = 1:n
tt    = abs(yvf(j) - yevf(j)) ; 
evfli = max (evfli , tt ) ; 
evfl2 = evfl2 + tt**2 ; end 
evfl2 = sqrt ( dx * evfl2 ) ;  
evfli ; 


if (n==8) ; 
plot (xx,ye,'*k;solution analytique;' , xvf,yevf,'*k;;' , ... 
        xx,yy,'*r;differences finies;' , xvf , yvf , '*b;volumes finis;' ) 
 axis ( [ -0.01, 1.01 , -0.055  , 0.055 ]) ;
title ('Laplacien, Dirichlet homogene, second membre polynomial de degre 1, n=8') ;
print('01a-dd3-dirichlet-laplacien.jpg', '-djpg','-S891,630'); end 

if (n==16) ; 
plot (xx,ye,'*k;solution analytique;' , xvf,yevf,'*k;;' , ... 
        xx,yy,'*r;differences finies;' , xvf , yvf , '*b;volumes finis;' ) 
 axis ( [ -0.01, 1.01 , -0.055  , 0.055 ]) ;
title ('Laplacien, Dirichlet homogene, second membre polynomial de degre 1, n=16') ;
print('01b-dd3-dirichlet-laplacien.jpg', '-djpg','-S891,630'); end 

if (n==32) ; 
plot (xx,ye,'*k;solution analytique;' , xvf,yevf,'*k;;' , ... 
        xx,yy,'*r;differences finies;' , xvf , yvf , '*b;volumes finis;' ) 
 axis ( [ -0.01, 1.01 , -0.055  , 0.055 ]) ;
title ('Laplacien, Dirichlet homogene, second membre polynomial de degre 1, n=32') ;
print('01c-dd3-dirichlet-laplacien.jpg', '-djpg','-S891,630'); end 

if (n==64) ; 
plot (xx,ye,'*k;solution analytique;' , xvf,yevf,'*k;;' , ... 
        xx,yy,'*r;differences finies;' , xvf , yvf , '*b;volumes finis;' ) 
 axis ( [ -0.01, 1.01 , -0.055  , 0.055 ]) ;
title ('Laplacien, Dirichlet homogene, second membre polynomial de degre 1, n=64') ;
print('01d-dd3-dirichlet-laplacien.jpg', '-djpg','-S891,630'); end 

  pause (.5) 

nnn(lll) = log(n) / log(2)  ; 
erreur (lll,1) = log(errli) / log(2) ; 
erreur (lll,2) = log(errl2) / log(2) ; 
erreur (lll,3) = log(evfli) / log(2) ; 
erreur (lll,4) = log(evfl2) / log(2) ; 

n = 2*n ; 
end


hold off
%%%  plot (nnn,erreur(:,1),'*g;inf-df;' , nnn,erreur(:,2),'*m;L2-df;', ... 
%%%       nnn,erreur(:,3),'xr;inf-vf;' , nnn,erreur(:,4),'xb;L2-df;' )


  plot (nnn,erreur(:,3),'*k;volumes finis  norme L infini;', ... 
        nnn,erreur(:,4),'ok;volumes finis  norme L2;', ... 
        nnn,erreur(:,1),'*r;differences finies  norme L infini;', ... 
        nnn,erreur(:,2),'or;differences finies  norme L2;' )


  
axis ( [ 2.8 , 8.2 , -59 , -5 ]) ;
title('convergence differences finies et volumes finis')
%% print -dpdf dd3-erreurs.pdf
print('01e-dd3-erreurs.jpg', '-djpg','-S891,630');

hold off 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   

%%%  'ordres de convergence' 
%%%   [ toto1  , toto2 ]  = polyfit ( nnn , erreur(:,1) , 1) ; odfli=-toto1(1)  
%%%   [ toto1  , toto2 ]  = polyfit ( nnn , erreur(:,2) , 1) ; odfl2=-toto1(1)  
%%%   [ toto1  , toto2 ]  = polyfit ( nnn , erreur(:,3) , 1) ; ovfli=-toto1(1)  
%%%   [ toto1  , toto2 ]  = polyfit ( nnn , erreur(:,4) , 1) ; ovfl2=-toto1(1)  

fprintf ('ordres de convergence \n') 
 
fprintf ('differences finies, normes L infini et L2 \n') 
odfli = - codre(nbnn,nnn, erreur(:,1))
odfl2 = - codre(nbnn,nnn, erreur(:,2))
fprintf ('pour ce cas test, la methode des differences finies permet de calculer l interpole de la solution exacte \n') 
fprintf ('l erreur est simplement liee a la precision de la machine  \n') 
fprintf ('la notion d ordre de convergence n est pas significative \n \n') 

 
fprintf ('volumes finis, normes L infini et L2 \n') 
ovfli = - codre(nbnn,nnn, erreur(:,3))
ovfl2 = - codre(nbnn,nnn, erreur(:,4))
fprintf ('la methode des volumes finis est d ordre 2 \n') 
