%-----------------------------------------------------------
%    Geometricke znazorneni metody sdruzenych gradientu
%-----------------------------------------------------------

clc;
tisk=sprintf('Geometricke znazorneni metody sdruzenych gradientu');
disp(tisk);
disp(' ');
disp(' ');

a11=9;
a12=2;
a22=3;
b1=48;
b2=26;

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

x0=9;
y0=0;

iteraci=10;
epsilon=0.001;

presne=A\b;

figure(1);
clf;
hold on;
xmin=-1;
xmax=11;
ymin=-5;
ymax=13;
xlim=([xmin xmax]);
ylim=([ymin ymax]);
title('Metoda sdruzenych gradientu');
axis([xmin xmax ymin ymax]);
krok=0.005;

pomx=2.6:krok:6.1;
pomy=b1/a12-a11/a12*pomx;
p1=plot(pomx,pomy,'m');
pomx=0:krok:10;
pomy=b2/a22-a12/a22*pomx;
p2=plot(pomx,pomy,'b');
xx=xmin:krok:xmax;

rce1=strcat(int2str(a11),'x+',int2str(a12),'y=',int2str(b1));
rce2=strcat(int2str(a12),'x+',int2str(a22),'y=',int2str(b2));
text(0,-1,rce1);
text(0,-2,rce2);

% Metoda sdruzenych gradientu
X=x0;
Y=y0;
x=[X,Y]';
iter=0; 
g1=A*x-b;
citac=0;
nor=2*epsilon;
while (nor>epsilon) & (iter < iteraci)
   u1=g1'*g1;
   if (citac==0)  
             beta1=0;
             s1=-g1;
   else 
             beta1=u1/u0;
             s1=-g1+beta1*s0;
   end       
   z1=A*s1;
   pom=s1'*z1;
   alfa1=u1/pom;
   x=x+alfa1*s1;
   g1=g1+alfa1*z1;        
   s0=s1;
   u0=u1;
   citac=citac+1;         
   nor=norm(g1);
   X=[X,x(1)];
   Y=[Y,x(2)];
   iter=iter+1;
end

  iterace=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:iterace
 p3=plot(X(i),Y(i),'ro'); 
 p4=plot(0,0,'k');
 legend([p1 p2 p3 p4],{'1. rovnice','2. rovnice','iterace','hladina'})
 hladina=0.5*[X(i),Y(i)]*A*[X(i),Y(i)]'-b'*[X(i),Y(i)]';
 legend([p1 p2 p3 p4],{'1. rovnice','2. rovnice','iterace','hladina'})
 elipsa(xx,hladina,krok,a11,a12,a22,b1,b2);
 legend([p1 p2 p3 p4],{'1. rovnice','2. rovnice','iterace','hladina'})
 pause;
 plot([X(i),X(i+1)],[Y(i),Y(i+1)],'k--');
 legend([p1 p2 p3 p4],{'1. rovnice','2. rovnice','iterace','hladina'})
 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(iterace+1),Y(iterace+1),'ro'); 
legend([p1 p2 p3 p4],{'1. rovnice','2. rovnice','iterace','hladina'})
hladina=0.5*[X(iterace),Y(iterace)]*A*[X(iterace),Y(iterace)]'-b'*[X(iterace),Y(iterace)]';
legend([p1 p2 p3 p4],{'1. rovnice','2. rovnice','iterace','hladina'})
elipsa(xx,hladina,krok,a11,a12,a22,b1,b2);
legend([p1 p2 p3 p4],{'1. rovnice','2. rovnice','iterace','hladina'})

zoom on;





