%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright G. Allaire, Novembre 2000
%
% schemas numeriques a limiteur de flux (TVD d'ordre 2) 
% pour l'equation de transport lineaire 
% avec condition aux limites periodiques
% u,t + a u,x = 0
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% lg = longueur du domaine
lg = 1. ;
% nx = nombre de mailles
nx = 100 ;
% dx = pas d'espace
dx = lg/nx ;
% a = vitesse de transport (doit etre positive)
a = 1. ;
% dt = pas de temps
cfl = dx/a;
dt = 0.9*cfl ;
% nt = nombre de pas de temps pour faire un tour
nt = floor(lg/(a*dt)) ;
%
% initialisation
%
x=zeros(1,nx) ;
% uf = solution de Lax-Friedrichs
uf=zeros(1,nx) ;
% uw = solution de Lax-Wendroff
uw=zeros(1,nx) ;
% ua = solution upwind (amont)
ua=zeros(1,nx) ;
% ue = solution exacte (donnee initiale)
ue=zeros(1,nx) ;
v=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
arg=zeros(1,nx) ;
%
for i=1:nx
  x(i) = (i-1)*dx  ;
  ue(i) = 0. ;
  if abs(x(i)-0.3)<=0.1 
     ue(i) = 1. ;
  end
end
uf = ue ;
uw = ue ;
ua = ue ;

% affichage de la donnee initiale
plot(x,ue) ;
axis([0 lg -0.2 1.2]) ;
title('Lax-Friedrichs, Lax-Wendroff, upwind, cfl=0.9')

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Friedrichs, Lax-Wendroff, upwind
% cfl=0.9 comparaison sur des temps petits
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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:3*nt
%
% Lax-Friedrichs
up = wshift2('1D',uf,+1) ;
um = wshift2('1D',uf,-1) ;
v = (up+um)/2 - (0.5*dt/dx)*a*(up-um) ;
uf = v ; 
% Lax-Wendroff
up = wshift2('1D',uw,+1) ;
um = wshift2('1D',uw,-1) ;
v = uw - (0.5*dt/dx)*a*(up-um) + ... 
(0.5*dt*dt/(dx*dx))*( a*a*(up-uw) - a*a*(uw-um) ) ;
uw = v ; 
% upwind
um = wshift2('1D',ua,-1) ;
v = ua - (dt/dx)*a*(ua-um) ;
ua = v ; 
% affichage de la solution approchee et exacte
if mod(n,5) == 0
%
for i=1:nx
  ue(i) = 0. ;
  arg(i) = mod(x(i)-n*a*dt,lg) ;
  if abs(arg(i)-0.3)<=0.1 
     ue(i) = 1. ;
  end
end

  plot(x,ue,x,uf,':',x,uw,'d',x,ua,'o') ;
  axis([0 lg -0.2 1.2]) ;
  title('Lax-Friedrichs, Lax-Wendroff, upwind, cfl=0.9')
  legend('exact','Lax-Friedrichs','Lax-Wendroff','upwind') ;
  pause(0) ;
end
%
end

pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lax-Friedrichs, Lax-Wendroff, upwind
% cas tres particulier cfl=0.999
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% dt = pas de temps
dt = 0.999*cfl ;
% nt = nombre de pas de temps pour faire un tour
nt = floor(lg/(a*dt)) ;
%
% initialisation
%
x=zeros(1,nx) ;
% uf = solution de Lax-Friedrichs
uf=zeros(1,nx) ;
% uw = solution de Lax-Wendroff
uw=zeros(1,nx) ;
% ua = solution upwind (amont)
ua=zeros(1,nx) ;
% ue = solution exacte (donnee initiale)
ue=zeros(1,nx) ;
v=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
arg=zeros(1,nx) ;
%
for i=1:nx
  x(i) = (i-1)*dx  ;
  ue(i) = 0. ;
  if abs(x(i)-0.3)<=0.1 
     ue(i) = 1. ;
  end
end
uf = ue ;
uw = ue ;
ua = ue ;

% affichage de la donnee initiale
plot(x,ue) ;
axis([0 lg -0.2 1.2]) ;
title('Lax-Friedrichs, Lax-Wendroff, upwind, cfl=0.999')

pause(0) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Friedrichs, Lax-Wendroff, upwind
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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:3*nt
%
% Lax-Friedrichs
up = wshift2('1D',uf,+1) ;
um = wshift2('1D',uf,-1) ;
v = (up+um)/2 - (0.5*dt/dx)*a*(up-um) ;
uf = v ; 
% Lax-Wendroff
up = wshift2('1D',uw,+1) ;
um = wshift2('1D',uw,-1) ;
v = uw - (0.5*dt/dx)*a*(up-um) + ... 
(0.5*dt*dt/(dx*dx))*( a*a*(up-uw) - a*a*(uw-um) ) ;
uw = v ; 
% upwind
um = wshift2('1D',ua,-1) ;
v = ua - (dt/dx)*a*(ua-um) ;
ua = v ; 
% affichage de la solution approchee et exacte
if mod(n,5) == 0
%
for i=1:nx
  ue(i) = 0. ;
  arg(i) = mod(x(i)-n*a*dt,lg) ;
  if abs(arg(i)-0.3)<=0.1 
     ue(i) = 1. ;
  end
end

  plot(x,ue,x,uf,':',x,uw,'d',x,ua,'o') ;
  axis([0 lg -0.2 1.2]) ;
  title('Lax-Friedrichs, Lax-Wendroff, upwind, cfl=0.999')
  legend('exact','Lax-Friedrichs','Lax-Wendroff','upwind') ;
  pause(0) ;
end
%
end


pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lax-Friedrichs, Lax-Wendroff, upwind
% cfl=0.999 comparaison sur 50 tours
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% dt = pas de temps
dt = 0.999*cfl ;
% nt = nombre de pas de temps pour faire un tour
nt = floor(lg/(a*dt)) ;
%
% initialisation
%
x=zeros(1,nx) ;
% uf = solution de Lax-Friedrichs
uf=zeros(1,nx) ;
% uw = solution de Lax-Wendroff
uw=zeros(1,nx) ;
% ua = solution upwind (amont)
ua=zeros(1,nx) ;
% ue = solution exacte (donnee initiale)
ue=zeros(1,nx) ;
v=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
arg=zeros(1,nx) ;
%
for i=1:nx
  x(i) = (i-1)*dx  ;
  ue(i) = 0. ;
  if abs(x(i)-0.3)<=0.1 
     ue(i) = 1. ;
  end
end
uf = ue ;
uw = ue ;
ua = ue ;

% affichage de la donnee initiale
plot(x,ue) ;
axis([0 lg -0.2 1.2]) ;
title('Lax-Friedrichs, Lax-Wendroff, upwind, cfl=0.999')

pause(0) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Friedrichs, Lax-Wendroff, upwind
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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:50*nt
%
% Lax-Friedrichs
up = wshift2('1D',uf,+1) ;
um = wshift2('1D',uf,-1) ;
v = (up+um)/2 - (0.5*dt/dx)*a*(up-um) ;
uf = v ; 
% Lax-Wendroff
up = wshift2('1D',uw,+1) ;
um = wshift2('1D',uw,-1) ;
v = uw - (0.5*dt/dx)*a*(up-um) + ... 
(0.5*dt*dt/(dx*dx))*( a*a*(up-uw) - a*a*(uw-um) ) ;
uw = v ; 
% upwind
um = wshift2('1D',ua,-1) ;
v = ua - (dt/dx)*a*(ua-um) ;
ua = v ; 
% affichage de la solution approchee et exacte
if mod(n,2*nt) == 0
%
for i=1:nx
  ue(i) = 0. ;
  arg(i) = mod(x(i)-n*a*dt,lg) ;
  if abs(arg(i)-0.3)<=0.1 
     ue(i) = 1. ;
  end
end

  plot(x,ue,x,uf,':',x,uw,'d',x,ua,'o') ;
  axis([0 lg -0.2 1.2]) ;
  in = floor(n/nt) ;
  title('50 tours, cfl=0.999')
  legend('exact','Lax-Friedrichs','Lax-Wendroff','upwind') ;
  pause(0) ;
end
%
end


pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% schemas TVD
% cfl=0.9 comparaison sur 50 tours
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% dt = pas de temps
dt = 0.5*cfl ;
% nt = nombre de pas de temps pour faire un tour
nt = floor(lg/(a*dt)) ;
%
% initialisation
%
x=zeros(1,nx) ;
% phi = limiteur
phi=zeros(1,nx) ;
phim=zeros(1,nx) ;
% du = increment u(i+1)-u(i)
du=zeros(1,nx) ;
% ua = solution upwind (amont)
ua=zeros(1,nx) ;
% ur = solution TVD minmod
ur=zeros(1,nx) ;
% us = solution TVD superbee
us=zeros(1,nx) ;
% ub = solution TVD ultrabee
ub=zeros(1,nx) ;
% uk = solution ordre 2 espace, ordre 1 temps
uk=zeros(1,nx) ;
% ue = solution exacte (donnee initiale)
ue=zeros(1,nx) ;
v=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
f=zeros(1,nx) ;
fm=zeros(1,nx) ;
mi=zeros(1,nx) ;
ma=zeros(1,nx) ;
arg=zeros(1,nx) ;
%
for i=1:nx
  x(i) = (i-1)*dx  ;
  ue(i) = 0. ;
  if abs(x(i)-0.3)<=0.1 
     ue(i) = 1. ;
  end
end
ua = ue ;
ur = ue ;
us = ue ;
ub = ue ;
uk = ue ;

% affichage de la donnee initiale
plot(x,ue) ;
axis([0 lg -0.2 1.2]) ;
title('schemas TVD, cfl=0.9')

pause(0) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : TVD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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:50*nt
%
% upwind
um = wshift2('1D',ua,-1) ;
v = ua - (dt/dx)*a*(ua-um) ;
ua = v ; 
% minod
um = wshift2('1D',ur,-1) ;
up = wshift2('1D',ur,+1) ;
du = ur-um ;
du = du + 1.e-6*(1+sign(du)).*(1-sign(du)) ;
dum = wshift2('1D',du,-1) ;
phi = dum./du ;
phi = max(0.,min(1.,phi)) ;
phim = wshift2('1D',phi,-1) ;
v = ur - (dt/dx)*a*(ur-um) + 0.5*(dt/dx)*a*...
    ((dt/dx)*a-1.)*(phi.*(up-ur)-phim.*(ur-um)) ;
ur = v ; 
% superbee
um = wshift2('1D',us,-1) ;
up = wshift2('1D',us,+1) ;
du = us-um ;
du = du + 1.e-6*(1+sign(du)).*(1-sign(du)) ;
dum = wshift2('1D',du,-1) ;
phi = dum./du ;
phi = max(min(2*phi,1.),min(2.,phi)) ;
phi = max(phi,0.) ;
phim = wshift2('1D',phi,-1) ;
v = us - (dt/dx)*a*(us-um) + 0.5*(dt/dx)*a*...
    ((dt/dx)*a-1.)*(phi.*(up-us)-phim.*(us-um)) ;
us = v ; 
% ultrabee a la Despres-Lagoutiere
um = wshift2('1D',ub,-1) ;
up = wshift2('1D',ub,+1) ;
mi = min(ub,um) ;
ma = max(ub,um) ;
for i=1:nx
  binf = ma(i) + (ub(i)-ma(i))/((dt/dx)*a) ;
  bsup = mi(i) + (ub(i)-mi(i))/((dt/dx)*a) ;
  if up(i) <= binf 
     f(i) = binf ;
  elseif up(i) >= bsup
     f(i) = bsup ;
  else
     f(i) = up(i) ;
  end
end
fm = wshift2('1D',f,-1) ;
v = ub - (dt/dx)*a*(f-fm) ;
ub = v ;
% ordre 2 en x et 1 en t
um = wshift2('1D',uk,-1) ;
up = wshift2('1D',uk,+1) ;
du = up-uk ;
dum = wshift2('1D',du,-1) ;
phi = sign(du).*max(0.,min(abs(du),sign(du).*dum)) ;
f = uk + 0.5*phi ;
fm = wshift2('1D',f,-1) ;
v = uk - (dt/dx)*a*(f-fm) ;
uk = v ; 

% affichage de la solution approchee et exacte
if mod(n,nt) == 0
%
for i=1:nx
  ue(i) = 0. ;
  arg(i) = mod(x(i)-n*a*dt,lg) ;
  if abs(arg(i)-0.3)<=0.1 
     ue(i) = 1. ;
  end
end

  plot(x,ue,x,ua,':',x,ur,'d',x,us,'o',x,ub,'<',x,uk,'p') ;
  axis([0 lg -0.2 1.2]) ;
  in = floor(n/nt) ;
  title('creneau TVD 50 tours, cfl=0.9')
  legend('exact','upwind','minmod','superbee','ultrabee DL','ordre 2-1') ;
  pause(0) ;
end
%
end


pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% schemas TVD donnee reguliere
% cfl=0.9 comparaison sur 200 tours
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% dt = pas de temps
dt = 0.9*cfl ;
% nt = nombre de pas de temps pour faire un tour
nt = floor(lg/(a*dt)) ;
%
% initialisation
%
x=zeros(1,nx) ;
% phi = limiteur
phi=zeros(1,nx) ;
phim=zeros(1,nx) ;
% du = increment u(i+1)-u(i)
du=zeros(1,nx) ;
% ua = solution upwind (amont)
ua=zeros(1,nx) ;
% ur = solution TVD minmod
ur=zeros(1,nx) ;
% us = solution TVD superbee
us=zeros(1,nx) ;
% ub = solution TVD ultrabee
ub=zeros(1,nx) ;
% ue = solution exacte (donnee initiale)
ue=zeros(1,nx) ;
v=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
f=zeros(1,nx) ;
fm=zeros(1,nx) ;
mi=zeros(1,nx) ;
ma=zeros(1,nx) ;
arg=zeros(1,nx) ;
%
for i=1:nx
  x(i) = (i-1)*dx  ;
  ue(i)= sin(x(i)*2*pi) ;
end
ua = ue ;
ur = ue ;
us = ue ;
ub = ue ;

% affichage de la donnee initiale
plot(x,ue) ;
axis([0 lg -1.2 1.2]) ;
title('schemas TVD, cfl=0.9')

pause(0) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : TVD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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:200*nt
%
% upwind
um = wshift2('1D',ua,-1) ;
v = ua - (dt/dx)*a*(ua-um) ;
ua = v ; 
% minod
um = wshift2('1D',ur,-1) ;
up = wshift2('1D',ur,+1) ;
du = ur-um ;
du = du + 1.e-6*(1+sign(du)).*(1-sign(du)) ;
dum = wshift2('1D',du,-1) ;
phi = dum./du ;
phi = max(0.,min(1.,phi)) ;
phim = wshift2('1D',phi,-1) ;
v = ur - (dt/dx)*a*(ur-um) + 0.5*(dt/dx)*a*...
    ((dt/dx)*a-1.)*(phi.*(up-ur)-phim.*(ur-um)) ;
ur = v ; 
% superbee
um = wshift2('1D',us,-1) ;
up = wshift2('1D',us,+1) ;
du = us-um ;
du = du + 1.e-6*(1+sign(du)).*(1-sign(du)) ;
dum = wshift2('1D',du,-1) ;
phi = dum./du ;
phi = max(min(2*phi,1.),min(2.,phi)) ;
phi = max(phi,0.) ;
phim = wshift2('1D',phi,-1) ;
v = us - (dt/dx)*a*(us-um) + 0.5*(dt/dx)*a*...
    ((dt/dx)*a-1.)*(phi.*(up-us)-phim.*(us-um)) ;
us = v ; 
% ultrabee a la Despres-Lagoutiere
um = wshift2('1D',ub,-1) ;
up = wshift2('1D',ub,+1) ;
mi = min(ub,um) ;
ma = max(ub,um) ;
for i=1:nx
  binf = ma(i) + (ub(i)-ma(i))/((dt/dx)*a) ;
  bsup = mi(i) + (ub(i)-mi(i))/((dt/dx)*a) ;
  if up(i) <= binf 
     f(i) = binf ;
  elseif up(i) >= bsup
     f(i) = bsup ;
  else
     f(i) = up(i) ;
  end
end
fm = wshift2('1D',f,-1) ;
v = ub - (dt/dx)*a*(f-fm) ;
ub = v ;

% affichage de la solution approchee et exacte
if mod(n,4*nt) == 0
%
for i=1:nx
  ue(i) = 0. ;
  arg(i) = mod(x(i)-n*a*dt,lg) ;
  ue(i) = sin(arg(i)*2*pi) ;
end

  plot(x,ue,x,ua,':',x,ur,'d',x,us,'o',x,ub,'<') ;
  axis([0 lg -1.2 1.2]) ;
  in = floor(n/nt) ;
  title('sinus TVD 200 tours, cfl=0.9')
  legend('exact','upwind','minmod','superbee','ultrabee DL') ;
  pause(0) ;
end
%
end


pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% schemas TVD donnee initiale reguliere
% cfl=0.66 comparaison sur 20 tours
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% dt = pas de temps
dt = 0.66*cfl ;
% nt = nombre de pas de temps pour faire un tour
nt = floor(lg/(a*dt)) ;
%
% initialisation
%
x=zeros(1,nx) ;
% phi = limiteur
phi=zeros(1,nx) ;
phim=zeros(1,nx) ;
% du = increment u(i+1)-u(i)
du=zeros(1,nx) ;
% ua = solution upwind (amont)
ua=zeros(1,nx) ;
% ur = solution TVD minmod
ur=zeros(1,nx) ;
% us = solution TVD superbee
us=zeros(1,nx) ;
% ub = solution TVD ultrabee
ub=zeros(1,nx) ;
% uk = solution ordre 2 espace, ordre 1 temps
uk=zeros(1,nx) ;
% ue = solution exacte (donnee initiale)
ue=zeros(1,nx) ;
v=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
f=zeros(1,nx) ;
fm=zeros(1,nx) ;
mi=zeros(1,nx) ;
ma=zeros(1,nx) ;
arg=zeros(1,nx) ;
%
for i=1:nx
  x(i) = (i-1)*dx  ;
  ue(i) = 0.5*(1.+sin(2*pi*x(i))) ;
%  ue(i) = exp(-30.*(0.5-x(i))^2) ;
end
ua = ue ;
ur = ue ;
us = ue ;
ub = ue ;
uk = ue ;

% affichage de la donnee initiale
plot(x,ue) ;
axis([0 lg -0.2 1.2]) ;
title('schemas TVD, cfl=0.66')

pause(0) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : TVD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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:20*nt
%
% upwind
um = wshift2('1D',ua,-1) ;
v = ua - (dt/dx)*a*(ua-um) ;
ua = v ; 
% minod
um = wshift2('1D',ur,-1) ;
up = wshift2('1D',ur,+1) ;
du = ur-um ;
du = du + 1.e-6*(1+sign(du)).*(1-sign(du)) ;
dum = wshift2('1D',du,-1) ;
phi = dum./du ;
phi = max(0.,min(1.,phi)) ;
phim = wshift2('1D',phi,-1) ;
v = ur - (dt/dx)*a*(ur-um) + 0.5*(dt/dx)*a*...
    ((dt/dx)*a-1.)*(phi.*(up-ur)-phim.*(ur-um)) ;
ur = v ; 
% superbee
um = wshift2('1D',us,-1) ;
up = wshift2('1D',us,+1) ;
du = us-um ;
du = du + 1.e-6*(1+sign(du)).*(1-sign(du)) ;
dum = wshift2('1D',du,-1) ;
phi = dum./du ;
phi = max(min(2*phi,1.),min(2.,phi)) ;
phi = max(phi,0.) ;
phim = wshift2('1D',phi,-1) ;
v = us - (dt/dx)*a*(us-um) + 0.5*(dt/dx)*a*...
    ((dt/dx)*a-1.)*(phi.*(up-us)-phim.*(us-um)) ;
us = v ; 
% ultrabee a la Despres-Lagoutiere
um = wshift2('1D',ub,-1) ;
up = wshift2('1D',ub,+1) ;
mi = min(ub,um) ;
ma = max(ub,um) ;
for i=1:nx
  binf = ma(i) + (ub(i)-ma(i))/((dt/dx)*a) ;
  bsup = mi(i) + (ub(i)-mi(i))/((dt/dx)*a) ;
  if up(i) <= binf 
     f(i) = binf ;
  elseif up(i) >= bsup
     f(i) = bsup ;
  else
     f(i) = up(i) ;
  end
end
fm = wshift2('1D',f,-1) ;
v = ub - (dt/dx)*a*(f-fm) ;
ub = v ;
% ordre 2 en x et 1 en t
um = wshift2('1D',uk,-1) ;
up = wshift2('1D',uk,+1) ;
du = up-uk ;
dum = wshift2('1D',du,-1) ;
phi = sign(du).*max(0.,min(abs(du),sign(du).*dum)) ;
f = uk + 0.5*phi ;
fm = wshift2('1D',f,-1) ;
v = uk - (dt/dx)*a*(f-fm) ;
uk = v ; 

% affichage de la solution approchee et exacte
if mod(n,20) == 0
%
  plot(x,ua,':',x,ur,'d',x,us,'o',x,ub,'<',x,uk,'p') ;
  axis([0 lg -0.2 1.2]) ;
  in = floor(n/nt) ;
  title('creneau TVD 50 tours, cfl=0.66')
  legend('upwind','minmod','superbee','ultrabee DL','ordre 2x-1t') ;
  pause(0) ;
end
%
end


pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% schemas TVD donnee initiale reguliere
% cfl=0.5 comparaison sur 20 tours
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% dt = pas de temps
dt = 0.5*cfl ;
% nt = nombre de pas de temps pour faire un tour
nt = floor(lg/(a*dt)) ;
%
% initialisation
%
x=zeros(1,nx) ;
% phi = limiteur
phi=zeros(1,nx) ;
phim=zeros(1,nx) ;
% du = increment u(i+1)-u(i)
du=zeros(1,nx) ;
% ua = solution upwind (amont)
ua=zeros(1,nx) ;
% ur = solution TVD minmod
ur=zeros(1,nx) ;
% us = solution TVD superbee
us=zeros(1,nx) ;
% ub = solution TVD ultrabee
ub=zeros(1,nx) ;
% uk = solution ordre 2 espace, ordre 1 temps
uk=zeros(1,nx) ;
% ue = solution exacte (donnee initiale)
ue=zeros(1,nx) ;
v=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
f=zeros(1,nx) ;
fm=zeros(1,nx) ;
mi=zeros(1,nx) ;
ma=zeros(1,nx) ;
arg=zeros(1,nx) ;
%
for i=1:nx
  x(i) = (i-1)*dx  ;
  ue(i) = 0.5*(1.+sin(2*pi*x(i))) ;
%  ue(i) = exp(-30.*(0.5-x(i))^2) ;
end
ua = ue ;
ur = ue ;
us = ue ;
ub = ue ;
uk = ue ;

% affichage de la donnee initiale
plot(x,ue) ;
axis([0 lg -0.2 1.2]) ;
title('schemas TVD, cfl=0.66')

pause(0) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : TVD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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:50*nt
%
% upwind
um = wshift2('1D',ua,-1) ;
v = ua - (dt/dx)*a*(ua-um) ;
ua = v ; 
% minod
um = wshift2('1D',ur,-1) ;
up = wshift2('1D',ur,+1) ;
du = ur-um ;
du = du + 1.e-6*(1+sign(du)).*(1-sign(du)) ;
dum = wshift2('1D',du,-1) ;
phi = dum./du ;
phi = max(0.,min(1.,phi)) ;
phim = wshift2('1D',phi,-1) ;
v = ur - (dt/dx)*a*(ur-um) + 0.5*(dt/dx)*a*...
    ((dt/dx)*a-1.)*(phi.*(up-ur)-phim.*(ur-um)) ;
ur = v ; 
% superbee
um = wshift2('1D',us,-1) ;
up = wshift2('1D',us,+1) ;
du = us-um ;
du = du + 1.e-6*(1+sign(du)).*(1-sign(du)) ;
dum = wshift2('1D',du,-1) ;
phi = dum./du ;
phi = max(min(2*phi,1.),min(2.,phi)) ;
phi = max(phi,0.) ;
phim = wshift2('1D',phi,-1) ;
v = us - (dt/dx)*a*(us-um) + 0.5*(dt/dx)*a*...
    ((dt/dx)*a-1.)*(phi.*(up-us)-phim.*(us-um)) ;
us = v ; 
% ultrabee a la Despres-Lagoutiere
um = wshift2('1D',ub,-1) ;
up = wshift2('1D',ub,+1) ;
mi = min(ub,um) ;
ma = max(ub,um) ;
for i=1:nx
  binf = ma(i) + (ub(i)-ma(i))/((dt/dx)*a) ;
  bsup = mi(i) + (ub(i)-mi(i))/((dt/dx)*a) ;
  if up(i) <= binf 
     f(i) = binf ;
  elseif up(i) >= bsup
     f(i) = bsup ;
  else
     f(i) = up(i) ;
  end
end
fm = wshift2('1D',f,-1) ;
v = ub - (dt/dx)*a*(f-fm) ;
ub = v ;
% ordre 2 en x et 1 en t
um = wshift2('1D',uk,-1) ;
up = wshift2('1D',uk,+1) ;
du = up-uk ;
dum = wshift2('1D',du,-1) ;
phi = sign(du).*max(0.,min(abs(du),sign(du).*dum)) ;
f = uk + 0.5*phi ;
fm = wshift2('1D',f,-1) ;
v = uk - (dt/dx)*a*(f-fm) ;
uk = v ; 

% affichage de la solution approchee et exacte
if mod(n,20) == 0
%
  plot(x,ua,':',x,ur,'d',x,us,'o',x,ub,'<',x,uk,'p') ;
  axis([0 lg -0.2 1.2]) ;
  in = floor(n/nt) ;
  title('creneau TVD 50 tours, cfl=0.5')
  legend('upwind','minmod','superbee','ultrabee DL','ordre 2x-1t') ;
  pause(0) ;
end
%
end