Contents
syms x y
(part-a) The original silly objective of NW-2.2
f1 = ( 8*x + 12*y + x.^2 - 2*y^2 )
f1_grad = [ diff(f1,x); diff(f1,y) ]
f1_hess = [ diff(diff(f1,x),x) diff(diff(f1,x),y); ...
diff(diff(f1,y),x) diff(diff(f1,y),y) ]
sp_1 = eval([ solve(f1_grad(1)); solve(f1_grad(2)) ])
f1_i = inline(char(f1))
f1_grad_i = ...
inline( [ '[' char(f1_grad(1)) ';' char(f1_grad(2)) ']' ] )
f1_hess_i = ...
inline( [ '[[' char(f1_hess(1,1)) ' ' char(f1_hess(1,2)) '];' ...
'[' char(f1_hess(2,1)) ' ' char(f1_hess(2,2)) ']]' ], ...
'x','y')
hess1_sp1 = f1_hess_i(sp_1(1),sp_1(2))
xv=(-1:0.05:1);
yv=(-1:0.05:1);
[X,Y] = ndgrid(xv,yv);
Z = zeros(size(X));
for i=1:size(X,1)
for j=1:size(X,2)
Z(i,j) = ...
f1_i(sp_1(1),sp_1(2)) + ...
1/2 * [xv(i) yv(j)] * f1_hess_i(sp_1(1),sp_1(2)) * [xv(i);yv(j)];
end
end
figure(101)
contourf(X+sp_1(1),Y+sp_1(2),Z,25)
axis('square')
hold on;
plot(sp_1(1),+sp_1(2),'ro')
hold off
colorbar
xa1 = sp_1(1) + 0.25;
ya1 = sp_1(2) + 0.5;
ga1 = f1_grad_i(xa1,ya1)
ha1 = f1_hess_i(xa1,ya1)
sd1 = -ga1/norm(ga1);
nd1 = -ha1\ga1;
hold on
plot([xa1 xa1+sd1(1)],[ya1 ya1+sd1(2)],'ko-')
plot([xa1 xa1+nd1(1)],[ya1 ya1+nd1(2)],'bo-')
hold off
legend('Model','Stat. Pt.','SD-direction','Newton-direction',...
'Location','NorthEastOutside')
title(['f(x,y) = ' char(f1)])
f1 =
8*x+12*y+x^2-2*y^2
f1_grad =
8+2*x
12-4*y
f1_hess =
[ 2, 0]
[ 0, -4]
sp_1 =
-4
3
f1_i =
Inline function:
f1_i(x,y) = 8*x+12*y+x^2-2*y^2
f1_grad_i =
Inline function:
f1_grad_i(x,y) = [8+2*x;12-4*y]
f1_hess_i =
Inline function:
f1_hess_i(x,y) = [[2 0];[0 -4]]
hess1_sp1 =
2 0
0 -4
ga1 =
0.5000
-2.0000
ha1 =
2 0
0 -4
(part-b) Modified objective
f2 = ( 8*x + 12*y + x.^2 - 2*y.^2 + (1/6)*y.^4 )
f2_grad = [ diff(f2,x); diff(f2,y) ]
f2_hess = [ diff(diff(f2,x),x) diff(diff(f2,x),y); ...
diff(diff(f2,y),x) diff(diff(f2,y),y) ]
sp_2_x = eval([ solve(f2_grad(1)) ])
sp_2_y = eval([ solve(f2_grad(2)) ])
sp_2 = [ sp_2_x ; sp_2_y(1) ];
f2_i = inline(char(f2))
f2_grad_i = ...
inline( [ '[' char(f2_grad(1)) ';' char(f2_grad(2)) ']' ] )
f2_hess_i = ...
inline( [ '[[' char(f2_hess(1,1)) ' ' char(f2_hess(1,2)) '];' ...
'[' char(f2_hess(2,1)) ' ' char(f2_hess(2,2)) ']]' ], ...
'x','y')
hess2_sp2 = f2_hess_i(sp_2(1),sp_2(2))
xv=(-1:0.05:1);
yv=(-1:0.05:1);
[X,Y] = ndgrid(xv,yv);
Z = zeros(size(X));
for i=1:size(X,1)
for j=1:size(X,2)
Z(i,j) = ...
f1_i(sp_2(1),sp_2(2)) + ...
1/2 * [xv(i) yv(j)] * f2_hess_i(sp_2(1),sp_2(2)) * [xv(i);yv(j)];
end
end
figure(201)
contourf(X+sp_2(1),Y+sp_2(2),Z,25)
axis('square')
hold on;
plot(sp_2(1),+sp_2(2),'ro')
hold off
colorbar
xa2 = sp_2(1) + 0.25;
ya2 = sp_2(2) + 0.5;
ga2 = f2_grad_i(xa2,ya2)
ha2 = f2_hess_i(xa2,ya2)
sd2 = -ga2/norm(ga2);
nd2 = -ha2\ga2;
hold on
plot([xa2 xa2+sd2(1)],[ya2 ya2+sd2(2)],'ko-')
plot([xa2 xa2+nd2(1)],[ya2 ya2+nd2(2)],'ro-')
hold off
legend('Model','Stat. Pt.','SD-direction','Newton-direction',...
'Location','NorthEastOutside')
title(['f(x,y) = ' char(f2)])
f2 =
8*x+12*y+x^2-2*y^2+1/6*y^4
f2_grad =
8+2*x
12-4*y+2/3*y^3
f2_hess =
[ 2, 0]
[ 0, -4+2*y^2]
sp_2_x =
-4
sp_2_y =
-3.3681
-63.2485
66.6166
f2_i =
Inline function:
f2_i(x,y) = 8*x+12*y+x^2-2*y^2+1/6*y^4
f2_grad_i =
Inline function:
f2_grad_i(x,y) = [8+2*x;12-4*y+2/3*y^3]
f2_hess_i =
Inline function:
f2_hess_i(x,y) = [[2 0];[0 -4+2*y^2]]
hess2_sp2 =
2.0000 0
0 18.6885
ga2 =
0.5000
7.7435
ha2 =
2.0000 0
0 12.4522