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

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

a11=40;
a12=0.1;
a22=41;
b1=-8;
b2=2;

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

x0=27;
y0=0.6;

iteraci=20;
epsilon=0.001;

presne=A\b;

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

krok=0.02;
xx=xmin:krok:xmax;

pomx=-0.29:0.001:-.11;
pomy=b1/a12-a11/a12*pomx;
plot(pomx,pomy,'m')
pomx=-35:34;
pomy=b2/a22-a12/a22*pomx;
plot(pomx,pomy,'b')

%plot([0 0], [-7 7],'m');
%plot([-40, 40], [0 0],'b');
rce1=strcat(int2str(a11),'x+',num2str(a12,1),'y=',int2str(b1));
rce2=strcat(num2str(a12,1),'x+',int2str(a22),'y=',int2str(b2));
text(-30,-29,rce1);
text(-30,-32,rce2);

% Metoda nejvetsiho spadu
X=x0;
Y=y0;
x=[X,Y]';
iter=0; 
g1=-A*x+b;
nor=2*epsilon;
while (nor>epsilon) & (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;
   nor=norm(g1);
end

  iteraci=iter;
  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;





