% #!/usr/bin/env octave
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Modélisation numérique d'une équation différentielle non-linéaire d'ordre 1
% Ce code a pour but de mettre en oeuvre un schéma numérique simple d'ordre 1 pour réroudre numériquement une équation différentielle non-linéaire d'ordre 1.
% Programmée en Matlab v11+ par Lilian Besson (ENS Cachan), pour une collaboration avec Nicolas Vernier (Université Paris-Sud), à Mahindra École Centrale, (C) 2014.
% Distribué sous les termes de la licence GPLv3.

% @date Dimanche 21 Septembre 2014
% @author Lilian BESSON
% @see <http://perso.crans.org/besson/publis/PDE_Nicolas/> Les sources en lignes sur ma page.
% @see <http://math.stackexchange.com/q/897299> La question que j'avais posée sur Math StackExchange.
clear all; clc;

% 1 si on veut afficher symétriquement de -1 à 1, 0 si de 0 à 1.
printTout = 1
% 1 si on veut afficher la moyenne de Fréchet, 0 sinon.
printFenetre = 1

%
% Définition des constantes
%
k = 1.
x0 = 0.;
x1 = 1.;
lx = x1 - x0;

n = 10000
h = lx / n;

dt = 0.00005
c = h * dt * sqrt(k)

jmax = 10000
jskip = 10

% Applatir le profil initial
beta = 0.05
beta = 1.0

%
% Initialisation des tableaux
%


if (1 == printTout)
    % Afficher de -1 à 1 (double le temps de calcul)
    x = -x1 : h : x1;
else
    % Afficher seulement de 0 à 1 (plus rapide, moins visuel)
    x = x0 : h : x1;
    xmoins = -x1 : h : x0;
end;
printf('\nTaille de x :\n');
length(x)


% Profil initial affine : translation (comme le montre la méthode des caractéristiques)
y0 = beta .* (1. - x);
% Profil initial parabolique : évolution d'un pic central
y0 = beta .* (1. - (x .^ 2));

N = length(y0);
imoinsun = @(j) max(1,j-1);
iplusun = @(j) min(N,j+1);


% Largeur de la fenêtre utilisée pour réaliser la moyenne de Fréchet
fenetre = 50
hFenetre = h*fenetre

% Taille du vecteur moyenné
nFenetre = floor((N-1) / fenetre)
% Initialisation de la matrice de moyenne
iFenetre = 1:(N-1);
length(iFenetre)
jFenetre = ceil(iFenetre ./ fenetre);
length(jFenetre)
% Grâce à une matrice creuse (gain de 500/n % !)
matFenetre = sparse(jFenetre, iFenetre, 1/fenetre);
length(matFenetre)

if (1 == printTout)
    xFenetre = -x1 : hFenetre : (x1 - hFenetre);
else
    xFenetre = x0 : hFenetre : (x1 - hFenetre);
end;
printf('\nTaille de xFenetre :\n');
length(xFenetre)

%
% Lancement de la simulation
%

hold off;
% hold on;

% Changer ici le type de schéma (1, 2, ou 3, ou 1:3)
for schemaN = 3
    % Initialisation de la simulation pour ce schéma
    switch schemaN
        case 1
            schema = 'droit'
        case 2
            schema = 'centre'
        otherwise
            schema = 'gauche'
    end
    clf;
    % Init y
    y = y0;
    printf('\nTaille de y :\n');
    length(y)
    % Init yFenetre
    yFenetre = matFenetre * (y0(1:(N-1))');
    printf('\nTaille de yFenetre :\n');
    length(yFenetre)
    % Option des graphiques
    if (1 == printFenetre)
        % On affiche les moyennes locales
        HandlePlot = plot(xFenetre, yFenetre, '-r', 'YDataSource', 'yFenetre');
        ylabel(sprintf('Position y moyennée du front de matière, yFenetre > 0 (adimensionée), avec %i valeurs.', nFenetre));
        if (1 == printTout)
            xlabel(sprintf('Position x, -1 < xFenetre < 1 (adimensionée), avec %i points.', nFenetre));
        else
            xlabel(sprintf('Position x, 0 < xFenetre < 1 (adimensionée), avec %i points.', nFenetre));
        end;
    else
        % On affiche le vecteur entier
        HandlePlot = plot(x, y, '-r', 'YDataSource', 'y');
        ylabel(sprintf('Position y du front de matière, y > 0 (adimensionée), avec %i valeurs.', N));
        if (1 == printTout)
            xlabel(sprintf('Position x, -1 < x < 1 (adimensionée), avec %i points.', n));
        else
            xlabel(sprintf('Position x, 0 < x < 1 (adimensionée), avec %i points.', n));
        end;
    end;
    grid on;
    % Début de la simulation pour ce schéma
    for j = 0:jmax
        % Actualisation du vecteur y (en une seule étape vectorielle)
        switch schemaN
            case 1
                % 1) Schéma gauche
                y = y + c ./ (sqrt(h.^2 + ( y - y(imoinsun(1:N)) ).^2));
            case 2
                % 2) Schéma droit
                y = y + c ./ (sqrt(h.^2 + ( y(iplusun(1:N)) - y ).^2));
            case 3
                % 3) Schéma centré
                y = y + c ./ (sqrt(h.^2 + 0.25 .* ( y(iplusun(1:N)) - y(imoinsun(1:N)) ).^2));
            otherwise
                % 4) Ne fait rien.
        end;
        %
        % Actualisation du vecteur y (en une boucle lente et lourdre)
        % for i = 1:N
            % 1) Schéma gauche
            % y(i) = y(i) + c ./ (sqrt(h.^2 + ( y(i) - y(imoinsun(i)) ).^2));
            % 2) Schéma droit
            % y(i) = y(i) + c ./ (sqrt(h.^2 + ( y(iplusun(i)) - y(i) ).^2));
            % 3) Schéma centré
            % y(i) = y(i) + c ./ (sqrt(h.^2 + 0.25*( y(iplusun(i)) - y(imoinsun(i)) ).^2));
        % end;
        % Affichage du y courant (de temps en temps)
        if (mod(j,jskip)==0)
            title(sprintf('Simulation du déplacement du front de matière au cours du temps (étape j=%i), via un schéma %s \n(C) Lilian Besson et Nicolas Vernier, 2014.\n', j, schema));
            if (1 == printFenetre)
                % Régularisation de la solution (fondée sur une moyenne locale de Fréchet)
                yFenetre = matFenetre * (y(1:(N-1))');
                set(HandlePlot, 'YData', yFenetre);
            else
                set(HandlePlot, 'YData', y);
            end;
            drawnow;
            if (mod(j,100*jskip)==0)
                saveas(HandlePlot, sprintf('plot_%s_%i.png', schema, j), 'png' );
            end;
        end;
    end;
    printf('\n\nFin de la simulation pour le schéma %s.\n', schema);
    % http://www.mathworks.in/help/matlab/ref/saveas.html
    saveas(HandlePlot, sprintf('plot_%s_last.png', schema) );
    pause();
end;

fprintf('\n\nFin de la simulation numérique.\nProgramme Matlab v2011+ par Lilian Besson (ENS Cachan), pour une collaboration avec Nicolas Vernier (Université Paris-Sud), à Mahindra École Centrale, (C) 2014.\n');
% pause();
