% laplacien04.m
% fd, 04 octobre 2009. 
% fd, 01 juillet 2026
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear , format long
   
printf("laplacien02 ; probleme de Neumann par differences finies  \n"); 
printf("              iterations de la methode du gradient  \n"); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

printf("\n");

nex= 70  , n= 30 , rho = .01 
xex  = [0:1/nex:1] ; 
yex = xex ; 
for j = 1:nex+1 ; 
yex(j)  = 1 + 12*xex(j)**2-8*xex(j)**3; 
end 

hold off
%%%%% axis ( [-.02, 1.02 , .9  , 5.7 ] ) ; 
%%%%% plot ( xex ,  yex , 'b;;') 
%%%%% plot ( xex ,  yex , 'b;solution exacte;') 
 
%          maillage  
h = 1/n ;
x  = [0:h:1] ; y = [h/2:h:1-h/2] ; 

%          second membre
f = y ; 
for j = 1:n ; 
f(j)  = -23 + 48*y(j) + 12*y(j)**2 - 8*y(j)**3 ; 
end 
% plot ( y ,  f , '*r;second membre de l edp;') 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%         methode du gradient simple 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                initialisation 
u = x ; 
for j = 1:n+1 ; 
%%  u(j)  = 1 + 4*x(j) ;
u(j)  = 1 ;
end 

%%% plot ( x ,  u , 'g;;')
  
plot ( xex ,  yex , 'b;;' ,  x ,  u , '*k;;')
%%% plot ( xex ,  yex , '*b;;' , xex ,  yex , 'b;;')
title ('initialisation par la fonction constante egale a 1') ;
axis ( [ -.025, 1.025, .9, 5.1 ]) 
print('04a-probleme-de-neumann-pour-un-laplacien-1d', '-djpg','-S891,630');  
hold on 

for lll = 1 : 10 
for ll = 1 : 2000 

%          calcul du gradient de la fonctionnelle 
gr = x ; 
gr(1) = (u(1)-u(2))/h + h*(u(1)+u(2))/4 - h*f(1) / 2 ; 
for j = 2:n ; 
gr(j) = (-u(j-1) +2*u(j)-u(j+1))/h + h*(u(j-1) +2*u(j)+u(j+1))/4 \
               - h*(f(j-1)+f(j))/2 ;  
end 
gr(n+1) = (u(n+1)-u(n))/h + h*(u(n)+u(n+1))/4 - h*f(n) / 2 ; 

% plot ( x ,  gr , '*g;gradient;') 

%          incrementaion de u  
for j = 1:n+1 ; 
u(j) = u(j) - rho*gr(j) ; 
end 
% plot ( x ,  u , 'r;;') 
end

if (lll==1) 
plot ( x ,  u , 'b*;;') 
title ('apres 2000 iterations de la methode du gradient') ;
axis ( [ -.025, 1.025, .9, 5.1 ]) 
print('04b-probleme-de-neumann-pour-un-laplacien-1d', '-djpg','-S891,630');  
pause (.5) ; hold on ; endif 

if (lll==2) 
plot ( x ,  u , 'm*;;') 
title ('apres 4000 iterations de la methode du gradient') ;
print('04c-probleme-de-neumann-pour-un-laplacien-1d', '-djpg','-S891,630');  
pause (.5) ; hold on ; endif 


if (lll==3) 
plot ( x ,  u , 'g*;;') 
title ('apres 6000 iterations de la methode du gradient') ;
print('04d-probleme-de-neumann-pour-un-laplacien-1d', '-djpg','-S891,630');  
pause (.5) ; hold on ; endif 

if (lll==4) 
plot ( x ,  u , 'r*;;') 
title ('apres 8000 iterations de la methode du gradient') ;
print('04e-probleme-de-neumann-pour-un-laplacien-1d', '-djpg','-S891,630');  
pause (.5) ; hold on ; endif 

end

plot ( x ,  u , 'k*;;') 
title ('apres 20000 iterations de la methode du gradient...') ;
print('04f-probleme-de-neumann-pour-un-laplacien-1d', '-djpg','-S891,630');  

hold off 








%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%     resolution des equations d'Euler 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



aa = zeros(1,n+1) ; bb = aa ; cc = aa ; sm = aa ; 

%                      second membre 
sm(1)   = h*f(1) / 2 ; 
for   j = 2:n ; 
sm(j)   =  h*(f(j-1)+f(j))/2 ;  
end 
sm(n+1) =  h*f(n) / 2 ;

%                      systeme lineaire  
aa(1)   = 1/h + h/4 ; 
bb(2)   = -1/h + h/4 ; 
for   j = 2:n ;  
aa(j)   =  2/h + h/2 ; 
cc(j)   = -1/h + h/4 ; 
bb(j+1) = -1/h + h/4 ; 
end 
aa(n+1) = 1/h + h/4 ; 
cc(n+1) = -1/h + h/4 ; 
 

 xx=tridia(n+1,aa,bb,cc,sm) ; 

%%%%  plot ( x ,  xx , 'r*;solution discrete;') 



