////////////////////////////////////////////////////////////////
// Copyright G. Allaire, Octobre 2017
//
// Schéma centré instable
// pour l'equation de la chaleur:
//
////////////////////////////////////////////////////////////////
//
// equation de la chaleur
// u,t - u,xx = 0
// schema centre instable
//
////////////////////////////////////////////////////////////////
//
lg = 10. ; // lg = demie longueur du domaine
dx = 0.2 ; // dx = pas d'espace
nx = lg/dx ; // nx = demi nombre de mailles
cfl = 0.1 ; // cfl
dt = dx*dx*cfl ; // dt = pas de temps
nt = 20 ; // 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 ;
u1= u0 ; // solution au temps n
u2 = u0 ; // solution au temps n-1

id=color(255,0,0) ;

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,-1;lg,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 centré"],[[id;1],[-5;1]],opt="ur",font_size=6) ;
title ('EQUATION DE LA CHALEUR, SCHEMA CENTRE, CFL=0.1' ,'fontsize',6) ;


halt() ;


////////////////////////////////////////////////////////////////
// schema centre instable
////////////////////////////////////////////////////////////////
//
// premier pas de temps
//
up = shift('+1',u0) ;
um = shift('-1',u0) ;

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

//
// pas de temps suivants
//
for n=2:nt
//
up = shift('+1',u1) ;
um = shift('-1',u1) ;

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

if modulo(n,1) == 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,-1;lg,2];
        a.font_size=6;
        a.thickness=2;
        plot(x,u,'b-d',x,uexacte,'r-')
	legends(["solution exacte";"schéma centré"],[[id;1],[-5;1]],opt="ur",font_size=6) ;
          title ('EQUATION DE LA CHALEUR, SCHEMA CENTRE, CFL=0.1' ,'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.5;lg,2];
        a.font_size=6;
        a.thickness=2;
        plot(x,u,'b-d',x,uexacte,'r-')
	      legends(["solution exacte";"schéma centré"],[[id;1],[-5;1]],opt="ur",font_size=6) ;
        title ('EQUATION DE LA CHALEUR, SCHEMA CENTRE, CFL=0.1' ,'fontsize',6) ;