// On cherche à résoudre u_t + c u_x =0 pour c>0
stacksize(100000000);


c = 1; // valeur de la vitesse de transport
L = 10; // domaine [-L,L]
M = 100; // nombre de points de discrétisation en espace
dx = 2*L/M; // pas d'esapce
x = (-L:dx:L)'; 
Tf = 3; // temps final
dt = 0.01; // pas de temps
a = c*dt/dx;


u0 = exp(-x.^2); // condition initiale régulière
//u0 = (1-abs(1-x)).*(x>=0).*(x<=2); // condition initiale pas régulière
//u0 = (x>=-1).*(x<=1); // condition initiale discontinue

// On utilise un schéma décentré
band1 = sparse(diag( ones(M+1,1)); // construit une matrice diagonale
band2 = sparse(diag( ones(M,1),-1) ); // construit une matrice sous-diagonale
mat = (1-a)*band1 + a*band2;
mat(1,M+1) = a; // condition de bords périodique
t = 0;
while t<Tf
    t = t+dt;
    u0 = mat*u0; // solution approchée
    u1 = exp(-(x-c*t).^2); // sol exacte avec CI régulière
    //u1 = (1-abs(1-(x-c*t))).*(x-c*t>=0).*(x-c*t<=2); // sol exacte avec CI pas régulière
    //u1 = (x-c*t>=-1).*(x-c*t<=1); // sol exacte avec CI discontinue
    drawlater();
    clf();
    plot(x,u0,'d',x,u1,2)
    drawnow();
end


//// On utilise un schéma dit de Lax-Friedrichs
//band2 = sparse(diag( ones(M,1),-1)); // construit une matrice sous-diagonale
//band3 = sparse(diag( ones(M,1),1) ); // construit une matrice sur-diagonale
//mat = ((1+a)/2)*band2+((1-a)/2)*band3;
//mat(1,M+1) = (1+a)/2; // condition de bords périodique
//mat(M+1,1) = (1-a)/2; // condition de bords périodique
//t = 0;
//while t<Tf
//    t = t+dt;
//    u0 = mat*u0; // solution approchée
//    //u1 = exp(-(x-c*t).^2); // sol exacte avec CI régulière
//    //u1 = (1-abs(1-(x-c*t))).*(x-c*t>=0).*(x-c*t<=2); // sol exacte avec CI pas régulière
//    u1 = (x-c*t>=-1).*(x-c*t<=1); // sol exacte avec CI discontinue
//    drawlater()
//    clf;
//    plot2d(x,u0)
//    plot2d(x,u1,2)
//    drawnow()
//end

//// On utilise un schéma dit de Lax-Wendroff
//band1 = sparse(diag( ones(M+1,1))); // construit une matrice diagonale
//band2 = sparse(diag( ones(M,1),-1)); // construit une matrice sous-diagonale
//band3 = sparse(diag( ones(M,1),1) ); // construit une matrice sur-diagonale
//mat = (1-a^2)*band1+a*((1+a)/2)*band2-a*((1-a)/2)*band3;
//mat(1,M+1) = a*(1+a)/2; // condition de bords périodique
//mat(M+1,1) = -a*(1-a)/2; // condition de bords périodique
//t = 0;
//while t<Tf
//    t = t+dt;
//    u0 = mat*u0; // solution approchée
//    //u1 = exp(-(x-c*t).^2); // sol exacte avec CI régulière
//    u1 = (1-abs(1-(x-c*t))).*(x-c*t>=0).*(x-c*t<=2); // sol exacte avec CI pas régulière
//    //u1 = (x-c*t>=-1).*(x-c*t<=1); // sol exacte avec CI discontinue 
//    drawlater()
//    clf;
//    plot2d(x,u0)
//    plot2d(x,u1,2)
//    drawnow()
//end
