%  test5.m    demo1 pour l X    demo1_2 pour octave 6.4  
%  advection avec conditions limites periodiques
%  schema instable ;  profil regulier,
%  fd,   12 novembre, 14 decembre  1999, 4 janvier 2K, 16 et ?? fevrier 2k, 
%  15 janvier 2002, 17 janvier 2005, 05 fevrier 2006, 24, 25 fevrier 2025
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cfl=0.3;  tmax = 0.6 ;  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 ;

plot(x,u0(:,1),"ob;condition initiale;" )
pause (1.5)

%% gset yrange[-0.6:1.6] ; gset xrange[0.: 1.0] ; 
%%%   axis ([0 1 -.6 1.6]) 
%%% title ('Schema instable, conditions limites periodiques') 
%%% xlabel (titx ) ; plot(x,u0,'*;condition initiale discretisee;') ; pause 
%%% 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); zzz=V(n:n,1)-V(1:1,1) ; 
U=V;W=V;  for i=1:m ; 
W(2:n-1,1) = U(2:n-1,1) -(0.5*cfl)* (U(3:n,1)-U(1:n-2,1)) ;
W(n:n,1)   = U(n:n,1) -(0.5*cfl)* (U(2:2,1)-U(n-1:n-1,1)) ;
W(1:1,1)   = W(n:n,1) ;   U=W ;
% 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,ue(:,1),"ob;;"  ,   x,U(:,1),'*;;' ) 
end ;

plot(x,ue(:,1),"b;solution exacte;", x,U(:,1),'*;solution approchee;',x,U(:,1),"k;;") ; 
title('schema  explicite centre instable')
axis ([0 1 -.4  2.2])
print -djpg figure-1.jpg
print -dpdf figure-1.pdf

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


