//
////////////////////////////////////////////////////////////
//          probleme de la canette (+ exemple quadratique)
// 
//   méthode d'Uzawa 
// 
///////////////////////////////////////////////////////////////
//
//
function [Xt,mut]=Uzawa(J,gradJ,cI,gradcI,X0,mu0,N,rho1,rho2)
   XN=X0;muN=mu0;Xt=XN;mut=muN;
   for i=1:N
      d=-gradJ(XN)-muN*gradcI(XN)
      XN=XN+rho1*d
      muN=max(0,muN+rho2*cI(XN))
      Xt=[Xt,XN];mut=[mut,muN]
   end      
endfunction
//
//////////////////
// exemple simple: quad+contrainte lineaire//
//   A=[2,-1,0;-1,2,-1;0,-1,2];b=[-1;1;6]
//
//function y=J(x)
//       y=1/2*x'*A*x-b'*x
//endfunction   
//function y=gradJ(x)
//       y=A*x-b
//endfunction   
//function y=cI(x)
//       y=x(1)
//endfunction 
//function y=gradcI(x)
//       y=[1;0;0]
//endfunction 
//X0=[-3;-5;10];mu0=-2;rho1=0.1;rho2=0.1;N=150;
//[Xt,mut]=Uzawa(J,gradJ,cI,gradcI,X0,mu0,N,rho1,rho2)
////
//////////////////
//figure(1)
//clf()
// plot(0:N,Xt(1,:),'b-')
// plot(0:N,Xt(2,:),'r-')
// plot(0:N,Xt(3,:),'g-')
//
//figure(2)
//clf()
//plot(0:N,mut,'r-')
//
//////// canette/////////////
////contrainte inégalité: volume cannette >330 cm3
/////////////////////////////
V0=330;
Vrmin=(V0/2/%pi)^(1/3)
Vhmin=V0/(%pi*rmin^2)

function s=surfcan(X)   
    r=X(1,:);h=X(2,:);
    s=2*%pi*r.^2+2*%pi*r.*h
endfunction
//
function s=cI(X)   
    r=X(1,:);h=X(2,:);
    s=V0-%pi*r.^2.*h
endfunction
//
function s=gradsurfcan(X)
    r=X(1);h=X(2);
    s=[4*%pi*r+2*%pi*h;2*%pi*r]
endfunction
//
function s=gradcI(X)
    r=X(1);h=X(2);
    s=-[2*%pi*r*h;%pi*r^2]
endfunction
X0=[6;8]
mu0=1;
N=2000;
rho1=1E-3
rho2=2E-2
[Xt,mut]=Uzawa(surfcan,gradsurfcan,cI,gradcI,X0,mu0,N,rho1,rho2)
//////////////////
figure(1)
clf()
 plot(0:N,Xt(1,:),'r-')
 plot(0:N,Xt(2,:),'r-')
 plot(0:N,rmin*ones(N+1,1),'b-')
plot(0:N,hmin*ones(N+1,1),'b-')
figure(2)
clf()
plot(0:N,mut,'r-')
