////////////////////////////////////////////////////////////////////////////////
// résolution numérique d'une équation de Burger 1D 
// du/dt+u du/dx=0   u=u(t,x)    u|[t=0]=u0
////////////////////////////////////////////////////////////////////////////////
driver("X11");  
////////////////////////////////////////////////////////////////////////////////
// discrétisation espace 
///////////////////////////////////////////////////////////////////////////////
X=20;// 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
///////////////////////////////////////////////////////////////////////////////
// la donnée initiale u0    
////////////////////////////////////////////////////////////////////////////// 
u=1-(2/%pi)*atan(x);fin=u(1);//choc
//u=bool2s(x>0)+0.5;fin=u(1);//détente
//////////////////////////////////////////////////////////////////////////////
//lancement affichage
//////////////////////////////////////////////////////////////////////////////
Y=2;rect=[-X,-0.5,X,Y];
titlepage(["résolution numérique d EDP";" ";"équation de Burger 1D";" ";"du/dt+u du/dx=0  u0=1-(2/%pi)*arctan(x)"])
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 burger 1D (cliquer pour lancer)","x","amplitude")
plot2d(x,u, 2,"111","du/dt+u du/dx=0",rect)
xset("wshow")
[c_i,c_x,c_y,c_w]=xclick();
/////////////////////////////////////////////////////////////////////////////////
//calcul récursif de l'évolution 
/////////////////////////////////////////////////////////////////////////////////
for i=1:1000
rightu=[0;u(1:K-1)];
// calcul de u au pas de temps suivant 
u=u+((dt*u)./(dx)).*(rightu-u);
u(1)=fin;// condition au bord à droite
//affichage du résultat
xclea(-X,Y,2*X,Y+0.5)
plot2d(x,u,2,"000")
xstring(15,1,"t="+string(i*dt));
xset("wshow")
end
xset("pixmap",0)