///////////////////////////////////////
//  Résolution de l'équation de la   //
// chaleur 1D par différences finies //
///////////////////////////////////////


////////////Choix des paramètres///////////////////
N=50;//Nombre de pas de discrétisation en espace

dx=1/(N+1);//Pas d'espace
dt=dx^2*0.1;//Pas de temps
a=0;b=1;//Conditions au bord


lambda=dt/dx^2;//Le paramètre de stabilité

if lambda<0.5//Test de la stabilité du schéma
  messagebox('La quantité dt/dx^2='+string(dt/dx^2)+' est inférieure à 0.5. Le schéma d''Euler sera donc stable.');
else
  messagebox('La quantité dt/dx^2='+string(dt/dx^2)+' est supérieur à 0.5. Le schéma d''Euler explicite sera donc instable.');
end;

x=[0:N+1]*dx;//Grille en espace


///////Définition de la condition initiale//////////////
//On prend une condition initiale aléatoire
sol_euler_exp=[grand(N,1,'nor',0,1)];
sol_euler_imp=sol_euler_exp;
sol_CN=sol_euler_exp;


/////////La matrice du schéma////////////////
A=-2*diag(ones(1,N))+diag(ones(1,N-1),1)+diag(ones(1,N-1),-1);
A_euler_exp=eye(N,N)+lambda*A;
A_euler_imp=eye(N,N)-lambda*A;
A_CN_1=eye(N,N)+lambda*A/2;
A_CN_2=eye(N,N)-lambda*A/2;


//////Configuration des fenêtres graphiques///////////////////////
f=scf(0);//On ouvre une fenêtre
clf();
f.pixmap="on";//Pour faire une animation


//////////Le schéma à proprement parler//////////////////////////////:
iteration=0;
log_norme=log(dx*sum(sol_euler_imp.^2))/2;
while 1==1
  scf(0);
  clf();
  subplot(2,2,1);//Euler explicite
  title("Schéma d''Euler explicite");
  plot2d(x,[a ; sol_euler_exp ; b],5);
  subplot(2,2,2);//Euler implicite
  title("Schéma d''Euler implicite");
  plot2d(x,[a ; sol_euler_imp ; b],5);
  subplot(2,2,3);//Crank-Nicolson
  title("Schéma de Crank-Nicolson");
  plot2d(x,[a ; sol_CN ; b],5);
  show_pixmap();

  sol_euler_exp=A_euler_exp*sol_euler_exp+lambda*[a ; zeros(N-2,1) ; b];//Schéma d'Euler explicite
  sol_euler_imp=A_euler_imp\(sol_euler_imp+lambda*[a ; zeros(N-2,1) ; b]);//Schéma d'Euler implicite
  sol_CN=A_CN_2\(A_CN_1*(sol_CN+lambda*[a ; zeros(N-2,1) ; 1]/2))+lambda*[a ; zeros(N-2,1) ; b]/2;//Schéma de Crank-Nicolson
end;
