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


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

dx=1/(N+1);//Pas d'espace
dt=0.00002;//Pas de temps

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


x=[0:N+1]*dx;//Grile en espace
X=x'*ones(x);//Matrice des abscisses
Y=ones(x')*x;//Matrice des ordonnées
//(On se place sur le carré [0,1]x[0,1])



///////Définition de la condition initiale//////////////
//On définit la condition initiale comme une somme bruitée
//d'indicatrices de disques
A=0.6*((X-0.8).^2+(Y-0.2).^2<0.05);
A=A+0.2*((X-0.5).^2+(Y-0.5).^2<0.05);
A=A+0.7*((X-0.8).^2+(Y-0.9).^2<0.01);
A=A+0.7*((X-0.1).^2+(Y-0.9).^2<0.004);

A=A+grand(N+2,N+2,'nor',0.5,0.1);//On rajoute du bruit

A(1,:)=zeros(1,N+2);//Condition au bord
A($,:)=ones(1,N+2);
A(:,1)=x';
A(:,$)=x';

a=min(A);
b=max(A);



//////Configuration des fenêtres graphiques///////////////////////
f=scf(0);//On ouvre une fenêtre
clf();
nbre_coul=256;
f.color_map=jetcolormap(nbre_coul);
f.pixmap="on";//Pour faire une animation
subplot(1,2,1);
plot3d1(x,x,A);//Affichage de la solution (comme surface)
subplot(1,2,2);
Matplot((nbre_coul-1)*(A-a)/(b-a)+1);//Affichage de la solution (avec des couleurs)
show_pixmap();


//////////Le schéma à proprement parler//////////////////////////////:
while 1==1
  //On calcule le laplacien discret de A
  lapl_A=  (  A([2:$-1],[1:$-2]) + A([2:$-1],[3:$]) + A([3:$],[2:$-1]) + A([1:$-2],[2:$-1]) - 4*A([2:$-1],[2:$-1])  ) / (dx^2);
  //On préserve les conditions nulles au bord
  lapl_A=[ zeros(1,N+2) ; zeros(N,1) lapl_A zeros(N,1) ; zeros(1,N+2)];
  A=A+dt*lapl_A;//Schéma d'Euler explicite

  scf(0);
  clf();
  subplot(1,2,1);//Affichage en surface;
  plot3d1(x,x,A);
  subplot(1,2,2);//Affichage en couleurs
  Matplot((nbre_coul-1)*(A-a)/(b-a)+1);
  show_pixmap();
end;
