%  test4.m        demo 4 pour l X         demo4_2 pour octave 6.4
%  advection avec conditions limites periodiques
%  schema de Lax-Wendroff ;  profil regulier  ,  
%  fd, 12 nov  99, 4 janvier 2K, 18 fevrier 2K, 15 janvier 2002, 
%      17 janvier 2005, 05 fevrier 2006, fevrier 2025
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  cfl=1.1 ;  tmax = .4 ;  n1=101;   dx=1/(n1-1);
cfl=.8 ;  tmax = 1. ;  n1=101;   dx=1/(n1-1);
%%%%  [chcfl,msg] = sprintf (' nombre de Courant : %5.2f',cfl);
%%%%  [chn1,msg] = sprintf ('        nombre de mailles : %3i',n1-1);
%%%%  [chtps,msg] = sprintf ('        temps final  : %5.2f',tmax);chtps;
%%%%  titx = [chcfl , chn1] ; titx;
%%%%  titx1 = [chcfl , chn1 , chtps] ; titx1;
x=0:dx:1 ;   n=length(x) ; 
%           profil regulier de classe C2 comme condition initiale
pdebut = [16*314  -4*794  539  0 0 0 ] ; 
pfin   = [16*350   -4*878 589  0 0 0 ] ; 
p      = [1 -1 0 1 ] ;
u0=zeros(n,1);      u0(2:n1/2-1,1)=1;  
for j=1:n ; xx = (j-1)*dx  ;
if xx<0.25     ; u0(j:j,1)= polyval(pdebut,xx);
elseif xx>0.75 ; u0(j:j,1)= polyval(pfin,1-xx); 
else  u0(j:j,1)= polyval(p,xx-0.5) ;  end ; end ;  
%%% gset yrange[-0.3:1.3] ; gset xrange[0.: 1.0] ; 
%%% axis ([0 1 -.3 1.3]) 
%%%%  title ('Condition initiale discretisee ') ;  xlabel (titx ) ;

plot (x,u0,'*;condition initiale;') 
axis ([0 1 -.01  1.01])  
pause (1.5)


%%% xlabel (titx1 ) ;
ue = u0;    dt=cfl*dx;    t=0:dt:tmax;   m=length(t)-1 ; 
V=zeros(n,1);  V(1:n,1)=u0(1:n);
%%%%  title ('Advection d un profil regulier     Conditions limites  periodiques   Schema de Lax Wendroff') 
%%%%% gset yrange[-0.3:1.3]
U=V;  for i=1:m ; 
W(1:n-1,1) = 0.5*(U(1:n-1,1)+U(2:n,1))-(0.5*dt/dx)*   (U(2:n,1)-U(1:n-1,1));
% W(n:n,1)   = 0.5*(U(n:n,1)+U(1:1,1))  -(0.5*dt/dx)*   (U(1:1,1)-U(n:n,1));
W(n:n,1)   = W(1:1,1) ; 
V(2:n,1) = U(2:n,1) -dt/dx*(W(2:n,1)-W(1:n-1,1));
V(1:1,1) = V(n:n,1) ;     U=V;
%   solution exacte 
    for j=1:n ; zz =  (j*dx- i*dt) / dx ; mm = round (zz) ;
    if mm<=0 ;    p=mm+n-1; mu=mm ; 
    while p<=0 ; pp=mm ;  p=pp+n-1 ; mm=p ; end ; mm=p; 
    zz = zz + mm-mu ; end ; 
    if mm==1 ;  mmm1=n-1; mmp1 = 2 ; else ; if mm==n; mmm1=n-1; mmp1=2 ;
    else mmm1=mm-1 ; mmp1 = mm+1 ; end ; end ; 
    if mm <= zz ; ue(j:j,1) = (1-zz+mm)*u0(mm)+(zz-mm)*u0(mmp1);
    else  ue(j:j,1) = (mm-zz)*u0(mmm1)+ (1+zz-mm)*u0(mm) ; end ;  end; 
%   fin de la solution exacte
%   plot(x,V(:,1),'*',x,ue(:,1),"b") ,  end ;
%%%%        plot(x,ue(:,1),"ob;;" , x,U(:,1),'*;;' ) 
end ;

plot(x,ue(:,1),"b;solution exacte;", x,U(:,1),"*k;solution approchee;") ; 
axis ([0 1 -.01  1.01])
title('schema de Lax Wendroff')
print -djpg figure-4.jpg
print -dpdf figure-4.pdf



% plot(x,ue(:,1),"b;solution exacte;", x,U(:,1),"*;solution approchee;",x,U(:,1),"r;;")
%
%  gset terminal postscript ;  gset output "octave/laxwendroff.ps" ;
%  plot(x,ue(:,1),"b;solution exacte;", x,U(:,1),'*;solution approchee;',x,U(:,1),"r;;")

