/////////////////////////////////////////////////////////////
// Copyright G. Allaire, Octobre 2002, Janvier 2012
//
// schemas numeriques sur l'equation de transport lineaire 
// avec condition aux limites periodiques
// u,t + a u,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 ;
// a = vitesse de transport
a = 1. ;
cfl = dx/a ;
// dt = pas de temps
dt = 0.9*cfl ;
// nt = nombre de pas de temps effectues
nt = 800 ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=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 ('transport, Lax-Friedrichs, cfl=0.9' ,' ',' ') ;

halt() ;
/////////////////////////////////////////////////
// marche en temps : Lax-Friedrichs avec cfl=0.9
/////////////////////////////////////////////////
//
// w est la solution exacte
// si u est le vecteur des valeurs u(i) alors 
//   up est le vecteur des u(i+1) 
//   um est le vecteur des u(i-1)
// 
for n=1:nt
//
up = shiftp('+1',u) ;
um = shiftp('-1',u) ;
v = (up+um)/2 - (0.5*dt/dx)*a*(up-um) ;
u = v ; 
// affichage de la solution approchee et exacte
if modulo(n,5) == 0
  w = soltr(w,nx,dx,dt,a,n) ;
  clf() 
  plotframe([0,-1.1,lg,1.1],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Lax-Friedrichs")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('transport, Lax-Friedrichs, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// nouveau schema: Lax-Friedrichs avec cfl=1.1 (explose)
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 1.1*cfl ;
nt = 400 ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=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.3,lg,1.3],tics);
plot2d(x,u,1,"000")
xtitle ('transport, Lax-Friedrichs, cfl=1.1' ,' ',' ') ;

xpause(800000) ;
/////////////////////////////////////////
// marche en temps : Lax-Friedrichs
/////////////////////////////////////////
// 
for n=1:nt
//
up = shiftp('+1',u) ;
um = shiftp('-1',u) ;
v = (up+um)/2 - (0.5*dt/dx)*a*(up-um) ;
u = v ; 
if modulo(n,5) == 0
  w = soltr(w,nx,dx,dt,a,n) ;
  clf() 
  plotframe([0,-1.3,lg,1.3],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Lax-Friedrichs")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('transport, Lax-Friedrichs, cfl=1.1',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// nouveau schema: Lax-Friedrichs avec cfl=0.1 (tres dissipatif)
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 0.1*cfl ;
nt = 1000 ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=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.3,lg,1.3],tics);
plot2d(x,u,1,"000")
xtitle ('transport, Lax-Friedrichs, cfl=0.1' ,' ',' ') ;

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Lax-Friedrichs
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
up = shiftp('+1',u) ;
um = shiftp('-1',u) ;
v = (up+um)/2 - (0.5*dt/dx)*a*(up-um) ;
u = v ; 
if modulo(n,5) == 0
  w = soltr(w,nx,dx,dt,a,n) ;
  clf() 
  plotframe([0,-1.3,lg,1.3],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Lax-Friedrichs")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('transport, Lax-Friedrichs, cfl=0.1',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// nouveau schema: centre explicite avec cfl=0.5 (explose)
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 0.5*cfl ;
nt = 400 ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=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.3,lg,1.3],tics);
plot2d(x,u,1,"000")
xtitle ('transport, centré explicite, cfl=0.5' ,' ',' ') ;

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : centre explicite 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
up = shiftp('+1',u) ;
um = shiftp('-1',u) ;
v = u - (0.5*dt/dx)*a*(up-um) ;
u = v ; 
if modulo(n,5) == 0
  w = soltr(w,nx,dx,dt,a,n) ;
  clf() 
  plotframe([0,-1.3,lg,1.3],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","centré explicite")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('transport, centré explicite, cfl=0.5',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// nouveau schema: centre explicite avec cfl=0.1 (n'explose pas)
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 0.1*cfl ;
nt = 1000 ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=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.3,lg,1.3],tics);
plot2d(x,u,1,"000")
xtitle ('transport, centré explicite, cfl=0.1' ,' ',' ') ;

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : centre explicite 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
up = shiftp('+1',u) ;
um = shiftp('-1',u) ;
v = u - (0.5*dt/dx)*a*(up-um) ;
u = v ; 
if modulo(n,5) == 0
  w = soltr(w,nx,dx,dt,a,n) ;
  clf() 
  plotframe([0,-1.3,lg,1.3],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","centré explicite")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('transport, centré explicite, cfl=0.1',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// nouveau schema: centre implicite avec cfl=0.9
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 0.9*cfl ;
nt = 400 ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  u(i) = sin(x(i)*2*pi) ;
end
w = u ;
mat=zeros(nx,nx) ;
for i=2:nx-1
  mat(i,i) = 1. ;
  mat(i,i+1) = (0.5*dt/dx)*a ;
  mat(i,i-1) = -(0.5*dt/dx)*a ;
end
mat(1,1) = 1. ;
mat(1,2) = (0.5*dt/dx)*a ;
mat(1,nx) = -(0.5*dt/dx)*a ;
mat(nx,nx-1)= -(0.5*dt/dx)*a ;
mat(nx,1) = (0.5*dt/dx)*a ;
mat(nx,nx) = 1. ;
smat = sparse(mat) ;

clf()
tics=[4,10,4,10];
plotframe([0,-1.3,lg,1.3],tics);
plot2d(x,u,1,"000")
xtitle ('transport, centré implicite, cfl=0.9' ,' ',' ') ;

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : centre implicite 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
u0=u';
v=smat\u0 ;
u = v' ; 
if modulo(n,5) == 0
  w = soltr(w,nx,dx,dt,a,n) ;
  clf() 
  plotframe([0,-1.3,lg,1.3],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","centré implicite")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('transport, centré implicite, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// nouveau schema: centre implicite avec cfl=2.0
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 2.0*cfl ;
nt = 200 ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  u(i) = sin(x(i)*2*pi) ;
end
w = u ;
mat=zeros(nx,nx) ;
for i=2:nx-1
  mat(i,i) = 1. ;
  mat(i,i+1) = (0.5*dt/dx)*a ;
  mat(i,i-1) = -(0.5*dt/dx)*a ;
end
mat(1,1) = 1. ;
mat(1,2) = (0.5*dt/dx)*a ;
mat(1,nx) = -(0.5*dt/dx)*a ;
mat(nx,nx-1)= -(0.5*dt/dx)*a ;
mat(nx,1) = (0.5*dt/dx)*a ;
mat(nx,nx) = 1. ;
smat = sparse(mat) ;

clf()
tics=[4,10,4,10];
plotframe([0,-1.3,lg,1.3],tics);
plot2d(x,u,1,"000")
xtitle ('transport, centré implicite, cfl=2.0' ,' ',' ') ;

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : centre implicite 
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
u0=u';
v=smat\u0 ;
u = v' ; 
if modulo(n,5) == 0
  w = soltr(w,nx,dx,dt,a,n) ;
  clf() 
  plotframe([0,-1.3,lg,1.3],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","centré implicite")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('transport, centré implicite, cfl=2.0',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// nouveau schema: Lax-Wendroff avec cfl=0.9
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 0.9*cfl ;
nt = 800 ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=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.3,lg,1.3],tics);
plot2d(x,u,1,"000")
xtitle ('transport, Lax-Wendroff, cfl=0.9' ,' ',' ') ;

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Lax-Wendroff
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
up = shiftp('+1',u) ;
um = shiftp('-1',u) ;
v = u - (0.5*dt/dx)*a*(up-um) + ... 
(0.5*dt*dt/(dx*dx))*( a*a*(up-u) ...
        - a*a*(u-um) ) ;
u = v ; 
if modulo(n,5) == 0
  w = soltr(w,nx,dx,dt,a,n) ;
  clf() 
  plotframe([0,-1.3,lg,1.3],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Lax-Wendroff")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('transport, Lax-Wendroff, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// donnee initiale = creneau
// nouveau schema: Lax-Friedrichs avec cfl=0.9
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 0.9*cfl ;
nt = 400 ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
  u(i) = -1. ;
  else
  u(i) = 1. ;
  end
end
w = u ;

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

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Lax-Friedrichs
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
up = shiftp('+1',u) ;
um = shiftp('-1',u) ;
v = (up+um)/2 - (0.5*dt/dx)*a*(up-um) ;
u = v ; 
if modulo(n,5) == 0
  for i=1:nx
      xi = (i-1)*dx-a*n*dt  ;
      xint = floor(xi/lg) ;
      xi = xi - xint*lg ; 
      if xi < lg/2
         w(i) = -1. ;
      else
         w(i) = 1. ;
      end
  end
  clf() 
  plotframe([0,-1.3,lg,1.3],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Lax-Friedrichs")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('transport, Lax-Friedrichs, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end

halt() ;

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// donnee initiale = creneau
// nouveau schema: Lax-Wendroff avec cfl=0.9
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//
dt = 0.9*cfl ;
nt = 400 ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
w=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  if x(i) < lg/2
  u(i) = -1. ;
  else
  u(i) = 1. ;
  end
end
w = u ;

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

xpause(800000) ;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// marche en temps : Lax-Wendroff
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
up = shiftp('+1',u) ;
um = shiftp('-1',u) ;
v = u - (0.5*dt/dx)*a*(up-um) + ... 
(0.5*dt*dt/(dx*dx))*( a*a*(up-u) ...
        - a*a*(u-um) ) ;
u = v ; 
if modulo(n,5) == 0
  for i=1:nx
      xi = (i-1)*dx-a*n*dt  ;
      xint = floor(xi/lg) ;
      xi = xi - xint*lg ; 
      if xi < lg/2
         w(i) = -1. ;
      else
         w(i) = 1. ;
      end
  end
  clf() 
  plotframe([0,-1.3,lg,1.3],tics);
  xset("thickness",2) ;
  plot2d(x,u,[1,1],"100","Lax-Wendroff")
  xset("thickness",1) ;
  plot2d(x,w,[2,4],"100","exact")
  xtitle('transport, Lax-Wendroff, cfl=0.9',' ',' ');
  xpause(100000) ;
end
//
end