%-----------------------------------------------------------
%     Geometricke znazorneni metody nejvetsiho spadu
%-----------------------------------------------------------

clc;
tisk=sprintf(['Geometricke znazorneni metody nejvetsiho spadu']);
disp(tisk);
disp(' ');
disp(' ');

a11=3;
a12=0;
a22=200;
b1=-8;
b2=2;

A=[a11, a12; a12, a22];
b=[b1;b2];

x0=27;
y0=0.6;

iteraci=10;
epsilon=0.001;

presne=A\b;

figure(1);
clf;
hold on;
xmin=-50;
xmax=50;
ymin=-9;
ymax=10;
xlim=([xmin xmax]);
ylim=([ymin ymax]);
title('Metoda nejvetsiho spadu');
axis([xmin xmax ymin ymax]);
set(figure(1),'Units','centimeters');
set(figure(1),'Position',[2 6 32 10]);

krok=0.05;
xx=xmin:krok:xmax;

plot([-8/3 -8/3], [-7 7],'m');
plot([-40, 40], [0.01 0.01],'b');
rce1=strcat(int2str(a11),'x+',num2str(a12,1),'y=',int2str(b1));
rce2=strcat(num2str(a12,1),'x+',int2str(a22),'y=',int2str(b2));
text(-45,-4,rce1);
text(-45,-5,rce2);

% Metoda nejvetsiho spadu
X=x0;
Y=y0;
x=[X,Y]';
iter=0; 
g1=-A*x+b;
while (iter < iteraci+1),
   x=x+g1'*g1/(g1'*A*g1)*g1;
   g1=-A*x+b;  
   X=[X,x(1)];
   Y=[Y,x(2)];
   iter=iter+1;
end

  tisk=sprintf( ... 
  ['| k |   x(k)   |   y(k)   | x(k)-x(k-1) | y(k)-y(k-1) | norma chyby |']);
  disp(tisk);
  delka=length(tisk);
  retezec='';
  if (delka>80), delka=80; end;
  iter=0;
  while iter<delka
   retezec=['-',retezec];
   iter=iter+1;
  end;
  disp(retezec);

norma=norm(presne-[X(1), Y(1)]');
tisk=sprintf('| 0 |%9.6f |%9.6f |             |             |  %9.6f  |',...
	      X(1), Y(1), norma);
disp(tisk);

for i=1:iteraci
 plot(X(i),Y(i),'ro'); 
 hladina=0.5*[X(i),Y(i)]*A*[X(i),Y(i)]'-b'*[X(i),Y(i)]';
 elipsa(xx,hladina,krok,a11,a12,a22,b1,b2);
 pause;
 legend('1. rovnice','2. rovnice','iterace','hladina')
 plot([X(i),X(i+1)],[Y(i),Y(i+1)],'k--');
 norma=norm(presne-[X(i+1), Y(i+1)]');
 tisk=sprintf('|%2d |%9.6f |%9.6f |  %9.6f  |  %9.6f  |  %9.6f  |',...
	       i,X(i+1), Y(i+1),X(i+1)-X(i),Y(i+1)-Y(i), norma);
 disp(tisk);
end;
plot(X(iteraci+1),Y(iteraci+1),'ro'); 
hladina=0.5*[X(iteraci),Y(iteraci)]*A*[X(iteraci),Y(iteraci)]'-b'*[X(iteraci),Y(iteraci)]';
elipsa(xx,hladina,krok,a11,a12,a22,b1,b2);

zoom on;





