/////////////////////////////////////////////////////////////
// Copyright G. Allaire, Octobre 2002, Janvier 2012
//
// schemas numeriques sur l'equation de Burgers avec 
// conditions aux limites periodiques
// u,t + (u^2/2),x = 0 
// 
/////////////////////////////////////////////////////////////
//
xset("font",2,3) ; 
xset("thickness",2) ;
//
pi = %pi ;
//
// lg = longueur du domaine
lg = 1. ;
// nx = nombre de mailles
nx = 100 ;
// dx = pas d'espace
dx = lg/nx ;
cfl = dx ;
// dt = pas de temps
// (on sait que la vitesse maximale sera 1.)
dt = 0.9*cfl ;
// nt = nombre de pas de temps effectues
nt = 300 ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
// x = points du maillage
// u = donnee initiale
for i=1:nx
  x(i) = (i-1)*dx  ;
  u(i) = sin(x(i)*2*pi) ;
end
w = u ;

// affichage de la donnee initiale
clf()
tics=[4,10,4,10];
plotframe([0,-1.1,lg,1.1],tics);
plot2d(x,u,1,"000")
xtitle ('Lax-Friedrichs, cfl=0.9' ,' ',' ') ;

halt() ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Lax-Friedrichs
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:40
//
up = shiftp('+1',u) ;
um = shiftp('-1',u) ;
v  = (up+um)/2 - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
u = v ; 
if modulo(n,2) == 0
  clf() 
  plotframe([0,-1.3,lg,1.3],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Lax-Friedrichs")
  xtitle('Lax-Friedrichs, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

for n=41:nt
//
up = shiftp('+1',u) ;
um = shiftp('-1',u) ;
v  = (up+um)/2 - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
u = v ; 
if modulo(n,2) == 0
  clf() 
  plotframe([0,-1.1,lg,1.1],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Lax-Friedrichs")
  xtitle('Lax-Friedrichs, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// nouveau schema: Lax-Wendroff avec cfl=0.9
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 0.9*cfl ;
nt = 300 ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  u(i) = sin(x(i)*2*pi) ;
end
w = u ;

clf()
tics=[4,10,4,10];
plotframe([0,-1.1,lg,1.1],tics);
plot2d(x,u,1,"000")
xtitle ('Lax-Wendroff, cfl=0.9' ,' ',' ') ;

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Lax-Wendroff
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:40
//
up = shiftp('+1',u) ;
um = shiftp('-1',u) ;
v = u - (0.5*dt/dx)*(up.^2-um.^2)/2 + ... 
(0.5*dt*dt/(dx*dx))*( (u+up)*0.5.*(up.^2-u.^2) ...
                    - (u+um)*0.5.*(u.^2-um.^2) )/2 ;
u = v ; 
if modulo(n,2) == 0
  clf() 
  plotframe([0,-1.1,lg,1.1],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Lax-Wendroff")
  xtitle('Lax-Wendroff, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

for n=41:nt
//
up = shiftp('+1',u) ;
um = shiftp('-1',u) ;
v = u - (0.5*dt/dx)*(up.^2-um.^2)/2 + ... 
(0.5*dt*dt/(dx*dx))*( (u+up)*0.5.*(up.^2-u.^2) ...
                    - (u+um)*0.5.*(u.^2-um.^2) )/2 ;
u = v ; 
if modulo(n,2) == 0
  clf() 
  plotframe([0,-1.1,lg,1.1],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Lax-Wendroff")
  xtitle('Lax-Wendroff, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// nouveau schema: centre explicite avec cfl=0.3 (explose)
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 0.3*cfl ;
nt = 300 ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  u(i) = sin(x(i)*2*pi) ;
end
w = u ;

clf()
tics=[4,10,4,10];
plotframe([0,-1.1,lg,1.1],tics);
plot2d(x,u,1,"000")
xtitle ('centre explicite, cfl=0.3' ,' ',' ') ;

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : centre explicite
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:120
//
up = shiftp('+1',u) ;
um = shiftp('-1',u) ;
v = u - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
u = v ; 
if modulo(n,2) == 0
  clf() 
  plotframe([0,-1.1,lg,1.1],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","centre explicite")
  xtitle('centre explicite, cfl=0.3',' ',' ');
  xpause(100000) ;
end
//
end