%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright G. Allaire, Octobre 2000
%
% schemas numeriques sur un systeme lineaire hyperbolique ou non
% comparaison equation des ondes - systeme de Cauchy-Riemann 
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% equation des ondes
% u,t + v,x = 0
% v,t + u,x = 0
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% lg = longueur du domaine
lg = 1. ;
% nx = nombre de mailles
nx = 500 ;
% dx = pas d'espace
dx = lg/nx ;
cfl = dx ;
% dt = pas de temps
dt = 0.9*cfl ;
% nt = nombre de pas de temps effectues
nt = 1000 ;
%
% initialisation
%
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
us=zeros(1,nx) ;
vs=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
vp=zeros(1,nx) ;
vm=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  u(i) = sin(x(i)*2*pi) ;
  v(i) = 0. ;
end

%M = moviein(1);

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

%M(:,1)=getframe;
%mm=2;

pause ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Friedrichs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
up = wshift2('1D',u,+1) ;
um = wshift2('1D',u,-1) ;
vp = wshift2('1D',v,+1) ;
vm = wshift2('1D',v,-1) ;
vs = (vp+vm)/2 - (0.5*dt/dx)*(up-um) ;
us = (up+um)/2 - (0.5*dt/dx)*(vp-vm) ;
v = vs ;
u = us ;
if mod(n,5) == 0
plot(x,u) ;
axis([0 lg -1.1 1.1]) ;
title('Lax-Friedrichs, cfl=0.9')
%M(mm)=getframe;
%mm=mm+1;
pause(0) ;
end
%
end

pause ;
%newplot;
%movie(M,3);
%pause;
% movie2avi(M,'animhyperb') ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% initialisation
%
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
us=zeros(1,nx) ;
vs=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
vp=zeros(1,nx) ;
vm=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  u(i) = 0. ;
  v(i) = 0. ;
end
ii= round(nx/2) ;
u(ii-4)=0.4 ;
u(ii-3)=0.8 ;
u(ii-2)=1. ;
u(ii-1)=1. ;
u(ii)=1. ;
u(ii+1)=1. ;
u(ii+2)=1. ;
u(ii+3)=0.8 ;
u(ii+4)=0.4 ;

plot(x,u) ;
axis([0 lg -1.1 1.1]) ;
title('Ondes, 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) ;
vp = wshift2('1D',v,+1) ;
vm = wshift2('1D',v,-1) ;
vs = (vp+vm)/2 - (0.5*dt/dx)*(up-um) ;
us = (up+um)/2 - (0.5*dt/dx)*(vp-vm) ;
v = vs ;
u = us ;
if mod(n,5) == 0
plot(x,u) ;
axis([0 lg -1.1 1.1]) ;
title('Ondes, Lax-Friedrichs, cfl=0.9')
pause(0) ;
end
%
end

pause ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Cauchy Riemann
% u,t + v,x = 0
% v,t - u,x = 0
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%parametres
%
lg = 1. ;
nx = 100 ;
dx = lg/nx ;
cfl = dx ;
dt = 0.5*cfl ;
nt = 200 ;
%
% initialisation par un etat constant legerement perturbe
%
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
us=zeros(1,nx) ;
vs=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
vp=zeros(1,nx) ;
vm=zeros(1,nx) ;
rr=rand(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  u(i) = 1. + 0.001*rr(i) ;
%  u(i) = 1. + 0.01*sin(x(i)*8*pi) ;
  v(i) = 0. ;
end

%M = moviein(1);

plot(x,u) ;
axis([0 lg -4.1 4.1]) ;
title('Cauchy-Riemann, Lax-Friedrichs, cfl=0.5')

%M(:,1)=getframe;
%mm=2;

pause ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Friedrichs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
up = wshift2('1D',u,+1) ;
um = wshift2('1D',u,-1) ;
vp = wshift2('1D',v,+1) ;
vm = wshift2('1D',v,-1) ;
vs = (vp+vm)/2 + (0.5*dt/dx)*(up-um) ;
us = (up+um)/2 - (0.5*dt/dx)*(vp-vm) ;
v = vs ;
u = us ;
if mod(n,5) == 0
plot(x,u) ;
axis([0 lg -4.1 4.1]) ;
title('Cauchy-Riemann, Lax-Friedrichs, cfl=0.5')
%M(mm)=getframe;
%mm=mm+1;
pause(0) ;
end
%
end

pause ;
%newplot;
%movie(M,3);
%pause;
% movie2avi(M,'animhyperb2') ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% initialisation
%
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
us=zeros(1,nx) ;
vs=zeros(1,nx) ;
up=zeros(1,nx) ;
um=zeros(1,nx) ;
vp=zeros(1,nx) ;
vm=zeros(1,nx) ;
for i=1:nx
  x(i) = (i-1)*dx  ;
  u(i) = 0. ;
  v(i) = 0. ;
end
ii= round(nx/2) ;
u(ii-4)=0.4 ;
u(ii-3)=0.8 ;
u(ii-2)=1. ;
u(ii-1)=1. ;
u(ii)=1. ;
u(ii+1)=1. ;
u(ii+2)=1. ;
u(ii+3)=0.8 ;
u(ii+4)=0.4 ;

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

pause(1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% marche en temps : Lax-Friedrichs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
for n=1:nt
%
up = wshift2('1D',u,+1) ;
um = wshift2('1D',u,-1) ;
vp = wshift2('1D',v,+1) ;
vm = wshift2('1D',v,-1) ;
vs = (vp+vm)/2 + (0.5*dt/dx)*(up-um) ;
us = (up+um)/2 - (0.5*dt/dx)*(vp-vm) ;
v = vs ;
u = us ;
if mod(n,5) == 0
plot(x,u) ;
axis([0 lg -4.1 4.1]) ;
title('Cauchy-Riemann, Lax-Friedrichs, cfl=0.5')
pause(0) ;
pause(0) ;
end
%
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

