//
// Comparaison des méthodes Jacobi-Gauss Seidel et SOR
// sur le cas de la matrice du Laplacien
//
n=evstr(x_dialog('taille de la matrice','50'));
eps=evstr(x_dialog('precision souhaitee','1E-10'));
//A=toeplitz([4,1,zeros(1,n-2)]);  // autre matrice importante a tester
A=toeplitz([2,-1,zeros(1,n-2)]);  
D=diag(diag(A));
E=-triu(A)+D;
F=-tril(A)+D;
uexact=ones(n,1);uinit=rand(n,1);errinit=max(abs(uinit-uexact));
b=A*uexact;
//
//Jacobi
//
M1=inv(D);
N1=D-A;
u=uinit;errj=errinit;
while(errj($)>eps)
u=M1*(N1*u+b);
errj=[errj;max(abs(u-uexact))];
end
clf();
plot2d((1:prod(size(errj)))',log(errj),1)
//
//Gauss Seidel
//
M2=inv(D-E);
N2=F;
u=uinit;errgs=errinit;
while(errgs($)>eps)
u=M2*(N2*u+b);
errgs=[errgs;max(abs(u-uexact))];
end
plot2d((1:prod(size(errgs)))',log(errgs),2)
//
//relaxation optimale
//
rhoj=max(abs(spec(M1*N1)));
omeg=2/(1+sqrt(1-rhoj^2));
M3=inv(1/omeg*D-E);
N3=(1-omeg)/omeg*D+F;
u=uinit;errrel=errinit;
compt=0;
while(errrel($)>eps)
compt=compt+1;
u=M3*(N3*u+b);
errrel=[errrel;max(abs(u-uexact))];
end
plot2d((1:prod(size(errrel)))',log(errrel),3)
xtitle('comparaison des vitesses de convergence Jacobi(noire), gauss Seidel(bleue), relaxation(verte)','n','erreur');
