% graconju03.m
% fd, 10, 13 novembre 2006. 
% fd, 28, 29 juin 2026
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   clear , format long
   
printf("graconju03 ; gradient conjugue non preconditionne \n"); 
printf("             matrice penta-diagonale du laplacien a 2 dimensions\n"); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

n = 30 ; nbiter = n ; dx = 1 / n ; 
%                                            grille 
xx = zeros(1,n+1) ;  yy = zeros(1,n+1) ; 
for i = 1 : n+1 ;  xx(1,i:i) = (i-1) * dx ; end ;  
for j = 1 : n+1 ;  yy(1,j:j) = (j-1) * dx ; end ; 
%
uu=zeros(n+1,n+1) ;  bb=uu ; dx = 1 / n ; 
sgn = -1 ; %  sgn = 0 ; 
for i = 1 : n+1 ;  for j = 1 : n+1 ; 
bb(i:i,j:j) = 1 ; 
uu(i:i,j:j) = sgn ; sgn = -sgn ;  end ;  end ;  
for j = 1 : n+1 ; 
uu(1:1,j:j) = 0 ; uu(n+1:n+1,j:j) = 0 ;  
uu(j:j,1:1) = 0 ; uu(j:j,n+1:n+1) = 0 ;  end ; 
%

hold off 
mesh ( xx , yy , uu ) ;
title ('solution initiale de l algorithme de gradient conjugue') ;
print('03-a-gradient-conjugue-matrice-penta-diagonale', '-djpg','-S891,630'); 
mesh ( xx , yy , uu ) ;
pause (.5) 

gr = uu ;  for i = 2 : n ; for j = 2 : n ;
gr(i:i,j:j) = ( 4*uu(i:i,j:j) - uu(i:i,j-1:j-1) - uu(i:i,j+1:j+1) \
      - uu(i-1:i-1,j:j) - uu(i+1:i+1,j:j) ) / dx**2  - bb(i:i,j:j) ; 
end ; end ; 

% mesh ( xx , yy , gr  ) ; pause
 
ww = gr ;   %  preconditionnement a prevoir ici ; 


gg = 0 ; for i = 2 : n ; for j = 2 : n ;
gg = gg + gr(i:i,j:j) * gr(i:i,j:j) ; end ; end ; 

zco = zeros(1,nbiter) ;  
for iter = 1:nbiter  ; 
waw = 0 ; gw = 0 ;   
 for i = 2 : n ;  for j = 2 : n ; 
waw = waw +  ww(i:i,j:j) *  \
( 4*ww(i:i,j:j) - ww(i:i,j-1:j-1) - ww(i:i,j+1:j+1) \
      - ww(i-1:i-1,j:j) - ww(i+1:i+1,j:j) ) / dx**2 ; 
gw = gw + gr(i:i,j:j) * ww(i:i,j:j) ; 
end ; end ;  
zco(1,iter:iter) = log(gw) / (2*log(10)) ;  
ro = - gw / waw ; 
for i = 2 : n ;   for j = 2 : n ; 
uu(i:i,j:j) = uu(i:i,j:j) + ro*ww(i:i,j:j) ; 
end  ; end ; 
mesh ( xx , yy , uu  ) ;


hold off 
mesh ( xx , yy , uu ) ; zlim([0 0.08])    % set limits of axis

pause (.1) 

%                               calcul de la fonction cout ; 
cout = 0 ;  for i = 2 : n ;  for j = 2 : n ;  
cout = cout + .5 * uu(i:i,j:j) * \
( 4*uu(i:i,j:j) - uu(i:i,j-1:j-1) - uu(i:i,j+1:j+1) \
      - uu(i-1:i-1,j:j) - uu(i+1:i+1,j:j) ) / dx**2 ; 
cout = cout -  bb(i:i,j:j) * uu(i:i,j:j) ; end ; end ; 
% zco(1,iter:iter) = cout ;  cout
% 
grn = gr ;  for i = 2 : n ; for j = 2 : n ;
grn(i:i,j:j) = ( 4*uu(i:i,j:j) - uu(i:i,j-1:j-1) - uu(i:i,j+1:j+1) \
      - uu(i-1:i-1,j:j) - uu(i+1:i+1,j:j) ) / dx**2  - bb(i:i,j:j) ; 
end ; end ;  ggn = 0 ;
for i = 2 : n ;  for j = 2 : n ;  
ggn = ggn + grn(i:i,j:j) * grn(i:i,j:j) ; end ; end ; 
ggn ;  alpha = ggn / gg ; 
gr = grn ; gg = ggn ; 
  for i = 2 : n ;  for j = 2 : n ; 
ww(i:i,j:j) = gr(i:i,j:j)  + alpha *  ww(i:i,j:j) ; end ;  end ;  

% mesh ( xx , yy , ww  ) ; 
end ; 


hold off 
mesh ( xx , yy , uu ) %%% ; zlim([0 0.08])    
title ('solution finale de l algorithme de gradient conjugue') ;
print('03-b-gradient-conjugue-matrice-penta-diagonale', '-djpg','-S891,630'); 
mesh ( xx , yy , uu ) %% ; zlim([0 0.8])  




iters = zeros(1,nbiter) ; for j = 1:nbiter ; 
iters(1,j:j) = j ; 
%  zco (1,j:j) = zco (1,j:j) - zco(1,nbiter:nbiter) -1 ;
   end ; 

pause (1)
plot ( zco  , '*b;logarithme de la norme du gradient;' ) ;

title ('convergence de l algorithme') ;
print('03-c-gradient-conjugue-matrice-penta-diagonale', '-djpg','-S891,630'); 



