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

printf("\n"); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nex= 30  , 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
plot ( xex ,  yex , '*b;;' , xex ,  yex , 'b;;')
title ('solution exacte u(x) = 1 + 12 x**2 - 8*x**3') ;
axis ( [ -.025, 1.025, .9, 5.1 ]) 
print('02a-probleme-de-neumann-pour-un-laplacien-1d', '-djpg','-S891,630');  
pause (.5) ; hold off 


 
%          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;') 

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

hold off 
plot ( xex ,  yex , 'b;;') ; hold on 
plot ( x ,  u , 'r*;;')
title ('champ initial pour les iterations de la methode du gradient') ;
axis ( [ -.025, 1.025, .9, 5.1 ]) 
print('02b-probleme-de-neumann-pour-un-laplacien-1d', '-djpg','-S891,630');  
pause (.5) ; hold on


 
for lll = 1 : 3

for ll = 1 : 200 
%          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 , 'm*;;') 
title ('apres 200 iterations de la methode du gradient') ;
print('02c-probleme-de-neumann-pour-un-laplacien-1d', '-djpg','-S891,630');  
pause (.5) ; hold on ; endif 

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

end

plot ( x ,  u , 'k*;;') 
title ('solution approchee apres 600 iterations de la methode du gradient') ;
print('02e-probleme-de-neumann-pour-un-laplacien-1d', '-djpg','-S891,630');  

%%% plot ( x ,  u , 'r;solution approchee finale;') 
