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

   clear, format long 


fprintf ('dn4 \n') ;
fprintf ('probleme de Dirichlet-Neumann pour le Laplacien \n') ;
fprintf ('differences finies, elements finis et volumes finis  \n') 
fprintf ('solution exacte = x (1-x) (1/3-x) (2/3 - x) \n') 
fprintf ('octave 6.4 \n')



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

n = 8   ;  

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  ;  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   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) = 22/9 - 12*xa + 12*xa**2 ; xa = xa + dx ;    end   
bdf(n) = -2/9 ;   
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
xm = 0 ; xa = dx/2 ;  xp = dx ;     
for  j = 1:n ;       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    point milieu
bvf(j) =  22/9 - 12*xa + 12*xa**2  ;  xa = xa + dx ;   end  
% bvf(j) =  22/9 - 12*(xm/6 + (2*xa/3) + xp/6) + \
%              12* ( (xm**2)/6 + (2*(xa**2)/3) + (xp**2)/6 )   ;  
% xm = xm + dx ;  xa = xa + dx ; xp = xp + dx ;   end  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    Simpson 
bvf(n) = bvf(n) - (2/9)/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) =  22/9 - 12*xa + 12*xa**2 ; 
xa = xa + dx ; end 

bef = mef*ff ; 
bef(n) = bef(n) -2/9 ; 


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

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%   erreurs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    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/3-xx(j))*(2/3-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 ; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%       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/3-xa) * (2/3-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 ;   


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%      graphiques    

if (n==8) ; 
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.011  , 0.018 ]) ;
title ('Laplacien, Dirichlet Neumann, second membre polynomial de degre 2, n=8') ;
print('04a-dn4-dirichlet-neumann-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;' , ...
      xx,yef,'*m;elements finis;' ) 
 axis ( [ -0.01, 1.01 , -0.011  , 0.018 ]) ;
title ('Laplacien, Dirichlet Neumann, second membre polynomial de degre 2, n=16') ;
print('04b-dn4-dirichlet-neumann-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;' , ...
      xx,yef,'*m;elements finis;' ) 
 axis ( [ -0.01, 1.01 , -0.011  , 0.018 ]) ;
title ('Laplacien, Dirichlet Neumann, second membre polynomial de degre 2, n=32') ;
print('04c-dn4-dirichlet-neumann-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;' , ...
      xx,yef,'*m;elements finis;' ) 
 axis ( [ -0.01, 1.01 , -0.011  , 0.018 ]) ;
title ('Laplacien, Dirichlet Neumann, second membre polynomial de degre 2, n=64') ;
print('04d-dn4-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 -.02 .02])
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) ; 

n = 2*n ; 
end 

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;'  )
axis ( [ 2.8 , 8.2 , -19 , -6.5 ]) ;
  title('convergence differences finies, volumes finis, elements')
print('04e-dn4-dirichlet-neumann-laplacien.jpg', '-djpg','-S891,630');
%% print('04-dn4-erreurs.jpg', '-djpg','-S891,630');


  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%  ordre de convergence    
 
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 2 \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 ('la methode des elements finis est d ordre 2 \n') 
 


 

