%
% Steepest Descent Method
%

clear all;
close all;

tol = 1e-6; 
IterMax = 20;

A = [15 2; 2 15];  % 等高线不是很扁 
xt = [1; 1];
b = A*xt;   

% initial guess
x0 = [-0.5; 0]; 

% SD iteration 
n = length(b);
x = x0; 
xx = x; % store the iterate solution

r = b - A*x; 
r0norm = norm(r);
if abs(r0norm)<1e-16, return, end 
fprintf('k=0, x=[%.4f,%.4f], r0=%.2e\n', x(1),x(2),r0norm);
q = A * r;
for k = 1 : IterMax
    alpha = r'*r / (r'*q);
    x = x + alpha * r;
    xx = [xx,x];
    r = r - alpha*q;
    relres = norm(r) / r0norm;
    fprintf('k=%d, x=[%.8f,%.8f], relres=%.2e\n', k,x(1),x(2),relres);
    if relres < tol
        break
    end
    q = A * r;
end

figure
% f = @(x,y) 1/2*[x y] * A * [x; y] - [x y]*b;
f = @(x,y) 15/2*x.*x + 2*x.*y + 15/2*y.*y - 17*x - 17*y;
hold on 
set(gcf,'position',[50,50,900,900]);
fc = fcontour(f,[-4,4,-4,4],'b-','linewidth',1);
fc.LevelList = [-16 -12 0 15 45 80 120];
plot([-4,4],[0,0],'r-', [0,0],[-4,4],'r-')
% plot(1,1,'sb','markersize',10,'linewidth',2)
% text(-0.4,-0.2,'\it O','fontsize',20)
text([-3:-1]-0.2,[zeros(3,1)]-0.2, num2cell(-3:-1),'fontsize',20)
text([1:3]-0.1,[zeros(3,1)]-0.2, num2cell(1:3),'fontsize',20)
text([zeros(3,1)]-0.3,[-3:-1], num2cell([-3:-1]),'fontsize',20)
text([zeros(3,1)]-0.2,[1:3], num2cell([1:3]),'fontsize',20)
text(-0.2,-0.2, {0}, 'fontsize', 20)
set(gca,'fontsize',20)
grid on
axis equal
xticks([-4:4]); yticks([-4:4]);
set(gca,'YTickLabel',[],'XTickLabel',[]);

plot(xx(1,:),xx(2,:),'*--r','linewidth',1,'markersize',10)
hold off

    
