//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////    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  
       //   [QQ Exemples Attraction   :  (beta,N,alpha) = 
       //                                 (1,10, 0) (1,10,10) (5,10,0) (5,10,1) (5,10,100)
       //                                 (1,10,0), (1,10,1), (1,10,10), (1,10,100) 
       //                                 (1,100,0), (1,100,1), (1,100,10), (1,100,100)
       //                                 (2,100,1), (5,100,1), (5,100,10), (5,100,100), ...]
       //    Attraction/Répulsion   :    (epsilon,beta,N,alpha) = 
       //                                 (1,1,10, 0) (-1,1,10,10) (1,1,10,10) (1,5,10,10) (-1,5,10,10) (-1,1,10,100) (1,1,10,100)
       //    Effet des températures       (+/-1,0.1, 10, 10 ou 100) (+1/-1, 5 ou 10, 10, 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;



for t=1:n
  
Mean_X=[Mean_X;mean(X(t,:))];
X=[X;X(t,:)-beta*dv(X(t,:))*(T/n)+epsilon*alpha*(mean(X(t,:))-X(t,:))*(T/n)+sigma*dW(t,:)];

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)