%  test11.m           demo2 pour l X       demo2_2 pour octave 6.4 
%  advection avec conditions limites periodiques
%  schema implicite centre d ordre un en temps  ;  profil regulier  , 
%  fd, 18 fevrier 2k, 15 janvier 2002, 17 janvier 2005
%      05 fevrier 2006, 24, 25 fevrier 2025.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cfl=0.8 ;  tmax = 1.0 ;  n1=101;   dx=1/(n1-1);
% cfl=1.2 ;  tmax = 1.0 ;  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 ; ; 
%  u0(j:j,1) = sin(2*xx*pi) ; 
end  ;
% label (chcfl) at 0.1,-0.1
% gset yrange[-0.2:1.2] ; gset xrange[-0: 1.0] ; 
%%% axis ([0 1 -.2 1.2]) 
%%% title (' Advection  profil regulier     Cns lims periodiques    Schema impli centre d ordre un en temps  ') ; 
%%% xlabel (titx ) ;

plot (x,u0,'*;condition initiale;') 
pause (1.5)

%%% xlabel (titx1 ) ;
ue = u0;    dt=cfl*dx;    t=0:dt:tmax;   m=length(t)-1 ; 
A=zeros(n-1,n-1) ;
for i=1:n-1 ; A(i:i,i:i)=1 ; end ; 
for i=1:n-2 ; A(i:i,i+1:i+1)=0.5*cfl ; A(i+1:i+1,i:i)=-0.5*cfl ; end ; 
A(1:1,n-1:n-1)=-0.5*cfl ; A(n-1:n-1,1:1)=0.5*cfl ; A ;
V=zeros(n-1,1);  V(1:n-1,1)=u0(1:n-1) ; 
U=zeros(n,1) ; U(1:n-1,1)=  V(1:n-1,1) ; U(n:n)=U(1:1) ; W=U ; 
%  debut de la boucle en temps
% m = 3
for i=1:m ; 
W(2:n-1,1) = U(2:n-1,1)    ;
W(n:n,1)   = U(n:n,1)   ;
W(1:1,1)   = W(n:n,1) ; V(1:n-1,1)=W(1:n-1,1) ;
%   fin du predicteur
VVV=A\V ; ee=V-A*VVV ;  
for bi = 1:10 ; 
VV=A\ee ; VV1 = VVV + VV ; ee=V-A*VVV ;
 VVV = VV1 ;  end ;
U(1:n-1,1)=VVV(1:n-1,1) ;  U(n:n,1)=U(1:1,1) ;
%   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  implicite centre')
axis ([0 1 -.01  1.01])

print -djpg figure-2.jpg
print -dpdf figure-2.pdf





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