clear;
for ni=1:4 // boucle sur la taille n des matrices 



n = 10*ni // taille de la matrice
nn(ni) = n
h = 1/(n+1) 

u = 2*ones(n,1);
v = - ones(n-1,1)
A = diag(u)+diag(v,1)+diag(v,-1);
A = (1/h**2)*A 


b = rand(n,1)  

y=A\b

// algorithme de GAUSS SEIDEL
x = zeros(n,1);
eps = 1.d-10 ;
kmax = 100000 ;
r = b-A*x
nr0 = norm(r)
nr = nr0
k = 1

while (nr/nr0>eps)&(k<kmax) then 
    
    k = k +1 
    res(k) = log(nr/nr0)/log(10);
     
   
 //   for i=1:n
 //       s = 0
 //       for j=1:i-1
 //           s = s + A(i,j)*x(j)
 //       end
 //       for j=i+1:n
 //           s = s + A(i,j)*x(j)
 //       end
 //       
 //       x(i) = (b(i)-s)/A(i,i)
 //   end 
   
    for i=1:n 
     dxi = ( b(i) - A(i,:)*x )/A(i,i)
     x(i) = x(i) + dxi
    end
    
    r = b-A*x
    nr = norm(r)
    
end

nit1(ni) = k

disp('erreur ',norm(x-y))
disp('nit',k)

scf(1)
xset("font size",5)
plot(res,'-xb')


// Algorithme SOR

x = zeros(n,1);
eps = 1.d-10 ;
kmax = 100000 ;
r = b-A*x
nr = norm(r)
nr0 = nr
k = 1
omega = 1.2
while (nr/nr0>eps)&(k<kmax) then 
     
    k = k +1 
    res1(k) = log(nr/nr0)/log(10);
    
    for i=1:n 
     dxi = ( b(i) - A(i,:)*x )/A(i,i) 
     xti = x(i) + dxi
     x(i) = omega*xti + (1-omega)*x(i) 
    end    
   
 
    r = b-A*x
    nr = norm(r)
    
end

nit2(ni) = k

disp('erreur ',norm(x-y))
disp('nit',k)

scf(1)
xset("font size",5)
plot(res1,'-xr')



// Methode de Richarson a pas variable 
eps = 1.d-10 ;
kmax = 100000 ;
x=zeros(n,1) // choix initial de x 
r = b-A*x; // residu initial 
nr0 = norm(r); // norme du residu initial 
nr = nr0;

k = 0
while (nr/nr0>eps)&(k<kmax) do 
    k = k + 1;
    res2(k) = log(nr/nr0)/log(10);
    p = A*r;
    alpha = (r'*r)/(r'*p);
    x = x + alpha*r;
    r = r - alpha*p;
    nr = norm(r);   
end

nit3(ni) = k

disp('erreur ',norm(x-y))
disp('nit',k)

scf(1)
xset("font size",5)
plot(res2,'-xm')


// la pente beta des courbes suivantes donne le comportement en C.n^beta du nombre d'iterations des algorithmes Richardson a pas variable, GS et SOR (on trouve grosso modo beta=2 pour les trois algorithmes. 

// Par contre au vu du nombre d'itérations, la constante C est plus petite pour SOR (avec omega=1.2) que pour GS et plus petite pour GS que pour Richardson à pa variable )
// le meilleur algorithme est donc SOR avec omega=1.2 pour cette famille de matrices SDP tridiagonales 
scf(5)
xset("font size",5)
plot(log(nn/nn(1)),log(nit1/nit1(1)),'-xb')
plot(log(nn/nn(1)),log(nit2/nit2(1)),'-xr')
plot(log(nn/nn(1)),log(nit3/nit3(1)),'-xm')

end // fin boucle sur la taille n des matrices 
