%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright G. Allaire, Octobre 2000
%
% schemas numeriques sur 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
a = 1. ;
if a==0 
 'pas interessant !'
 return
else
 cfl = dx/a ;
end
% 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
plot(x,u) ;
axis([0 lg -1.1 1.1]) ;
title('Lax-Friedrichs, cfl=0.9')

pause ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 = wshift2('1D',u,+1) ;
um = wshift2('1D',u,-1) ;
v = (up+um)/2 - (0.5*dt/dx)*a*(up-um) ;
u = v ; 
% affichage de la solution approchee et exacte
if mod(n,5) == 0
  w = soltr(w,nx,dx,dt,a,n) ;
  plot(x,u,'o',x,w) ;
  axis([0 lg -1.1 1.1]) ;
  title('Lax-Friedrichs, cfl=0.9')
  pause(0) ;
end
%
end

pause ;

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

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

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

pause ;

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

plot(x,u) ;
axis([0 lg -1.1 1.1]) ;
title('Lax-Friedrichs, cfl=0.1')

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Friedrichs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
up = wshift2('1D',u,+1) ;
um = wshift2('1D',u,-1) ;
v = (up+um)/2 - (0.5*dt/dx)*a*(up-um) ;
u = v ; 
if mod(n,5) == 0
w = soltr(w,nx,dx,dt,a,n) ;
plot(x,u,'o',x,w) ;
axis([0 lg -1.1 1.1]) ;
title('Lax-Friedrichs, cfl=0.1')
pause(0) ;
end
%
end

pause ;


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

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

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : centre explicite 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
up = wshift2('1D',u,+1) ;
um = wshift2('1D',u,-1) ;
v = u - (0.5*dt/dx)*a*(up-um) ;
u = v ; 
if mod(n,5) == 0
w = soltr(w,nx,dx,dt,a,n) ;
plot(x,u,'o',x,w) ;
axis([0 lg -1.1 1.1]) ;
title('centre explicite, cfl=0.5')
pause(0) ;
end
%
end

pause ;

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

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

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : centre explicite 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
up = wshift2('1D',u,+1) ;
um = wshift2('1D',u,-1) ;
v = u - (0.5*dt/dx)*a*(up-um) ;
u = v ; 
if mod(n,5) == 0
w = soltr(w,nx,dx,dt,a,n) ;
plot(x,u,'o',x,w) ;
axis([0 lg -1.1 1.1]) ;
title('centre explicite, cfl=0.1')
pause(0) ;
end
%
end

pause ;


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

plot(x,u) ;
axis([0 lg -1.1 1.1]) ;
title('centre implicite, cfl=0.9')

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : centre implicite 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
u0=u';
v=smat\u0 ;
u = v' ; 
if mod(n,5) == 0
w = soltr(w,nx,dx,dt,a,n) ;
plot(x,u,'o',x,w) ;
axis([0 lg -1.1 1.1]) ;
title('centre implicite, cfl=0.9')
pause(0) ;
end
%
end

pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% nouveau schema: centre implicite avec cfl=2.0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
dt = 2.0*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) ;

plot(x,u) ;
axis([0 lg -1.1 1.1]) ;
title('centre implicite, cfl=2.0')

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : centre implicite 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
u0=u';
v=smat\u0 ;
u = v' ; 
if mod(n,5) == 0
w = soltr(w,nx,dx,dt,a,n) ;
plot(x,u,'o',x,w) ;
axis([0 lg -1.1 1.1]) ;
title('centre implicite, cfl=2.0')
pause(0) ;
end
%
end

pause ;


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

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:nt
%
up = wshift2('1D',u,+1) ;
um = wshift2('1D',u,-1) ;
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 mod(n,5) == 0
w = soltr(w,nx,dx,dt,a,n) ;
plot(x,u,'o',x,w) ;
axis([0 lg -1.1 1.1]) ;
title('Lax-Wendroff, cfl=0.9')
pause(0) ;
end
%
end

pause ;

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

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:nt
%
up = wshift2('1D',u,+1) ;
um = wshift2('1D',u,-1) ;
v = (up+um)/2 - (0.5*dt/dx)*a*(up-um) ;
u = v ; 
if mod(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
plot(x,u,'o',x,w) ;
axis([0 lg -1.1 1.1]) ;
title('Lax-Friedrichs, cfl=0.9')
pause(0) ;
end
%
end

pause ;



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

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

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


