clear


// DONNEES 

// EXEMPLE DU LAPLACIEN 1D 
L = 1
n = 40 
dx = L/(n+1) 

c = dx*(1:1:n)'

u = 2*ones(n,1);
v = - ones(n-1,1)


A = diag(u)+diag(v,1)+diag(v,-1);
A = (1/dx**2)*A // matrice discretisant le Laplacien 1D


b = 1000*exp(-100*(c/L-1/2)**2) // terme source 

y=A\b



// METHODE DE RICHARDON A PAS OPTIMAL (CAS SDP)

x=zeros(n,1) // choix initial de x 
r = b-A*x; // residu initial 
nr0 = norm(r); // norme du residu initial 
nr = nr0;

eps = 1.0E-10; // precision sur le residu relatif 

lambda = spec(A)
alpha = 2/(lambda(1)+lambda(n));


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

scf(1)
xset("font size",5)
plot(res,'-xb')
disp('richardson a pas fixe')
disp('x',x)
disp('nr',nr)
disp('nr/nr0',nr/nr0)
disp('k',k)

// METHODE DE RICHARDSON A PAS VARIABLE (CAS SDP)

x=zeros(n,1) // choix initial de x 
r = b-A*x; // residu initial 
nr0 = norm(r); // norme du residu initial 
nr = nr0;

eps = 1.0E-10; // precision sur le residu relatif 

k = 0
while (nr/nr0>eps)&(k<10000) do 
    k = k + 1;
    res1(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

plot(res1,'-xr')
disp('richardson a pas variable')
disp('x',x)
disp('nr',nr)
disp('nr/nr0',nr/nr0)
disp('k',k)

scf(2)
xset("font size",5)
plot(c,y,'-x')

scf(3)
xset("font size",5)
plot(c,b,'-x')



