clear all;close all;
% example of fill in in gaussian elimination

n=30;
G=numgrid('S',n);
A=delsq(G);
spy(A)
[L,U]=lu(A);
spy(L)
spy(U)
A(1,:)=ones(1,1:n^2)
A(:,1)=ones(1:n2,1)
spy(A)
[L,U]=lu(A);
spy(L)
spy(U)

% example of improvement for SOR, 
clear all
tol=1e-6;
n=30;
G=numgrid('S',n);
A=delsq(G);
lmax=8*(cos(pi/2/(n-1)))^2;
lmin=8*(sin(pi/2/(n-1)))^2;
rhoJ=(lmax-4)/4;rhoGS=rhoJ^2;
wopt=2/(1+sqrt(1-rhoJ^2));
rhoSOR=wopt-1;
[rhoJ rhoGS rhoSOR]
b=rand((n-2)^2,1);
x0=zeros(size(b));
[x,resJ]=Jacobi(A,b,x0,tol,20000);
[x,resGS]=SOR(A,b,1,x0,tol,10000);
[x,resSOR]=SOR(A,b,wopt,x0,tol,1000);
%maintenant Richardson
[x,resR]=Richardson(A,b,1/4,x0,tol,10000);
aopt=2/(lmin+lmax);
rhoR=(lmax-lmin)/(lmax+lmin);
[x,resRopt]=Richardson(A,b,aopt,x0,tol,10000);

semilogy(0:length(resJ)-1,resJ,'-',0:length(resGS)-1,resGS,'-',...
    0:length(resSOR)-1,resSOR,'-',0:length(resRopt)-1,resRopt,'-')
xlabel('iteration');
ylabel('residual');
legend('Jacobi','Gauss Seidel','SOR','Richardson');

