////////////////////////////////////////////////////////////////
// Copyright G. Allaire, Octobre 2017
//
// Schéma implicite
// pour l'equation de la chaleur:
//
////////////////////////////////////////////////////////////////
//
// equation de la chaleur
// u,t - u,xx = 0
// schema implicite: cfl=2.
//
////////////////////////////////////////////////////////////////
//
lg = 10. ; // lg = demie longueur du domaine
dx = 0.2 ; // dx = pas d'espace
nx = lg/dx ; // nx = demi nombre de mailles
cfl = 2.0 ; // cfl
dt = dx*dx*cfl ; // dt = pas de temps
nt = 30 ; // 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 ;
mat=zeros(2*nx+1,2*nx+1) ;

for i=2:2*nx
  mat(i,i) = 1. + 2*dt/(dx*dx) ;
  mat(i,i+1) =  -dt/(dx*dx) ;
  mat(i,i-1) = -dt/(dx*dx) ;
end
mat(1,1) = 1. + 2*dt/(dx*dx) ;
mat(1,2) =  -dt/(dx*dx) ;
mat(2*nx+1,2*nx) = -dt/(dx*dx);
mat(2*nx+1,2*nx+1) = 1. + 2*dt/(dx*dx) ;
smat = sparse(mat) ;
[R,P] = spchol(smat) ; // factorisation de Cholesky


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,-0.4;lg,1]; //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 implicite"],[[id;1],[-5;1]],opt="ur",font_size=6) ;
title ('EQUATION DE LA CHALEUR, SCHEMA IMPLICITE, CFL=2.0' ,'fontsize',6) ;

halt() ;

////////////////////////////////////////////////////////////////
// schema implicite: cfl=2.0
////////////////////////////////////////////////////////////////
//
for n=1:nt
//
//u = chsolve(spcho,u) ; // resolution du systeme lineaire
b=u';
b = P*(R'\(R\(P'*b)));
u=b';

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