let k = 1.0;;
let x0 = 0. and x1 = 1.;;
let lx = x1 -. x0;;
let n = 1500;;
let h = lx /. (float n);;
let x_i i = x0 +. ((float i) *. h);;
let dt = 0.0001;;
let t_j j = ((float j) *. dt);;
let jmax = 3000;;
let print = Format.printf;;
let sprint = Format.sprintf;;
let fprint = Printf.fprintf;;
let dftlext = "svg";;
let dftlext = "png";;
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;
;;
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;
;;
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);;
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);;
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 ());
;;
let beta = 0.1
let beta = 1.0
let y_0 x = (1. -. x) *. beta ;;
let y_0 x = (1. -. (x ** 2.)) *. beta ;;
let makeym () = Array.init (n+1) (fun i -> y_0 (x_i i));;
let c = h *. dt *. sqrt(k);;
let iMoins1 = function 0 -> 1 | i -> i-1;;
let iPlus1 ?(n=n) = function i when i=n -> i | i -> i+1;;
let newnextY ?(schema="gauche") ?(n=n) ?(c=c) ?(h=h) y =
match schema with
| "gauche" ->
Array.init (n+1) (fun i ->
y.(i) +. (c /. ( sqrt(h**2. +. (( y.(i) -. y.(iMoins1 i) )**2.) ) ))
)
| "droit" ->
Array.init (n+1) (fun i ->
y.(i) +. (c /. ( sqrt(h**2. +. (( y.(iPlus1 i) -. y.(i) )**2.) ) ))
)
| "centre" ->
Array.init (n+1) (fun i ->
y.(i) +. (c /. ( sqrt(h**2. +. 0.5*.(( y.(iPlus1 i) -. y.(iMoins1 i) )**2.) ) ))
)
| _ -> y;
;;
let nextY ?(schema="gauche") ?(n=n) ?(c=c) ?(h=h) yold ynew =
for i = 0 to n do
match schema with
| "gauche" ->
ynew.(i) <- yold.(i) +. (c /. ( sqrt(h**2. +. (( yold.(i) -. yold.(iMoins1 i) )**2.) ) ))
| "droit" ->
ynew.(i) <- yold.(i) +. (c /. ( sqrt(h**2. +. (( yold.(iPlus1 i) -. yold.(i) )**2.) ) ))
| "centré" ->
ynew.(i) <- yold.(i) +. (c /. ( sqrt(h**2. +. 0.5*.(( yold.(iPlus1 i) -. yold.(iMoins1 i) )**2.) ) ))
| _ -> ();
done
;;
let sleep time = assert(0 == Sys.command (sprint "sleep %g" time));;
let jskip = 5;;
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
;;
simulation ~schema:(try Sys.argv.(1) with _ -> "gauche") jmax ();;
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";;