%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright G. Allaire, Octobre 2000
%
% schemas numeriques sur l'equation de Burgers avec conditions aux limites periodiques
% u,t + (u^2/2),x = 0 
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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
plot(x,u,'*') ;
axis([0 lg -1.1 1.1]) ;
title('Lax-Friedrichs, cfl=0.9')

pause ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Friedrichs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:40
%
up = wshift2('1D',u,+1) ;
um = wshift2('1D',u,-1) ;
v  = (up+um)/2 - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
u = v ; 
if mod(n,2) == 0
plot(x,u,'*',x,u,':') ;
axis([0 lg -1.1 1.1]) ;
title('Lax-Friedrichs, cfl=0.9')
pause(0) ;
end
%
end

pause ;

for n=41:nt
%
up = wshift2('1D',u,+1) ;
um = wshift2('1D',u,-1) ;
v  = (up+um)/2 - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
u = v ; 
if mod(n,2) == 0
plot(x,u,'*',x,u,':') ;
axis([0 lg -1.1 1.1]) ;
title('Lax-Friedrichs, cfl=0.9')
pause(0) ;
end
%
end

pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 ;

plot(x,u,'*') ;
axis([0 lg -1.1 1.1]) ;
title('Lax-Wendroff, cfl=0.9')

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Wendroff
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:40
%
up = wshift2('1D',u,+1) ;
um = wshift2('1D',u,-1) ;
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 mod(n,2) == 0
plot(x,u,'*',x,u,':') ;
axis([0 lg -1.1 1.1]) ;
title('Lax-Wendroff, cfl=0.9')
pause(0) ;
end
%
end

pause ;

for n=41:nt
%
up = wshift2('1D',u,+1) ;
um = wshift2('1D',u,-1) ;
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 mod(n,2) == 0
plot(x,u,'*',x,u,':') ;
axis([0 lg -1.1 1.1]) ;
title('Lax-Wendroff, cfl=0.9')
pause(0) ;
end
%
end

pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 ;

plot(x,u,'*') ;
axis([0 lg -1.1 1.1]) ;
title('centre explicite, cfl=0.3')

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : centre explicite
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:120
%
up = wshift2('1D',u,+1) ;
um = wshift2('1D',u,-1) ;
v = u - (0.5*dt/dx)*(up.^2-um.^2)/2 ;
u = v ; 
if mod(n,2) == 0
plot(x,u,'*',x,u,':') ;
axis([0 lg -1.1 1.1]) ;
title('centre explicite, cfl=0.3')
pause(0) ;
end
%
end