////////////////////////////////////////////////////////////////
// Copyright G. Allaire, Octobre 2017
//
// Schéma explicite avec CFL=0.51
// pour l'equation de la chaleur:
//
////////////////////////////////////////////////////////////////
//
// equation de la chaleur
// u,t - u,xx = 0
// schema explicite: cfl=0.51
//
////////////////////////////////////////////////////////////////
//
lg = 10. ; // lg = demie longueur du domaine
dx = 0.2 ; // dx = pas d'espace
nx = lg/dx ; // nx = demi nombre de mailles
cfl = 0.51 ; // cfl
dt = dx*dx*cfl ; // dt = pas de temps
nt = 100 ; // nt = nombre de pas de temps effectues
//
// initialisation
//
x=zeros(1,2*nx+1) ;
u0=zeros(1,2*nx+1) ;
for i=1:2*nx+1
  x(i) = (i-nx-1)*dx  ;
  u0(i) = max(0.,1.-x(i)**2) ;
end
u=u0 ;
up=u0 ;
um=u0 ;
uexacte=u0 ;
clf()

a=get("current_axes"); //get the handle of the newly created axes
a.axes_visible="on"; // makes the axes visible
a.font_size=6; //set the tics label font size
a.x_location="bottom"; //set the x axis position
a.data_bounds=[-lg,-0.4;lg,1.2]; //set the boundary values for the x, y and z coordinates.
a.sub_tics=[5,0];
a.thickness=2;
//a.labels_font_color=5;

plot(x,u,'b-d')
legends(["solution exacte";"schéma explicite"],[[id;1],[-5;1]],opt="ur",font_size=6) ;
title ('EQUATION DE LA CHALEUR, SCHEMA EXPLICITE, CFL=0.51' ,'fontsize',6) ;

halt() ;

////////////////////////////////////////////////////////////////
// schema explicite: cfl=0.51
////////////////////////////////////////////////////////////////
//
for n=1:nt
//
up = shift('+1',u) ;
um = shift('-1',u) ;

u=u + dt/(dx*dx)*(up+um-2*u) ;

if modulo(n,5) == 0
// comparaison avec la solution exacte
        uexacte = zeros(1,2*nx+1) ;
        noyau = zeros(1,2*nx+1) ;
        for i=1:2*nx+1
            noyau(i) = exp(-((i-1)*dx)**2/(4*n*dt)) ;
            if noyau(i) < 1.e-14 then
            kmax = i ;
            break
            end
        end
        for i=1:2*nx+1
            jmin = max(1, i-kmax) ;
            jmax = min(2*nx+1,i+kmax) ;
            for j=jmin:jmax
                 k =1 + int(abs(i-j)) ;
                 uexacte(i) = uexacte(i) + u0(j)*dx*noyau(k) ;
            end
            uexacte(i) = uexacte(i)/sqrt(4*%pi*n*dt) ;
        end
        clf()
        a=get("current_axes");
        a.data_bounds=[-lg,-0.4;lg,1.2];
        a.font_size=6;
        a.thickness=2;
        plot(x,u,'b-d',x,uexacte,'r-')
	      legends(["solution exacte";"schéma explicite"],[[id;1],[-5;1]],opt="ur",font_size=6) ;
        title ('EQUATION DE LA CHALEUR, SCHEMA EXPLICITE, CFL=0.51' ,'fontsize',6) ;
        sleep(300) ;
end
//
end

// comparaison au temps final
        uexacte = zeros(1,2*nx+1) ;
        noyau = zeros(1,2*nx+1) ;
        for i=1:2*nx+1
            noyau(i) = exp(-((i-1)*dx)**2/(4*nt*dt)) ;
            if noyau(i) < 1.e-14 then
            kmax = i ;
            break
            end
        end
        for i=1:2*nx+1
            jmin = max(1, i-kmax) ;
            jmax = min(2*nx+1,i+kmax) ;
            for j=jmin:jmax
                 k =1 + int(abs(i-j)) ;
                 uexacte(i) = uexacte(i) + u0(j)*dx*noyau(k) ;
            end
            uexacte(i) = uexacte(i)/sqrt(4*%pi*nt*dt) ;
        end
        clf()
        a=get("current_axes");
        a.data_bounds=[-lg,-0.4;lg,1.2];
        a.font_size=6;
        a.thickness=2;
        plot(x,u,'b-d',x,uexacte,'r-')
        legends(["solution exacte";"schéma explicite"],[[id;1],[-5;1]],opt="ur",font_size=6) ;
        title ('EQUATION DE LA CHALEUR, SCHEMA EXPLICITE, CFL=0.51' ,'fontsize',6) ;