////////////////////////////////////////////////////////////////////////////////
// résolution numérique d'une équation de transport 1D 
// du/dt+v du/dx=0   u=u(t,x)    u|[t=0]=u0
////////////////////////////////////////////////////////////////////////////////
driver("X11");  
//chdir('C:\Documents and Settings\philippe\Mes documents\recherche\simulations\transport\gif');
////////////////////////////////////////////////////////////////////////////////
// discrétisation espace 
///////////////////////////////////////////////////////////////////////////////
X=10;// intervalle [-X,X]
dx=0.02;// pas en espace 
K=2*(X/dx)+1;//nombre de points
x=[-X:dx:X]';// vecteur contenant les abcisses
dt=0.01;//pas de temps
///////////////////////////////////////////////////////////////////////////////
// le champ des vitesses v(x) et la donnée initiale u0    
////////////////////////////////////////////////////////////////////////////// 
u=%e.^(-(x+8).^2); //densité initiale : une gaussienne à droite de la fenêtre
w=%e.^(-(x-4).^2)-0.5*%e.^(-(x+3).^2);v=1+w; w=0.5*w;//vitesse variable 
//////////////////////////////////////////////////////////////////////////////
//lancement affichage
//////////////////////////////////////////////////////////////////////////////
Y=1.5;rect=[-X,-0.5,X,Y];
titlepage(["résolution numérique d EDP";" ";"équation de transport 1D";" ";"du/dt+v(x) du/dx=0    u|[t=0]=exp(-(x+8)^2)"])
xtitle("(cliquer dans la fenêtre pour lancer l animation)")
[c_i,c_x,c_y,c_w]=xclick();
xbasc();
xset("pixmap",1)
xtitle("équation de transport 1D (cliquer pour lancer)","x","amplitude")
plot2d([x x],[u w], [2 3],"111","du/dt+v(x) du/dx=0@ différentiel de vitesse :v(x)-1",rect)
xset("wshow")
[c_i,c_x,c_y,c_w]=xclick();
/////////////////////////////////////////////////////////////////////////////////
//calcul récursif de l'évolution 
/////////////////////////////////////////////////////////////////////////////////
for i=1:1600
rightu=[0;u(1:K-1)];
// calcul de u au pas de temps suivant  
u=u+((dt/dx)*v).*(rightu-u);
u(K)=0;u(1)=0;
//affichage du résultat
xclea(-X,Y,2*X,Y+0.5);
plot2d([x x],[u w], [2 3],"000");
xstring(1,1.3,"t="+string(i*dt));
xset("wshow")
end
xset("pixmap",0)