function [x,it,Err]=grad1(x,a,b,nmax,toll)
%
% metodo del gradiente con disegno delle direzioni di discesa
%
% E' una function di tipo esplicativo, funziona per sistemi 2x2
%
% Author: Paola Gervasio 23/10/2017
% paola-gervasio.unibs.it

[n,m]=size(a);
xmin=-55;xmax=55;
ymin=-55;ymax=55;
% disegno le isolinee e le direzioni di discesa
xx=a\b;
h=.5;
xr=xmin:h:xmax;
yr=ymin:h:ymax;
[x1,y]=meshgrid(xr,yr);
z=0.5*(a(1,1)*x1.^2+2*a(1,2)*x1.*y+a(2,2)*y.^2)-b(1)*x1-b(2)*y;
f=@(x,y)0.5*(a(1,1)*x.^2+2*a(1,2)*x.*y+a(2,2)*y.^2)-b(1)*x-b(2)*y;
%xx1=x(1);yy1=x(2);
f0=f(x(1),x(2));

%
    figure(1);clf
    title('Isolinee della forma quadratica associata al sistema Ax=b');
   

    axis([xmin,xmax,ymin,ymax])
    axHndl=gca;
    set(axHndl, ...
        'XLim',[xmin xmax],'YLim',[ymin ymax], ...
        'View',[0,90]);
    axis equal
    xlabel('x');
    ylabel('y');

contour(x1,y,z,50)
title('isolinee della forma quadratica associata al sistema Ax=b')
fprintf('\n Si risolve ora il sistema lineare con il metodo del GRADIENTE\n ')
disp('Premere un tasto per continuare')
axis equal
pause
hold on
h1=plot(x(1),x(2),'mo','Markerfacecolor','m');
set(h1,'MarkerSize',6)
pause
contour(x1,y,z,[f0,f0],'k-','Linewidth',1)
axis equal
pause
%
re=[];
r=b-a*x;
bb=norm(b,2);
k=0;err=norm(r)/bb;
Err=[err];
while k<nmax & err> toll
q=a*r;
nrk=r'*r;
alpha=nrk/(r'*q);
%
h1=plot(x(1),x(2),'mo','Markerfacecolor','m');
set(h1,'MarkerSize',6)
xv=x;
%
x=x+alpha*r;
%
if k< 10
z=0.5*(a(1,1)*x1.^2+2*a(1,2)*x1.*y+a(2,2)*y.^2)-b(1)*x1-b(2)*y;
f0=f(x(1),x(2));
plot([x(1),xv(1)],[x(2),xv(2)],'m','Linewidth',2);
h1=plot(x(1),x(2),'m.');
set(h1,'MarkerSize',10)
axis equal
pause(1)
contour(x1,y,z,[f0,f0],'k-','Linewidth',1)
axis equal
pause(1)
else
plot([x(1),xv(1)],[x(2),xv(2)],'m','Linewidth',2);
h1=plot(x(1),x(2),'mo','Markerfacecolor','m');
set(h1,'MarkerSize',6)
end
%
r=r-alpha*q;

res=norm(r,2);
re=[re;res];
err=res/bb;
Err=[Err;err];
k=k+1;
end
it=k;
fprintf('\n Il metodo del GRADIENTE e` arrivato a convergenza in %3g ',it)
fprintf('iterazioni,\n fissata la tolleranza = %6.3g\n', toll)
return
