function [X,lamb]=uzawa(J, nablaJ,cE,nablacE,X0,lamb0,alpha1,alpha2,N)
    X=X0;lamb=lamb0;
    for i=1:N
        X=X-alpha1*(nablaJ(X)+lamb*nablacE(X))
        lamb=lamb+alpha2*cE(X)
    end
endfunction
//////////////////
// exemple 1
//////////////////
A=[2,-1,0;-1,2,-1;0,-1,2];b=[-1;1;6];
function y=quad(x)
    y=1/2*(A*x)'*x-b'*x
endfunction
function y=nablaquad(x)
    y=A*x-b
endfunction
function y=cEex(x)
    y=x(1)
endfunction        
function y=nablacEex(x)
    y=[1;0;0]
endfunction  
////////////////// 
alpha1=0.01;
alpha2=0.01;
X0=[1;1;1];
lamb0=1;
N=10000;
[X,lamb]=uzawa(quad,nablaquad,cEex,nablacEex,X0,lamb0,alpha1,alpha2,N)
disp('solution exemple 1 ');
disp(X)
//////////////////
// exemple 2
//////////////////
function y=rosen(x)
    y=100*(x(2)-x(1)^2)^2+(x(1)-1)^2
endfunction
//
function y=rosengrad(x)
y=[-400*x(1)*(x(2)-x(1)^2)+2*(x(1)-1);200*(x(2)-x(1)^2)]
endfunction
//
function y=cEex2(x)
    y=2*x(1)+x(2)-4
endfunction        
function y=nablacEex2(x)
    y=[2;1]
endfunction  
//////////////////////
alpha1=0.001;
alpha2=0.001;
X0=rand(2,1);
lamb0=0;
N=100000;
[X,lamb]=uzawa(rosen,rosengrad,cEex2,nablacEex2,X0,lamb0,alpha1,alpha2,N)
disp('solution exemple 2');
disp(X)
