//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////    Gradient stochastique avec intéraction : Equation de Langevin avec intéraction champ moyen                        ////
////////    Modèles de S. Herrmann et J. Tugaut : interactions (attraction/répulsion) autour des moyennes                     ////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


clf();
stacksize("max");


T=10;
n=1000; //Pas de temps T/n
N=100;  //Nombre de particules (N T/n<1/10); [Rmq : N=1 = discrétisation de Langevin]
alpha=-10; // taux d'interaction autour de la moyenne des particules (
           // alpha=0 => N iid gradients stoch.)           
epsilon=1; //Attraction +1 répulsion -1
sigma=1; // Ecart type des perturbations



beta=2; //Température de refroidissement : Choix max =>  beta * (T/n) < 5/100  
       


//Incréments browniens n pas de temps sur [0,T]
dW=grand(n,N,"nor",0,sqrt(T/n));
time_mesh=(T/n)*[0:n]';


// Fonction énergie double potentiel symétrique

//well=2;//Positionnement des puits

a=-1;
b=5; // choisir a<=b
maxi=5;
c=maxi*((b-a)/2)^(-4); //maxi entre les puits 

function y=v(x)
//y=(x^2-well^2)^2/(well^2);
y=c*(((x-a)^2).*((x-b)^2))
endfunction

//radient de la fonction d'énergie

function y=dv(x)
//y=4*x^3/(well^3)-4*x/(well^2);
y=(4*c)*((x-a).*(x-b).*(x-((a+b)/2)))
endfunction


//Initialisation des processus et de sa moyenne empirique
X=zeros(1,N);
Mean_X=0;



////////////////////////////////////////////////////////////////////////////
////// COMPLETER la ligne rouge :  EVOLUTION DU SCHEMA d'EULER          ////
////////////////////////////////////////////////////////////////////////////

for t=1:n
  
Mean_X=[Mean_X;mean(X(t,:))];
X=[X;	  ];

end


//AFFICHAGE DES RESULTATS 

clf(0);xset("window",0);show_window();

for NN=1:N
 plot2d(time_mesh,X(:,NN),NN)
end


clf(1);xset("window",1);show_window();

 plot2d(time_mesh,(Mean_X)',4)



clf(2);xset("window",2);show_window(); //tracé de v entre x_0< (a+b)/2< x1 pour v(x_0)=v(x_1)=v((a+b)/2)=maxi

u=linspace(((a+b)/2)-sqrt(((a+b)/2)^2+((b-a)/2)^2-a*b),((a+b)/2)+sqrt(((a+b)/2)^2+((b-a)/2)^2-a*b));
plot2d(u,v(u),5)