/// Résolution du problème elliptique
/// avec conditions de Neumann homogènes
/// -u''(x)+u(x) = f(x) sur (0,1)
///  u'(0) = u'(1) =0
/// avec terme source f(x) = cos(pi*x)


clear();
stacksize(100000000);

L = 1;            // taille du domaine [0;L]
J = 20;           // nombre de mailles en espace
dx= L/J;      // pas en espace
x = (0:J)'*dx;    // points de la grille 
lamb = dx^2;      


/// Assemblage de la matrice du Laplacien discret
/// avec conditions de bord de Neumann
function A = LaplacienN(J) 
    A = sparse(2*eye(J+1,J+1)-diag(ones(J,1),1)-diag(ones(J,1),-1));
    A(1,2) = -2;
    A(J+1,J) = -2;
endfunction

/// Terme source
F = cos(%pi*x);

/// Solution exacte associée
SolEx = cos(%pi*x)/(1+%pi^2);


/// Résolution du problème de Poisson
B = LaplacienN(J)+lamb*sparse(eye(J+1,J+1));
[BB,rk] = lufact(B);

Sol = B\(lamb*F);
SolLU = lusolve(BB,lamb*F);


/// Visualtion des solutions
figure(1)
subplot(211)
plot(x,Sol,'d',x,SolEx)
xlabel("x")
subplot(212)
plot(x,abs(Sol-SolEx))
xlabel("x")
title("Erreur F")

figure(2)
subplot(211)
plot(x,SolLU,'d',x,SolEx)
xlabel("x")
subplot(212)
plot(x,abs(SolLU-SolEx))
xlabel("x")
title("Erreur F")




/// Calcul de l'erreur en fonction de la discrétisation
JJ = [20 50 100 200 300];
ErrInf = zeros(1,5);
Err2 = zeros(1,5);

for k=1:5
    ddx= 1/(JJ(k));
    xk = (0:JJ(k))'*ddx;    
    lambk = ddx^2;  
    Ak = LaplacienN(JJ(k))+lambk*sparse(eye(JJ(k)+1,JJ(k)+1));
    Fk = cos(%pi*xk);
    Solk = Ak\(lambk*Fk);
    SolExk = cos(%pi*xk)/(1+%pi^2);
    ErrInf(1,k) = max(abs(Solk-SolExk));
    Err2(1,k) = norm(abs(Solk-SolExk),2);
end

/// Fonction qui permet de faire une meilleure approximation à un polynome
/// de degré donné en utilisant la méthode des moindres carrés
function cf = polyfit(x,y,n)
AA = ones(length(x),n+1)
for i=1:n
    AA(:,i+1) = x(:).^i
end
cf = lsq(AA,y(:))
endfunction

cfInf = polyfit(log(1 ./(JJ)),log(ErrInf),1)
cf2 = polyfit(log(1 ./(JJ)),log(Err2),1)
Xmin = min(log(1 ./(JJ)));
Xmax = max(log(1 ./(JJ)));
t = linspace(Xmin,Xmax);

figure(3)
subplot(211)
plot(log(1 ./(JJ)),log(ErrInf),'d',t,cfInf(1)+cfInf(2)*t)
subplot(212)
plot(log(1 ./(JJ)),log(Err2),'d',t,cf2(1)+cf2(2)*t)
