clear();

//====================================================================================
//Programme de résolution numérique de différentes EDP's via
// différentes méthodes numériques
//====================================================================================

       
"1 - Definition de fonctions"

//====================================================================================
//Choix de conditions initiales:
//      - u1: Second membre pour le problème 1 et condition initiale
// pour les problèmes 2 à 6
//      - u2: Condition initiale pour l'équation de la chaleur (problèmes 2 et 3)
//      - u4: Condition initiale pour l'équation de transport (problèmes 4 à 6)
//      - u7: Position initiale pour l'équation des ondes (problème 7)
//      - v7: Vitesse initiale pour l'équation des ondes (problème 7)
//      - u8: Condition initiale pour l'équation de Schrödinger (problème 8)
//====================================================================================

function y=u1(x)
    y = exp(-25*(x-0.5).^2);
endfunction

function y=u2(x)
    y = 100*exp(-25*(x-0.5).^2);
endfunction

function y=u4(x)
    y = exp(-25*(x-0.5).^2);
endfunction

function y=u7(x)
    y = exp(-100*(x-0.5).^2);
endfunction

function y=v7(x)
    y = 0;
endfunction

function y=u8(x)
    y = sqrt(sqrt(50/%pi))*exp(-25*(x-0.5).^2);
endfunction

"2 - Choix des problèmes modèles"

disp('== Différents problèmes modèles ==');
disp('1 - Problème elliptique 1D: -u'' + Du = u1 avec oonditions aux limites de Dirichlet');
disp('  homogènes: u(0) = u(1) = 0');
disp('2 - Equation de la chaleur avec schema explicite');
disp('3 - Equation de la chaleur avec schema implicite');
disp('4 - Equation de transport avec schema Upwind explicite');
disp('5 - Equation de transport avec schema Upwind implicite');
disp('6 - Equation de transport avec schema de Lax-Wendroff');
disp('7 - Equation des ondes avec schema implicite');
disp('8 - Equation de Schrödinger avec schema implicite');
choice = input('Choisissez un modèle SVP:');

select choice
case 1
J = input('Nombre de subdivisions en espace:');
D = input('Coefficient D:');
hx = 1/J; //Pas de subdivision en espace
XX= hx:hx:1-hx;
A = sparse([(1:J-1)',(1:J-1)'],-2*ones(J-1,1)) + sparse([(1:J-1)',[(2:J-1)';1]],[ones(J-2,1);0]) + sparse([[(2:J-1)';1],(1:J-1)'],[ones(J-2,1);0]) + D*speye(J-1,J-1);
F = u1(XX)';
U = A\F;
X = [0,XX,1]';
U = hx^2*[0;U;0];
clf();
set(gca(),"grid",[1 1]);
xlabel('$x$');
ylabel('$u(x)$');
title('$\text{Solution du problème elliptique }-u'' + Du=u_0$')
plot(X,U,'color','red','thickness',2);

case 2
//Selection des données
T = input('Durée de la simulation (en s):');
kappa = input('Coefficient de diffusion thermique (en m²/s):');
J = input('Nombre de subdivisions en espace:');
hx = 1/J; //Pas de subdivision en espace
disp('Nombre minimal de subsivisions en temps pour garantir la condition CFL:');
CFL = ceil(2*kappa*T/(hx^2));
disp(CFL);
N = input('Nombre de subdivisions en temps:');
ht = T/N; //Pas de subdivision en temps
disp('Nombre CFL (ne doit pas dépasser 0.5 pour garantir la stabilité du schéma):');
CFL2 = kappa*ht/hx^2;
disp(CFL2);
disp('Pourcentage d''erreur avec la CFL autorisée (en cas de dépassement):');
disp(abs(200*CFL2-100));

//Construction des matrices du schema
XX = 0:hx:1;
TT = 0:ht:T;
A = sparse([(1:J+1)',(1:J+1)'],-2*ones(J+1,1)) + sparse([(1:J+1)',[(2:J+1)';1]],ones(J+1,1)) + sparse([[(2:J+1)';1],(1:J+1)'],ones(J+1,1)); //Matrice d'approximation de la dérivée
//seconde en espace par différence finie
M = kappa*ht/hx^2*A + speye(J+1,J+1);
U = zeros(J+1,N+1);
U(:,1) = u2(XX)';

//Calcul de la solution
for n=1:N
    V = U(:,n);
    V = M*V;
    U(:,n+1) = V;
end
    
//Tracé de la solution
f = scf(1);
plot3d1(TT,XX,U',270,0,flag=[-6,2,4]);
f.color_map = jetcolormap(512);
title('$\text{Solution du problème parabolique }\frac{\partial u}{\partial t}=\kappa \frac{\partial^2 u}{\partial x^2}$')
xlabel('$t\text{ (en s)}$');
ylabel('$x\text{ (en m)}$');
colorbar(min(U),max(U));

xclick();

f0 = scf(2);
delete();
aa=get("current_axes")
aa.data_bounds=[0 min(U) ; 1 max(U)];
set(gca(),"grid",[1 1]);
title('Evolution de la température [°C] à t = '+string(0));
for n=1:N+1
    drawlater();
    b=get("current_entity");
    delete(b.children);
    plot(XX',U(:,n),'color','red','thickness',2);
    title('Evolution de la température [°C] à t = '+string(int(100*(n-1)*ht)/100));
    drawnow();
end

case 3
//Selection des données
T = input('Durée de la simulation (en s):');
kappa = input('Coefficient de diffusion thermique (en m²/s):');
J = input('Nombre de subdivisions en espace:');
hx = 1/J; //Pas de subdivision en espace
N = input('Nombre de subdivisions en temps:');
ht = T/N; //Pas de subdivision en temps

//Construction des matrices du schema
XX = 0:hx:1;
TT = 0:ht:T;
A = sparse([(1:J+1)',(1:J+1)'],-2*ones(J+1,1)) + sparse([(1:J+1)',[(2:J+1)';1]],ones(J+1,1)) + sparse([[(2:J+1)';1],(1:J+1)'],ones(J+1,1)); //Matrice d'approximation de la dérivée
//seconde en espace par différence finie
M = -kappa*ht/hx^2*A + speye(J+1,J+1);
U = zeros(J+1,N+1);
U(:,1) = u2(XX)';

//Calcul de la solution
for n=1:N
    V = U(:,n);
    V = M\V;
    U(:,n+1) = V;
end
    
//Tracé de la solution
clf();
f = gcf();
plot3d1(TT,XX,U',270,0,flag=[-6,2,4]);
f.color_map = jetcolormap(512);
title('$\text{Solution du problème parabolique }\frac{\partial u}{\partial t}=\kappa \frac{\partial^2 u}{\partial x^2}$')
xlabel('$t\text{ (en s)}$');
ylabel('$x\text{ (en m)}$');
colorbar(min(U),max(U));

xclick();

f0 = scf(2);
delete();
aa=get("current_axes")
aa.data_bounds=[0 min(U) ; 1 max(U)];
set(gca(),"grid",[1 1]);
title('Evolution de la température [°C] à t = '+string(0));
for n=1:N+1
    drawlater();
    b=get("current_entity");
    delete(b.children);
    plot(XX',U(:,n),'color','blue','thickness',2);
    title('Evolution de la température [°C] à t = '+string(int(100*(n-1)*ht)/100));
    drawnow();
end

case 4
//Selection des données
T = input('Durée de la simulation (en s):');
c = input('Célerité (en m/s, doit rester positive):');
J = input('Nombre de subdivisions en espace:');
hx = 1/J; //Pas de subdivision en espace
disp('Nombre minimal de subsivisions en temps pour garantir la condition CFL:');
CFL = ceil(c*T/hx);
disp(CFL);
N = input('Nombre de subdivisions en temps:');
ht = T/N; //Pas de subdivision en temps
disp('Nombre CFL (ne doit pas dépasser 1 pour garantir la stabilité du schéma):');
CFL2 = c*ht/hx;
disp(CFL2);
disp('Pourcentage d''erreur avec la CFL autorisée (en cas de dépassement):');
disp(abs(100*CFL2-100));

//Construction des matrices du schema
XX = 0:hx:1;
TT = 0:ht:T;
D2 = sparse([(1:J+1)',(1:J+1)'],ones(J+1,1)) - sparse([[(2:J+1)';1],(1:J+1)'],ones(J+1,1)); //Matrice d'approximation de la dérivée en espace par différence finie rétrograde
M = -c*ht/hx*D2 + speye(J+1,J+1);
U = zeros(J+1,N+1);
U(:,1) = u4(XX)';
UL2 = zeros(N+1,1); //Vecteur de l'évolution de la norme L2 (relativement à la solution initiale, en pourcentage)
UL2(1) = 100;

//Calcul de la solution
for n=1:N
    V = U(:,n);
    V = M*V;
    U(:,n+1) = V;
    UL2(n+1) = 100*sqrt(sum(U(:,n+1).^2)/sum(U(:,1).^2));
end

//Tracé de la solution
f = scf(1);
plot3d1(TT,XX,U',270,0,flag=[-6,2,4]);
f.color_map = jetcolormap(512);
title('$\text{Solution du problème hyperbolique }\frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0$')
xlabel('$t\text{ (en s)}$');
ylabel('$x\text{ (en m)}$');
colorbar(min(U),max(U));

g = scf(2);
xlabel('$t\text{ (en s))}$');
ylabel('$\text{Pourcentage de }\left|\left|u_0\right|\right|_{L^2}$')
title('Evolution relative de la norme L2 de la solution calculée du problème hyperbolique');
set(gca(),"grid",[1 1]);
plot(TT',UL2,'color','red','thickness',2);

xclick();

f0 = scf(3);
delete();
aa=get("current_axes")
aa.data_bounds=[0 min(U) ; 1 max(U)];
set(gca(),"grid",[1 1]);
title('Evolution de la solution à t = '+string(0));
for n=1:N+1
    drawlater();
    b=get("current_entity");
    delete(b.children);
    plot(XX',U(:,n),'color','blue','thickness',2);
    title('Evolution de la solution à t = '+string(int(100*(n-1)*ht)/100));
    drawnow();
end

case 5
//Selection des données
T = input('Durée de la simulation (en s):');
c = input('Célerité (en m/s):');
J = input('Nombre de subdivisions en espace:');
hx = 1/J; //Pas de subdivision en espace
N = input('Nombre de subdivisions en temps:');
ht = T/N; //Pas de subdivision en temps

//Construction des matrices du schema
XX = 0:hx:1;
TT = 0:ht:T;
D2 = sparse([(1:J+1)',(1:J+1)'],ones(J+1,1)) - sparse([[(2:J+1)';1],(1:J+1)'],ones(J+1,1)); //Matrice d'approximation de la dérivée en espace par différence finie rétrograde
M = c*ht/hx*D2 + speye(J+1,J+1);
U = zeros(J+1,N+1);
U(:,1) = u4(XX)';
UL2 = zeros(N+1,1); //Vecteur de l'évolution de la norme L2 (relativement à la solution initiale, en pourcentage)
UL2(1) = 100;

//Calcul de la solution
for n=1:N
    V = U(:,n);
    V = M\V;
    U(:,n+1) = V;
    UL2(n+1) = 100*sqrt(sum(U(:,n+1).^2)/sum(U(:,1).^2));
end

//Tracé de la solution
f = scf(1);
plot3d1(TT,XX,U',270,0,flag=[-6,2,4]);
f.color_map = jetcolormap(512);
title('$\text{Solution du problème hyperbolique }\frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0$')
xlabel('$t\text{ (en s)}$');
ylabel('$x\text{ (en m)}$');
colorbar(min(U),max(U));

g = scf(2);
xlabel('$t\text{ (en s))}$');
ylabel('$\text{Pourcentage de }\left|\left|u_0\right|\right|_{L^2}$')
title('Evolution relative de la norme L2 de la solution calculée du problème hyperbolique');
set(gca(),"grid",[1 1]);
plot(TT',UL2,'color','red','thickness',2);

xclick();

f0 = scf(3);
delete();
aa=get("current_axes")
aa.data_bounds=[0 min(U) ; 1 max(U)];
set(gca(),"grid",[1 1]);
title('Evolution de la solution à t = '+string(0));
for n=1:N+1
    drawlater();
    b=get("current_entity");
    delete(b.children);
    plot(XX',U(:,n),'color','blue','thickness',2);
    title('Evolution de la solution à t = '+string(int(100*(n-1)*ht)/100));
    drawnow();
end

case 6
//Selection des données
T = input('Durée de la simulation (en s):');
c = input('Célerité (en m/s, doit rester positive):');
J = input('Nombre de subdivisions en espace:');
hx = 1/J; //Pas de subdivision en espace
disp('Nombre minimal de subsivisions en temps pour garantir la condition CFL:');
CFL = ceil(c*T/hx);
disp(CFL);
N = input('Nombre de subdivisions en temps:');
ht = T/N; //Pas de subdivision en temps
disp('Nombre CFL (ne doit pas dépasser 1 pour garantir la stabilité du schéma):');
CFL2 = c*ht/hx;
disp(CFL2);
disp('Pourcentage d''erreur avec la CFL autorisée (en cas de dépassement):');
disp(abs(100*CFL2-100));

//Construction des matrices du schema
XX = 0:hx:1;
TT = 0:ht:T;
D = sparse([(1:J+1)',[(2:J+1)';1]],ones(J+1,1)) - sparse([[(2:J+1)';1],(1:J+1)'],ones(J+1,1)); //Matrice d'approximation de la dérivée en espace par différence finie centrée
A = sparse([(1:J+1)',(1:J+1)'],-2*ones(J+1,1)) + sparse([(1:J+1)',[(2:J+1)';1]],ones(J+1,1)) + sparse([[(2:J+1)';1],(1:J+1)'],ones(J+1,1)); //Matrice d'approximation de la dérivée
//seconde en espace par différence finie
M = speye(J+1,J+1) - (0.5*c*ht/hx)*D + 0.5*(c*ht/hx)^2*A;
U = zeros(J+1,N+1);
U(:,1) = u4(XX)';
UL2 = zeros(N+1,1); //Vecteur de l'évolution de la norme L2 (relativement à la solution initiale, en pourcentage)
UL2(1) = 100;

//Calcul de la solution
for n=1:N
    V = U(:,n);
    V = M*V;
    U(:,n+1) = V;
    UL2(n+1) = 100*sqrt(sum(U(:,n+1).^2)/sum(U(:,1).^2));
end

//Tracé de la solution
f = scf(1);
plot3d1(TT,XX,U',270,0,flag=[-6,2,4]);
f.color_map = jetcolormap(512);
title('$\text{Solution du problème hyperbolique }\frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0$')
xlabel('$t\text{ (en s)}$');
ylabel('$x\text{ (en m)}$');
colorbar(min(U),max(U));

g = scf(2);
xlabel('$t\text{ (en s))}$');
ylabel('$\text{Pourcentage de }\left|\left|u_0\right|\right|_{L^2}$')
title('Evolution relative de la norme L2 de la solution calculée du problème hyperbolique');
set(gca(),"grid",[1 1]);
plot(TT',UL2,'color','red','thickness',2);

xclick();

f0 = scf(3);
delete();
aa=get("current_axes")
aa.data_bounds=[0 min(U) ; 1 max(U)];
set(gca(),"grid",[1 1]);
title('Evolution de la solution à t = '+string(0));
for n=1:N+1
    drawlater();
    b=get("current_entity");
    delete(b.children);
    plot(XX',U(:,n),'color','blue','thickness',2);
    title('Evolution de la solution à t = '+string(int(100*(n-1)*ht)/100));
    drawnow();
end


case 7
//Selection des données
T = input('Durée de la simulation (en s):');
c = input('Célerité (en m/s):');
J = input('Nombre de subdivisions en espace:');
hx = 1/J; //Pas de subdivision en espace
N = input('Nombre de subdivisions en temps:');
ht = T/N; //Pas de subdivision en temps

//Données initiales

//Construction des matrices du schema
XX = 0:hx:1;
TT = 0:ht:T;
A = sparse([(1:J+1)',(1:J+1)'],-2*ones(J+1,1)) + sparse([(1:J+1)',[(2:J+1)';1]],ones(J+1,1)) + sparse([[(2:J+1)';1],(1:J+1)'],ones(J+1,1)); //Matrice d'approximation de la dérivée
//seconde en espace par différence finie;
M = sparse([[spzeros(J+1,J+1),speye(J+1,J+1)];[(c*ht/hx)^2*A-speye(J+1,J+1),2*speye(J+1,J+1)]]);
U = zeros(2*J+2,N+1);
U(:,1) = [u7(XX),u7(XX)+ht*v7(XX)]';

//Calcul de la solution
for n=1:N
    V = U(:,n);
    V = M\V;
    U(:,n+1) = V;
end

//Tracé de la solution
f = scf(1);
plot3d1(TT,XX,U(1:J+1,:)',270,0,flag=[-6,2,4]);
f.color_map = jetcolormap(512);
title('$\text{Solution du problème hyperbolique }\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}$')
xlabel('$t\text{ (en s)}$');
ylabel('$x\text{ (en m)}$');
colorbar(min(U),max(U));

xclick();

f0 = scf(2);
delete();
aa=get("current_axes")
aa.data_bounds=[0 min(U) ; 1 max(U)];
set(gca(),"grid",[1 1]);
title('Evolution de la solution à t = '+string(0));
for n=1:N+1
    drawlater();
    b=get("current_entity");
    delete(b.children);
    plot(XX',U(1:J+1,n),'color','blue','thickness',2);
    title('Evolution de la solution à t = '+string(int(100*(n-1)*ht)/100));
    drawnow();
end

case 8
//Selection des données
T = input('Durée de la simulation (en s):');
m = input('Masse de la particule (en kg):');
J = input('Nombre de subdivisions en espace:');
hx = 1/J; //Pas de subdivision en espace
N = input('Nombre de subdivisions en temps:');
ht = T/N; //Pas de subdivision en temps
hbar = 6.63*10**-34/2*%pi //Constante de Planck réduite

//Construction des matrices du schema
XX = 0:hx:1;
TT = 0:ht:T;
A = sparse([(1:J+1)',(1:J+1)'],-2*ones(J+1,1)) + sparse([(1:J+1)',[(2:J+1)';1]],ones(J+1,1)) + sparse([[(2:J+1)';1],(1:J+1)'],ones(J+1,1)); //Matrice d'approximation de la dérivée
//seconde en espace par différence finie;
M = -%i*hbar*ht/(2*m*hx^2)*A + speye(J+1,J+1);
U = zeros(J+1,N+1);
U(:,1) = u8(XX)';
UL2 = zeros(N+1,1);
UL2(1) = hx*sum(abs(U(:,1).^2));

//Calcul de la solution
for n=1:N
    V = U(:,n);
    V = M\V;
    U(:,n+1) = V;
    UL2(n+1) = hx*sum(abs(U(:,n+1).^2));
end

ReU=real(U);
ImU=imag(U);
AbsU=abs(U).^2;

//Tracé de la solution

f = scf(1);
plot3d1(TT,XX,AbsU',270,0,flag=[-6,2,4]);
f.color_map = jetcolormap(512);
title('$\text{Carré du module de la solution de l équation de Schrödinger }\frac{\partial u}{\partial t}=i\frac{\hbar}{2m} \frac{\partial^2 u}{\partial x^2}$')
xlabel('$t\text{ (en s)}$');
ylabel('$x\text{ (en m)}$');
colorbar(min(AbsU),max(AbsU));

xclick();

g = scf(2);
xlabel('$t\text{ (en s))}$');
ylabel('$\text{Valeur de}\left|\left|\Psi(\cdot,t)\right|\right|_{L^2}$')
title('Evolution du carré de la norme L2 de la fonction d''onde');
set(gca(),"grid",[1 1]);
plot(TT',UL2,'color','red','thickness',2);

xclick();

clf();
delete();
aa=get("current_axes")
aa.data_bounds=[0 min(AbsU) ; 1 max(AbsU)];
set(gca(),"grid",[1 1]);
title('Evolution de la solution (valeur absolue) à t = '+string(0));
for n=1:N+1
    drawlater();
    b=get("current_entity");
    delete(b.children);
    plot(XX',AbsU(:,n).^2,'color','green','thickness',2);
    title('Evolution de la solution (carré de la valeur absolue) à t = '+string(int(100*(n-1)*ht)/100));
    drawnow();
end

end
