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

   clear, format long 

fprintf ('dn3 \n') ;
fprintf ('probleme de Dirichlet-Neumann pour le Laplacien \n') ;
fprintf ('differences finies, volumes finis, elements finis  \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,6) ; 

n = 10   ;  

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        differences finies  

adf = zeros(n,n) ;          %  la matrice differences finies  
for j=1:n-1 ; 
adf(j,j) = 2 ;  adf(j,j+1) = -1 ;  end 
for j=1:n-1 ;   adf(j,j+1) = -1 ;  end
for j=2:n-1 ;   adf(j,j-1) = -1 ;  end
adf(n,n)  = 1 ;  adf(n,n-1) = -1 ;   %%%% schema simple a deux points 
                                       %%%% schema de Gear a trois points 
% adf(n,n)  = 3  ; adf(n,n-1) = -2 ;  adf(n,n-2) = 1/2  ;  


% adf(n,n-2) = 1/2; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   scaling 
for i = 1:n-1 , for j=1:n 
adf(i,j) =  adf(i,j) / (dx**2) ; end , end 
for j=1:n ; adf(n,j) =  adf(n,j) / dx ; end 


bdf = zeros(n,1) ;          %  second membre differences finies  
xa     = dx ;  
for  j = 1:n-1 ; 
bdf(j) = 3-6*xa ; xa = xa + dx ;    end
bdf(n) = 1/2 ;  
bdf'  ; 

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

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(n,n) = 1 ;  avf (n,n-1)  = -1 ;       
avf(1,1) = 3 ; avf (1,2) = -1 ;            %%%%%%%%%  premiere version 
%  avf(1,1) = 4 ; avf (1,2) = -4/3 ;       %%%%%%%%%  seconde version
avf  ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   scaling 
for i = 1:n , for j=1:n 
avf(i,j) =  avf(i,j) / (dx**2) ; end , end 


bvf = zeros(n,1) ;         %  second membre volumes  finis
xa = dx/2 ;     
for  j = 1:n ;    
bvf(j) =   3-6*xa ;         %%%%%%%%%%%%%%    point milieu
xa = xa + dx ;   end 
bvf(n) = bvf(n) + .5/dx ;  bvf'   ; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        elements finis  

aef = zeros(n,n) ;          %  la matrice elements finis  

for j=1:n ;   aef(j,j)   = 2 ;   end 
for j=1:n-1 ; aef(j,j+1) = -1 ;  end
for j=2:n ;   aef(j,j-1) = -1 ;  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   scaling 
for i = 1:n , for j=1:n 
aef(i,j) =  aef(i,j) /  dx  ; end , end 

mef = zeros(n,n+1) ;          %  matrice de masse  

for j=1:n ;   mef(j,j+1) = 2/3  ;   end 
for j=1:n-1 ; mef(j,j+2) = 1/6  ;   end
for j=1:n ;   mef(j,j)   = 1/6  ;   end
mef(n,n+1) = 1/3  ;   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   scaling 
for i = 1:n , for j=1:n+1 
mef(i,j) =  dx * mef(i,j)   ; end , end 

aef ; mef ; 

bef   = zeros(n,1) ;          %%%%%%%  second membre elements finis  
ff    =  zeros(n+1,1) ;  xa = 0 ;  
for j=1:n+1 ;  ff(j) = 3 - 6*xa ; xa = xa+dx ; end 

bef = mef*ff ; 
bef(n) = bef(n) + 1/2 ; 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%   algebre lineaire 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

udf = adf \ bdf ; 
uvf = avf \ bvf ; 
uef = aef \ bef ; 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%   graphique    %%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    solution exacte aux sommets  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    differences finies 
xx = zeros(1,n+1) ; yy = xx ;  ye = xx ; 
for j=1:n+1 ; xx(j) = (j-1)*dx ; end  

for j = 1:n+1 ; 
ye(j) = xx(j)*(1-xx(j))*(1/2-xx(j)) ; end 
for j = 2:n+1
yy(j) = udf(j-1) ; end 

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 ; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    solution exacte aux centres des mailles   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       volumes  finis
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  

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  ;  
 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        elements finis 
yef    =  zeros(1,n+1) ;     
for  j = 2:n+1 ; yef(j) = uef(j-1) ; end 


erefi = 0 ;     eref2  = 0 ; 
for j = 1:n
tt    = abs(yef(j) - ye(j)) ; 
erefi = max (erefi , tt ) ; 
eref2 = eref2  + tt**2 ; end 
eref2 = sqrt ( dx * eref2 ) ; erefi ;

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

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

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

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

if (n==160) ; 
plot (xx,ye,'*k;solution analytique;' , xvf,yevf,'*k;;' , ... 
      xx,yy,'*r;differences finies;' , xvf , yvf , '*b;volumes finis;' , ...
      xx,yef,'*m;elements finis;' ) 
 axis ( [ -0.01, 1.01 , -0.055  , 0.055 ]) ;
title ('Laplacien, Dirichlet Neumann, second membre polynomial de degre 1, n=160') ;
print('03e-dn3-dirichlet-neumann-laplacien.jpg', '-djpg','-S891,630'); end 



%%% plot (xx,yy,'*r;df;' ,  xvf , yvf , '*b;vf;' ,  xx,yef,'*m;ef;' , ...
%%%         xx,ye,'og;exact;' , xvf,yevf,'og;;' )
%%% axis ([0 1 -.06 .06])
pause(.5) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%       sorties synthetique des erreurs 


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) ; 
erreur (lll,5) = log(erefi) / log(2) ;
erreur (lll,6) = log(eref2) / log(2) ; 

erreur (lll,:) ; 

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-vf;' , ... 
%%%         nnn,erreur(:,5),'*c;inf-ef;' , nnn,erreur(:,6),'+r;L2-ef;'  ) 
%%%  title('convergence differences finies, volumes finis, elements finis')
%%%  print('03-dn3-erreurs.jpg', '-djpg','-S891,630');
%%% print -dpdf dn3-erreurs.pdf


  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;' , ... 
        nnn,erreur(:,5),'*b;elements finis  norme L infini;', ... 
        nnn,erreur(:,6),'ob;elements finis  norme L2;' )  

axis ( [ 3.2 , 8.8 ,  -57 , 10 ]) ;
  title('convergence differences finies, volumes finis, elements finis')
print('03f-dn3-erreurs.jpg', '-djpg','-S891,630');

hold off 




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%  ordre de convergence    

%%%%  '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)  
%%%%   [ toto1  , toto2 ]  = polyfit ( nnn , erreur(:,5) , 1) ; oefli=-toto1(1)  
%%%%   [ toto1  , toto2 ]  = polyfit ( nnn , erreur(:,6) , 1) ; oefl2=-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 ('la methode des differences finies est d ordre 1 \n') 
fprintf ('la mise en oeuvre de la condition de Neumann a droite n est pas assez precise \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') 


fprintf ('elements finis, normes L infini et L2 \n') 	  
oefli = - codre(nbnn,nnn, erreur(:,5))
oefl2 = - codre(nbnn,nnn, erreur(:,6))
fprintf ('pour ce cas test, la methode des elements finis 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') 

 

