disp(' '), disp('Ci-contre une illustration des trois méthodes.'), disp(' ')
disp('On dessine les approximations successives de la solution, avec :')
nmax=11; tol=1.e-6;
A=[2,1; 1,2], B=[1;2], X0=[0;0];

%Méthode du gradient à pas constant
alpha=0.5; X=X0; R=B-A*X; p1=[];
for n=0:nmax
  p1=[p1,X];
  if (norm(R)<=tol), break, end
  X=X+alpha*R;
  R=B-A*X;
end

%Méthode du gradient à pas variable
X=X0; R=B-A*X; p2=[];
for n=0:nmax
  p2=[p2,X];
  if (norm(R)<=tol), break, end
  AR=A*R;
  alpha=(R'*R)/(AR'*R);
  X=X+alpha*R;
  R=R-alpha*AR;
end

%Méthode du gradient conjugué
X=X0; R=B-A*X; P=R; p3=[];
r2=R'*R;
for n=0:nmax
  p3=[p3,X];
  if (r2<=tol^2), break, end
  AP=A*P;
  alpha=r2/(AP'*P);
  X=X+alpha*P;
  R=R-alpha*AP;
  s2=R'*R;
  P=R+(s2/r2)*P;
  r2=s2;
end

plot(p1(1,:)',p1(2,:)',':b',p2(1,:)',p2(2,:)','--g',p3(1,:)',p3(2,:)','-r')
%plot(p1(1,:)',p1(2,:)',p2(1,:)',p2(2,:)',p3(1,:)',p3(2,:)')
axis equal
debut1 = '   Résolution de AX=B avec A=[2,1; 1,2], B=[1;2] et X0=[0;0] : ';
debut2 = 'gradient à pas constant en';
fin1 = 'pointillés, gradient à pas variable en traits à demi, ';
%fin1 = 'bleu, gradient à pas variable en vert,';
fin2 = 'et gradient conjugué en traits pleins.';
%fin3 = 'et gradient conjugué en rouge.';
title(strvcat([debut1,debut2],[fin1,fin2]))
xlabel('x')
ylabel('y')