% graconju02.m
% fd, 10, 13 novembre 2006.
% fd, 28 juin 2026
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   clear , format long
   
printf("graconju02 ; gradient conjugue preconditionne \n"); 
printf("             cas de la matrice tridiagonale du laplacien\n"); 
printf("             preconditionnement diagonal\n"); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

n = 100 ;   %  la matrice 
aa = zeros(1,n-1) ; bb = aa ; cc = aa ; 
%  donnes du laplacien 
for j=1:n-1 ; 
aa(1,j:j) = 2 ;  bb(1,j:j) = -1 ;   cc(1,j:j) = -1 ; end ; 
sm = zeros(1,n-1) ; %  le second membre 
for j=1:n ; sm(1,j:j) = 1 ; end ; 
xx = zeros(1,n-1) ; %  la solution  
xx=tridia(n-1,aa,bb,cc,sm) ; 
%                               algo du gradient conjugue
uu=zeros(1,n+1) ;  bb=uu ; dx = 1 / n ; 
sgn = -1 ; 
for j = 1 : n+1 ; bb(1,j:j) = 1 ; 
uu(1,j:j) = sgn ; sgn = -sgn ;  end ;  
uu(1,1:1) = 0 ;  uu(1,n+1:n+1) = 0 ; 
gr = uu ; 
for j = 2 : n ;
gr(1,j:j) = ( 2*uu(1,j:j) - uu(1,j-1:j-1) - uu(1,j+1:j+1)) / dx**2 \
            - bb(1,j:j) ; end ; 
ww = gr ;   for j = 2 : n ; ww(1,j:j)  = .5 * gr(1,j:j)  ;  end ; 

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

% for iter = 1:n/2  ; 
nmax = 20 
nmax = 30 
nmax = 40 
nmax = 50 
for iter = 1:nmax ; 
waw = 0 ; gw = 0 ;   
for j = 2 : n ; 
waw = waw + \ 
ww(1,j:j) * ( 2*ww(1,j:j) - ww(1,j-1:j-1) - ww(1,j+1:j+1)) / dx**2 ; 
gw = gw + gr(1,j:j) * ww(1,j:j) ; 
end ;  % waw , gw , gg   , 
ro = - gw / waw ; 
for j = 2 : n ; 
uu(1,j:j) = uu(1,j:j) + ro*ww(1,j:j) ; 
end  ;    
%                               calcul de la fonction cout ; 
cout = 0 ;  for j = 2 : n ;  cout = cout + .5 * uu(1,j:j) * \
   ( 2*uu(1,j:j) - uu(1,j-1:j-1) - uu(1,j+1:j+1)) / dx**2  ; 
cout = cout -  bb(1,j:j) * uu(1,j:j) ; end ; cout 
% 
grn = gr ;  for j = 2 : n ;
grn(1,j:j) = ( 2*uu(1,j:j) - uu(1,j-1:j-1) - uu(1,j+1:j+1)) / dx**2 \
            - bb(1,j:j) ; end ;  ggn = 0 ;
for j = 2 : n ;  ggn = ggn + .5*grn(1,j:j) * grn(1,j:j) ; end ; 
ggn , alpha = ggn / gg 
gr = grn ; gg = ggn ; 
for j = 2 : n ; 
ww(1,j:j) = .5*gr(1,j:j)  + alpha *  ww(1,j:j) ; end ;  
%
%                               formattage des donnees 
%
xxx = zeros(1,n+1) ; uuu = xxx ; dx = 1 / n ; 
for j = 1 : n+1 ; xxx(1,j:j) = (j-1) * dx ; end ; 
for j = 2 : n ; uuu(1,j:j) = xx(1,j-1:j-1) * dx**2 ; end ;  
solex = uuu ;  
solgr = uu ; solgr(1,1:1) = 0 ;  solgr(1,n+1:n+1) = 0 ; 
for j = 2 : n ; x =  xxx(1,j:j) ;  solex(1,j:j) = x*(1-x) / 2 ; end ; 


plot (xxx,uuu,'*;;',xxx,uuu, 'b;;' , xxx, solex,  'g;;' , 
 xxx, solgr,  'b;;'  ) ; pause (.1) ; 
end ; 

hold on ;
plot(xxx,uuu,'*;solution exacte;', xxx, solgr, "b;algorithme de gradient;") ;

if nmax==20 ; title ('etat de l algorithme apres 20 iterations de gradient') ;
print('02-a-gradient-conjugue-preconditionne-diagonale', '-djpg','-S891,630');  endif

if nmax==30 ; title ('etat de l algorithme apres 30 iterations de gradient') ;
print('02-b-gradient-conjugue-preconditionne-diagonale', '-djpg','-S891,630');  endif

if nmax==40 ; title ('etat de l algorithme apres 40 iterations de gradient') ;
print('02-c-gradient-conjugue-preconditionne-diagonale', '-djpg','-S891,630');  endif

if nmax==50 ; title ('etat de l algorithme apres 50 iterations de gradient') ;
print('02-d-gradient-conjugue-preconditionne-diagonale', '-djpg','-S891,630');  endif
