function maratos

% Maratos effect for the l-1 penalty function for 
%
%   min 2(x^2 + y^2 - 1) + x st   x^2 + y^2 = 1
%
% 1/Jun/01
%

clear

h_space = 0.005;
%h_space = 0.05;
v_space = 0.2;
c_space = 1.0;
xx =  0.2:h_space:1.2 ;
yy =  -0.8:h_space:0.2 ;
[X,Y]=meshgrid(xx,yy);
[f,fc1,fint] = pf( X , Y ) ;

colormap gray

% This makes the background green (don't ask why!!)

contourf(X,Y,-fc1,[0 0], 'k' );

hold on

contourf(X,Y,fc1,[0 0], 'k' );

colormap summer

nz = -15:v_space:30 ;
nz2 = -15:c_space:2 ;

hold on

contour(X,Y,fint,nz,'y');

%contour(X,Y,-fc1,[0 0], 'g' );
contour(X,Y,-fc1,[0 0], 'b' );

% end

pi = 4 * atan(1.0) ; 
theta = - 0.2 * pi ;
%theta = - 0.05 * pi ;
ct = cos(theta) ;
st = sin(theta) ;
x1 = ct ;
y1 = st ;

sx1 = st*st ;
sy1 = - ct*st ;
x1p = x1 + sx1 ;
y1p = y1 + sy1 ;

x = [x1,x1p] ;
y = [y1,y1p] ;

plot(x,y,'ro-');

x1pm = x1 + 0.5 * sx1 ;
y1pm = y1 + 0.5 * sy1 ;
x = [x1,x1pm] ;
y = [y1,y1pm] ;

solx = [1] ;
soly = [0] ;

plot(solx,soly,'bo');


hold off

print -depsc eps/maratos

% penalty function

function [f,fc1,fint] = pf( x, y )
rho = 1.0 ;
f = obj( x, y ) ;
fc1= c1( x, y ) ;
fint = f+rho*abs(fc1);

function o = obj( x, y )
o = 2 * ( x.^2 + y.^2 - 1) - x;

function c = c1( x, y )
c = x.^2 + y.^2 - 1;

small = 0.01 ;
side = 2 * small ;

