%
%  quadratic penalty fn for an equality constrained problem
%
%  Nick Gould

%mu = 100.0;
%mu = 1.0;
%mu = 0.1;
mu = 0.01;
%mu = 0.0001;
inc = .005;
%inc = .1;
rho = 0.5 / mu ;
x = 0:inc:1;
y = 0:inc:1;
side = length( x );
[X,Y]  = meshgrid(x,y);
%ll0 = 0:0.05:5;
ll0 = 0:0.1:1;
%ll1 = 5:0.1:10;
ll1 = 1:0.5:3;
ll2 = 3:1:50;
ll3 = 50:10.:200;
levels = [ll0 ll1 ll2 ll3];
F = X.^2 + Y.^2 ;
C1 = X + Y.^2 - 1 ;
Q = C1.^2 ;
P = F + rho * Q ;
contour(x,y,P,levels);
%contour(x,y,P);
axis( [ 0 1 0 1 ] )
axis( 'square')

hold on

contour(x,y,C1,[0 0],'r');

xsol = 0.5 ;
ysol = sqrt( xsol ) ;

u = xsol:inc:xsol ;
v = ysol:inc:ysol ;

plot(u, v, 'ok', 'MarkerFaceColor','g', 'MarkerSize', 5 )

hold off

grid off

print -depsc eps/qpfp01

