////////////////////////////////////////////////////////////////
// Copyright G. Allaire, Janvier 2001
//
// Sur la necessite de l'hyperbolicite pour les systemes 
// de lois de conservation:
// calcul numerique sur des systemes lineaires hyperboliques 
// ou non, comparaison equation des ondes et de Cauchy-Riemann 
// (condition aux limites periodique)
//
////////////////////////////////////////////////////////////////
//
// exemple de la sinusoide (corde vibrante)
//
////////////////////////////////////////////////////////////////
//
// equation des ondes
// u,t + v,x = 0
// v,t + 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 ;
cfl = dx ;
// dt = pas de temps
dt = 0.9*cfl ;
// nt = nombre de pas de temps effectues
nt = 500 ;
//
// initialisation
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
us=zeros(1,nx) ;
vs=zeros(1,nx) ;
ue=zeros(1,nx) ;
ve=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

xbasc()
tics=[4,10,4,10];
plotframe([0,-1.1,lg,1.1],tics);
plot2d(x,u,1,"000")
xtitle ('ondes, Lax-Friedrichs' ,' ',' ') ;

halt() ;
////////////////////////////////////////////////////////////////
// marche en temps : Lax-Friedrichs
//
// solution exacte:
// u(t,x) = ( u0(x-t) + u0(x+t) )/2
// v(t,x) = ( u0(x-t) - u0(x+t) )/2
////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
up = shiftp('+1',u) ;
um = shiftp('-1',u) ;
vp = shiftp('+1',v) ;
vm = shiftp('-1',v) ;

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 modulo(n,5) == 0

for i=1:nx
    xi = (i-1)*dx  ;
    ue(i) = ( sin((xi-n*dt)*2*pi) + sin((xi+n*dt)*2*pi) )/2 ;
end
xbasc() 
plotframe([0,-1.1,lg,1.1],tics);
xset("thickness",2) ;
plot2d(x,u,[1,1],"100","Lax-Friedrichs")
xset("thickness",1) ;
plot2d(x,ue,[2,4],"100","exact")
xtitle('ondes, Lax-Friedrichs, cfl=0.9',' ',' ');
xpause(100000) ;

end
//
end


halt() ;
////////////////////////////////////////////////////////////////
//
// initialisation pour Cauchy Riemann
//
// dt = pas de temps
dt = 0.5*cfl ;
// nt = nombre de pas de temps effectues
nt = 200 ;
//
x=zeros(1,nx) ;
u=zeros(1,nx) ;
v=zeros(1,nx) ;
us=zeros(1,nx) ;
vs=zeros(1,nx) ;
ue=zeros(1,nx) ;
ve=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) ;
  v(i) = 0. ;
end

xbasc()
tics=[4,10,4,10];
plotframe([0,-1.1,lg,3.1],tics);
plot2d(x,u,1,"000")
xtitle ('Cauchy-Riemann, Lax-Friedrichs' ,' ',' ') ;

halt() ;
////////////////////////////////////////////////////////////////
// marche en temps : Lax-Friedrichs
//
////////////////////////////////////////////////////////////////
// 
for n=1:nt
//
up = shiftp('+1',u) ;
um = shiftp('-1',u) ;
vp = shiftp('+1',v) ;
vm = shiftp('-1',v) ;

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 modulo(n,5) == 0

xbasc() 
plotframe([0,-1.1,lg,3.1],tics);
xset("thickness",2) ;
plot2d(x,u,[1,1],"100","Lax-Friedrichs")
xtitle('Cauchy-Riemann, Lax-Friedrichs, cfl=0.9',' ',' ');
xpause(100000) ;

end
//
end