/// Résolution du problème de Poisson
/// avec conditions de Dirichlet homogènes
/// -u''(x) = f(x) sur (0,1)
///  u(0) = u(1) =0
/// avec divers termes sources


clear();
stacksize(100000000);

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


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

/// Termes sources
F1 = sin(%pi*x);
F2 = 1-abs(1-2*x);
F3 = (x>=1/4).*(x<=3/4);

/// Solutions exactes associées
SolEx1b = sin(%pi*xb)/(%pi^2);
SolEx2b = (xb/4-xb.^3/3).*(xb>=0).*(xb<=1/2)+(-1/12+3*xb/4-xb.^2+xb.^3/3).*(xb<=1).*(xb>=1/2);
SolEx3b = (xb/4).*(xb>=0).*(xb<=1/4)+(-1/32+xb/2-xb.^2/2).*(xb>=1/4).*(xb<=3/4)+(1/4-xb/4).*(xb<=1).*(xb>=3/4);

/// Résolution du problème de Poisson
A = Laplacien(J);
[AA,rk] = lufact(A);

Sol1 = A\(lamb*F1);
Sol1LU = lusolve(AA,lamb*F1);

Sol2 = A\(lamb*F2);
Sol2LU = lusolve(AA,lamb*F2);

Sol3 = A\(lamb*F3);
Sol3LU = lusolve(AA,lamb*F3);

/// Visualtion des solutions
figure(1)
subplot(211)
plot(xb,[0;Sol1;0],'d',xb,SolEx1b)
xlabel("x")
subplot(212)
plot(xb,abs([0;Sol1;0]-SolEx1b))
xlabel("x")
title("Erreur F1")

figure(2)
subplot(211)
plot(xb,[0;Sol1LU;0],'d',xb,SolEx1b)
xlabel("x")
subplot(212)
plot(xb,abs([0;Sol1LU;0]-SolEx1b))
xlabel("x")
title("Erreur F1")

figure(3)
subplot(211)
plot(xb,[0;Sol2;0],'d',xb,SolEx2b)
xlabel("x")
subplot(212)
plot(xb,abs([0;Sol2;0]-SolEx2b))
xlabel("x")
title("Erreur F2")


figure(4)
subplot(211)
plot(xb,[0;Sol3;0],'d',xb,SolEx3b)
xlabel("x")
subplot(212)
plot(xb,abs([0;Sol3;0]-SolEx3b))
xlabel("x")
title("Erreur F3")



/// 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)+1);
    xk = (1:JJ(k))'*ddx;    
    lambk = ddx^2;  
    Ak = Laplacien(JJ(k));
    F1k = sin(%pi*xk);
    Sol1k = Ak\(lambk*F1k);
    SolEx1k = sin(%pi*xk)/(%pi^2);
    ErrInf(1,k) = max(abs(Sol1k-SolEx1k));
    Err2(1,k) = norm(abs(Sol1k-SolEx1k),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+1)),log(ErrInf),1)
cf2 = polyfit(log(1 ./(JJ+1)),log(Err2),1)
Xmin = min(log(1 ./(JJ+1)));
Xmax = max(log(1 ./(JJ+1)));
t = linspace(Xmin,Xmax);

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