disp('Ci-contre une illustration des trois méthodes.')
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=zeros(2,1);

//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
p1=[X0*ones(1,nmax+1-size(p1,2)),p1];

//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
p2=[X0*ones(1,nmax+1-size(p2,2)),p2];

//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
p3=[X0*ones(1,nmax+1-size(p3,2)),p3];

xset("window",2), xbasc()
plot2d([p1(1,:)',p2(1,:)',p3(1,:)'],[p1(2,:)',p2(2,:)',p3(2,:)'],[2,3,5],'131',..
leg="gradient à pas constant@gradient à pas variable@gradient conjugué",..
rect=[0 0 0.5 1],nax=[1 15 1 10])
xtitle('Résolution de AX=B avec A=[2,1; 1,2], B=[1;2] et X0=[0;0]','x','y')