// TP EF
//H. Haddar
//
// Resolution numerique de l'equation de la chaleur //
// u,t - u,xx = 0 
// schéma explicite avec condensation de masse
// schéma explicite sans condensation                              
//////////////////////////////////////////////////////

// Paramètres du problème
lg = 10. ; // intervalle en x = [-lg lg]
dx = 0.1 ; // dx = pas d'espace
cfl = 1/6.   ;     // cfl = dt/dx^2
dt = dx*dx*cfl ; // dt = pas de temps.
Tfinal=4;      // Temps final souhaité.

////////////////////////////////////////////////////////////////
// initialisation
////////////////////////////////////////////////////////////////

nx = 2*lg/dx+1 ; // nx = nombre de noeuds
x=linspace(-lg,lg,nx) ; // x(i) = -lg + i*dx

u0=max(0.,(1.-x.^2)).^2; // donnée initiale

// On trace la donnée intiale

close(); xset("window",0); // créé une nouvelle figure 0
plot2d(x,u0); // trace u0 en fonction de x
xclick() ; // on attend un click sur la fenêtre graphique pour continuer

////////////////////////////////////////////////////////////////
//  Schémas numériques
////////////////////////////////////////////////////////////////
//
u=u0 ;
us=u0;
plot2d(x,[u' u0'],leg="schema expicite@donnee initiale")


// Construction de la matrice de masse 
Mjj  = 2./3.*ones(1,nx);// la digonale de Mh
Mjjp = 1./6.*ones(1,nx-1);// les sur et sous diagonales de Mh
Mh = diag(sparse(Mjj)) + diag(sparse(Mjjp),1) + diag(sparse(Mjjp),-1); 

Mhcho = chfact(Mh) ; // factorisation de Cholesky de Mh

nt = int(Tfinal/dt);  // nt = nombre de pas de temps effectues.
Tfinal=nt*dt;    // On corrige le temps final
up= zeros(1,nx);
um= zeros(1,nx);
for n=1:nt

    // Schéma avec condensation de masse
    up(1:nx-1) = u(2:nx) ; up(nx) = 0; //u(1); // up(i) = u(i+1);
    um(2:nx) = u(1:nx-1) ; um(1) = 0; //u(nx); // um(i) = u(i-1);
    
    u=u + cfl*(up-2*u+um) ; // schéma explicite

    // Schéma sans condensation de masse
    up(1:nx-1) = us(2:nx) ; up(nx) = 0; // up(i) = u(i+1);
    um(2:nx) = us(1:nx-1) ; um(1) = 0; // um(i) = u(i-1);
    us=us + cfl*chsolve(Mhcho,(up-2*us+um));; 
    //us(nx) = 0;
    //us(1)=0;
    
   // on affiche la solution tous les 5 pas de temps
   if modulo(n,50) == 0
       clf();
        subplot(211)
        plot2d(x,[u' u0'])
        xtitle("schema explicite avec condensation")
        subplot(212)
        plot2d(x,[us' u0'])
        xtitle("schema explicite sans condensation")
   end
   xpause(100000); // on attends 100000 microsecondes 
//
end

////////////////////////////////////////////////////////////////
// comparaison solution exacte solution numérique au temps final
////////////////////////////////////////////////////////////////
uexacte = zeros(1,nx) ;
deff('[y]=noyauc(x,t)','y=exp(-x**2/(4*t))/sqrt(4*%pi*t)');
// on calcule la solution exacte en utilisant la formule des rectangles
for i=1:nx
    for j=1:nx-1
        uexacte(i) = uexacte(i) + u0(j)*dx*noyauc((i-j)*dx,Tfinal) ;
    end
end

clf()
//plotframe(rect,tics,titre);
        subplot(211)
        plot2d(x,[u' uexacte'], leg="schema explicite avec condensation@solution exacte")
        xtitle("schema explicite avec condensation")
        subplot(212)
        plot2d(x,[us' uexacte'], leg="schema explicite sans condensation@solution exacte")
        xtitle("schema explicite sans condensation")
//plot2d(x,[u' us' uexacte'],leg="schema explicite avec condensation@schema explicite sans condensation@solution exacte")
