(**

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ésoudre numériquement une équation différentielle non-linéaire d'ordre 1.

Programmée en OCaml v4+ 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 Vendredi 12 Septembre 2014 @author Lilian BESSON @see <http://perso.crans.org/besson/publis/PDE_09_2014/> Les sources en lignes sur ma page. @see <http://math.stackexchange.com/q/897299> La question que j'avais posée sur Math StackExchange. *)



(**

Définition des constantes.

*)


(** La valeur réelle de k n'a pas d'importance ? *)

let k = 1.0;; (* FIXME *)

(** Intervalle de résolution. *)

let x0 = 0. and x1 = 1.;;
let lx = x1 -. x0;;


(** Nombre de points. *)

let n = 1500;; (* FIXME *)
let h = lx /. (float n);;

(** Calculer la valeur flottante du ième point x_i. *)

let x_i i = x0 +. ((float i) *. h);;


(** Cette valeur dt peut avoir une importance ! *)

let dt = 0.0001;; (* FIXME *)

(** Calculer la valeur flottante du jème instant t_i. *)

let t_j j = ((float j) *. dt);;

(** Nombre d'étape maximum. *)

let jmax = 3000;; (* FIXME *)


(**

Fonctions pour afficher.

*)


let print = Format.printf;;
let sprint = Format.sprintf;;
let fprint = Printf.fprintf;;

(** Extension par défaut utilisée pour produire les graphiques (soit PNG soit SVG). *)

let dftlext = "svg";;
let dftlext = "png";;

(**

Différentes fonctions d'écritures.

*)


(** Génération et écriture du fichier script Gnuplot. *)

let print_gnuplot ?(verb=false) ?(script="plotdata.gnuplot") ?(data="data.txt") ?(plot="plot") ?(ext=dftlext) ?(t=0.) ?(n=n) ?(ymin=0.) ?(ymax=0.02) () =
        let outchannel = open_out script in
        if verb then print "\n  1. Écriture du script GNUPlot v4.6+ dans le fichier '%s' (taille 1366×768, type %s, yrange=[%g:%g], fichier de données %s, image de sortie %s.%s)..." script ext ymin ymax data plot ext;
        fprint outchannel "#! /usr/bin/env gnuplot
# -*- coding:utf8; mode:gnuplot -*-
# Script GNUPlot v4.6+, (C) Lilian Besson, 2014, GPLv3 Licenced
# Modélisation numérique d'une équation différentielle non-linéaire d'ordre 1.
set terminal %s size 1366,768
set style line 1 linewidth 0.5
set output '%s.%s'
set grid
set title 'Graphe de y(., tj) (pour tj=%g, avec %i points). (C) Lilian Besson et Nicolas Vernier, 2014'
set xrange [-1:1]
set yrange [%g:%g]
set xlabel 'Position horizontale (-1 < x < 1) avec %i points.'
set ylabel 'Position verticale du front de matière magnétisée (0 < y < %g).'
set label 'Front de matière magnétisée (au temps t=%g).' at 0.4,%g
plot '%s' with lines ls 1
# Fin du script GNUPlot v4.6+.\n"
 ext plot ext t (n+1) ymin ymax n ymax t (0.9 *. ymax) data;
        flush outchannel;
        close_out outchannel;
;;

(** Génération et écriture du fichier de données brutes à partir du vecteur y. *)

let print_y ?(verb=false) ?(t=0.) ?(n=n) ?(dt=dt) ?(data="data.txt") y () =
        let outchannel = open_out data in
        if verb then print "\n  2. Écriture des données dans le fichier '%s' (n=%i, dt=%g, t=%g)..." data n dt t;
        fprint outchannel "#! /usr/bin/env cat
# -*- coding:utf8; mode:plain -*-
# Données numériques en mode texte, n=%i, dt=%g, t=%g
# (C) Lilian Besson, 2014, GPLv3 Licenced
# Modélisation numérique d'une équation différentielle non-linéaire d'ordre 1.\n"
 n t dt;
        for i = n downto 1 do
                fprint outchannel " \t%g \t%g\n" (-1. *. (x_i i)) y.(i);
        done;
        for i = 0 to n do
                fprint outchannel " \t%g \t%g\n" (x_i i) y.(i);
        done;
        fprint outchannel "# Fin du fichier de données numériques en mode texte.";
        flush outchannel;
        close_out outchannel;
;;

(** Exécuter le script Gnuplot (par un Sys.command). *)

let launch_gnuplot ?(verb=false) ?(inter=false) ?(script="plotdata.gnuplot") ?(t=0.) () =
        if verb then print "\n  3. Exécution du script '%s' par GNUPlot v4.6+..." script;
        if inter then
                Sys.command (sprint "/usr/bin/gnuplot \"%s\" && notify-send 'Succès gnuplot pour t=%g'" script t)
        else
                Sys.command (sprint "/usr/bin/gnuplot \"%s\"" script);;

(** Affichage de l'image produite, par /usr/bin/eog (par un Sys.command). *)

let view_eog ?(verb=false) ?(plot="plot") ?(ext=dftlext) () =
        if verb then print "\n  4. Affichage de l'image produit '%s.%s' par /usr/bin/eog..." plot ext;
        Sys.command (sprint "/usr/bin/eog \"%s.%s\" &" plot ext);;

(** Combine toutes les étapes précédentes. *)

let print_and_launch ?(schema="gauche") ?(verb=false) ?(inter=false) ?(ext=dftlext) ?(j=0) y () =
        print_gnuplot ~verb:false ~plot:(sprint "%s/plot_%05d" schema j) ~ext:ext ~t:(t_j j) ~ymax:(1.1 *. y.(0)) ();
        print_y ~verb:false ~t:(t_j j) y ();
        assert(0 = launch_gnuplot ());
        if inter then assert(0 = view_eog ~verb:verb ~plot:(sprint "%s/plot_%05d" schema j) ~ext:ext ());
;;


(**

Initialisation du vecteur y_i.

*)


(** Pour applatir la courbe au début. *)

let beta = 0.1
let beta = 1.0
(* FIXME *)

(** Profil initiale de la fonction y, affine, afin de vérifier que c'est le seul cas où la méthode des caractéristique permet de résoudre de façon formelle l'EDP. *)

let y_0 x = (1. -. x) *. beta ;;
(** Profil initiale de la fonction y, parabolique comme observé réellement. *)

let y_0 x = (1. -. (x ** 2.)) *. beta ;;

(** Créer le tableau initial. *)

let makeym () = Array.init (n+1) (fun i -> y_0 (x_i i));;

(** Précalcul de cette constante c = h × dt × sqrt(k). *)

let c = h *. dt *. sqrt(k);;

(** On utilise la symétrie naturelle du problème pour calculer y_{i-1} quand i=0. *)

let iMoins1 = function 0 -> 1 | i -> i-1;;

(** On utilise la symétrie naturelle du problème pour calculer y_{i+1} quand i=n. *)

let iPlus1 ?(n=n) = function i when i=n -> i | i -> i+1;;


(**

Schéma numérique.

*)


(** Crée un nouveau tableau. @param schema précise le type de schéma (décentré gauche, droit, ou centré). Apparemment, seul le gauche est physiquement valide. *)

let newnextY ?(schema="gauche") ?(n=n) ?(c=c) ?(h=h) y =
(*         Array.init (n+1) (fun i ->
                y.(i) +. (c /. ( sqrt(h**2. +. (( y.(i) -. y.(iMoins1 i) )**2.) ) ))
        ) *)

        match schema with
        | "gauche" ->
                
                (** Schéma gauche. *)

                Array.init (n+1) (fun i ->
                        y.(i) +. (c /. ( sqrt(h**2. +. (( y.(i) -. y.(iMoins1 i) )**2.) ) ))
                )
        | "droit" ->
                
                (** Schéma droit. *)

                Array.init (n+1) (fun i ->
                        y.(i) +. (c /. ( sqrt(h**2. +. (( y.(iPlus1 i) -. y.(i) )**2.) ) ))
                )
        | "centre" ->
                
                (** Schéma centré. @see <https://fr.wikipedia.org/wiki/Dérivée#D.C3.A9rivation_num.C3.A9rique> Cet article sur Wikipédia. *)

                Array.init (n+1) (fun i ->
                        y.(i) +. (c /. ( sqrt(h**2. +. 0.5*.(( y.(iPlus1 i) -. y.(iMoins1 i) )**2.) ) ))
                )
        | _ -> y;
;;

(** Stocke le résultat dans le tableau déjà existant (de taille n+1) ynew. @param schema précise le type de schéma (décentré gauche, droit, ou centré). Apparemment, seul le gauche est physiquement valide. *)

let nextY ?(schema="gauche") ?(n=n) ?(c=c) ?(h=h) yold ynew =
        for i = 0 to n do
                match schema with
                | "gauche" ->
                        
                        (** Schéma gauche. *)

                        ynew.(i) <- yold.(i) +. (c /. ( sqrt(h**2. +. (( yold.(i) -. yold.(iMoins1 i) )**2.) ) ))
                | "droit" ->
                        
                        (** Schéma droit. *)

                        ynew.(i) <- yold.(i) +. (c /. ( sqrt(h**2. +. (( yold.(iPlus1 i) -. yold.(i) )**2.) ) ))
                | "centré" ->
                        
                        (** Schéma centré. @see <https://fr.wikipedia.org/wiki/Dérivée#D.C3.A9rivation_num.C3.A9rique> Cet article sur Wikipédia. *)

                        ynew.(i) <- yold.(i) +. (c /. ( sqrt(h**2. +. 0.5*.(( yold.(iPlus1 i) -. yold.(iMoins1 i) )**2.) ) ))
                | _ -> ();
        done
;;


(**

Fonction principale de la simulation.

*)


(** Faire une pause pendant time seconde (par un appel Sys.command et pas Unix.sleep). *)

let sleep time = assert(0 == Sys.command (sprint "sleep %g" time));;

(** Petite optimisation : on affiche qu'une fois sur jskip. *)

let jskip = 5;;

(** Lance la simulation. @param schema précise le type de schéma (décentré gauche, droit, ou centré). Apparemment, seul le gauche est physiquement valide. @param inter affiche chaque étape si vrai. @param jmax durée maximum de la simulation (nombre d'étape, pas en temps effectivement calculé). *)

let simulation ?(schema="gauche") ?(inter=false) jmax () =
        let ym = ref (makeym ()) in
        print "Début de la simulation, pour n = %i points, h = %g, dt = %g, et k = %g.\n\n" n h dt k;
        for j = 0 to jmax do
                if 0 = (j mod jskip) then begin
                        print "\nSimulation, étape n° %i / %i.\n" j jmax;
                        print "  écriture du script dans 'plotdata.gnuplot'...\n  écriture des données dans 'data.txt'...\n  génération de l'image '%s/plot_%05d.png'." schema j;
                        print_and_launch ~schema:schema ~verb:true ~inter:inter ~j:j !ym ();
                end;
                ym := (newnextY ~schema:schema !ym);
                if inter then sleep 0.1;
                flush_all ();
        done
;;

(**

Lancement de la simulation.

*)


(* simulation ~schema:"droit" jmax ();; *)
(* simulation ~schema:"centre" jmax ();; *)
(* simulation ~schema:"gauche" jmax ();; *)

simulation ~schema:(try Sys.argv.(1) with _ -> "gauche") jmax ();;


(** Fin *)

print "\n\nFin de la simulation numérique.\nProgramme en OCaml v4+ par Lilian Besson (ENS Cachan), pour une collaboration avec Nicolas Vernier (Université Paris-Sud), à Mahindra École Centrale, (C) 2014.\n";;
print "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n";;