% uza7.m
% fd, 28, 30  janvier 2006, 08 juillet 2026 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear ; format long 

fprintf  ('\n') 
fprintf  ('uza7 \n')  
fprintf  ('opmisation avec deux contraintes inegalites \n') 
fprintf  ('algorithme de Uzawa a pas constant  \n  ') ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a=2 ; b=0 ; al0=2; be0=1 ;
%  axis ( [-2.1 , 2.1 , -2.1, 2.1],'square' );
   axis ( [ a-4.1 , a+4.1 , b-2.1, b+2.1] );

n=300 ; 
dth=2*pi/n ; 
xx=zeros(n+1,1); yy= zeros(n+1,1) ; 

th=0;  
hold off ; for j=1:n+1 ; 
y = 2 * sin(th) / sqrt(3) ; 
xx(j:j,1) = y/2 + cos(th) ; 
yy(j:j,1) = y ;  th=th+dth; end;
plot (xx,yy,'b;;') ; 
hold on

th=0;  for j=1:n+1 ; 
y = 2 * sin(th) / sqrt(3) ; 
xx(j:j,1) = -y/2 + cos(th) ; 
yy(j:j,1) = y ;  th=th+dth; end;
plot (xx,yy,'b;;')

axis ( [ -1.6 , 2.6 , -2.1, 2.1], 'square' );
title ('convexe des contraintes : une intersection d ellipses') ;
print('07a-uzawa-cas-de-deux-inegalites.jpg', '-djpg','-S891,630'); 
pause (.5) 
hold on


dr=.1; rr=.1 ; 
for j=1:20;
al=rr*al0 ; be=rr*be0 ; 
dth=ellips (a,b,al,be);
rr=rr+dr;
end ; 
title ('on rajoute les isovaleurs de la fonction cout') ;
print('07b-uzawa-cas-de-deux-inegalites.jpg', '-djpg','-S891,630');
pause (.5) 
hold on 

xx = -1.5  ; yy = 1 ; 
plot (xx,yy ,'r*;;')
title ('point initial') ;
print('07c-uzawa-cas-de-deux-inegalites.jpg', '-djpg','-S891,630');
pause (.5) 
hold on 


A = [ 2*xx+yy , 2*xx-yy ; xx+2*yy , -xx+2*yy ] ; 
bb= [ -(xx-a)/al0**2 ; -(yy-b)/be0**2 ] ;
xi = A \ bb ;  lambda = xi(1:1,1)  ;   mu = xi(2:2,1)  ;   

nituzawa = 20 
     for j = 1 :nituzawa %% 55 
A = [ 1/al0**2 + 2*(lambda+mu) ,  lambda-mu ;  lambda-mu , 1/be0**2 + 2*(lambda+mu) ] ; 
bb= [ a/al0**2 ; b/be0**2 ] ;
xi = A \ bb ;  xx1 = xi(1:1,1)   ; yy1 = xi(2:2,1) ;  
plot (xx1,yy1 ,'b*;;', [xx xx1], [yy yy1], 'r;;') ;;
pause (.2) 

xx=xx1 ; yy=yy1 ; 

dgx = xx**2 + xx*yy + yy**2 - 1 ; 
dgy = xx**2 - xx*yy + yy**2 - 1  ; 
rho = .1 ; 
lprov = lambda + rho * dgx ;  
mprov     = mu     + rho * dgy  ;  
[lambda , mu] = proj2(lprov , mprov ) ; 
end 
axis ( [ -1.6 , 2.6 , -2.1, 2.1], 'square' );
title ('iterations de l algorithme de Uzawa') ;
print('07d-uzawa-cas-de-deux-inegalites.jpg', '-djpg','-S891,630');
pause (.5)


axis ( [ -1.6 , 2.6 , -2.1, 2.1], 'square' );
title ('iterations de l algorithme de Uzawa') ;
print('07d-uzawa-cas-de-deux-inegalites.jpg', '-djpg','-S891,630');
pause (.5)

axis ( [ 0., 2., -1. , 1. ] , 'square'  );
title ('20 iterations de l algorithme de Uzawa') ;
print('07e-uzawa-cas-de-deux-inegalites.jpg', '-djpg','-S891,630');
pause (.5)


axis ( [ .5, 1.5, -.5 , .5 ] , 'square'  );
title ('20 iterations de l algorithme de Uzawa') ;
print('07f-uzawa-cas-de-deux-inegalites.jpg', '-djpg','-S891,630');
pause (.5)

eps=.1 ;  
for j = 1 : 2
axis ( [ 1 - eps  , 1 + eps  ,   - eps  ,  + eps ] , 'square' );
if j==1 
  title ('20 iterations de l algorithme de Uzawa : premier zoom') ;
print('07g-uzawa-cas-de-deux-inegalites.jpg', '-djpg','-S891,630'); end 
if j==2
  title ('20 iterations de l algorithme de Uzawa : second zoom') ;
print('07h-uzawa-cas-de-deux-inegalites.jpg', '-djpg','-S891,630'); end 

pause (.5)
eps = eps / 10 ; 
end ; 


hold off

