/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////                 Interprétation probabiliste des équations de Burgers (1d)                     //////////////
////////                 Modèles de McKean  : interactions avec les fonctions de répartition           ///////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


clf();
stacksize("max");


T=1;   // Modèle sur [0,T]
n=100; //Pas de temps T/n
s=.5;   // variance de la diffusion : s=sigma^2
N=100;  //Nombre de particules (N T/n<1/10), choisir N>=10.
N_MC=100; 

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





//Fonction Heavyside h(x)=1_{x>0}

function y=h(x)
H=(x>=0)    
y=double(H)
endfunction



//Solution Burgers 1d pour V(0,x)=1-h(x)



function y=VS(t,x)
e1=erfc(-x/sqrt(2*s*t));
e2=exp((t-2*x)/(2*s));
e3=erfc((t-x)/sqrt(2*s*t));
y=1-e1/(e1+e2*(2-e3));
endfunction


//Approximation empirique de V(t,x)=1-P(X_t<x)=P(X_t>x)

function y=VN(x)
    size_x=size(x);
    v=[];
    for i=1:size_x(1,2)
        u=x(i)-x;
        Hu=(u>=0);
        hu=double(Hu)
        v=[v,1/N*sum(hu)];
    end
    y=1-v;
endfunction

//Approximation Monte Carlo W^i N_MC iid N(moyenne = x, sigma^2= t)

function y=VMC(t,x,N_MC)
    Wx=grand(1,N_MC,"nor",x,sqrt(s*t));
    A=sum((double((Wx<0)).*(exp(-Wx/s))))/N_MC;
    B=1+sum((double((Wx<0)).*(exp(-Wx/s)-1)))/N_MC;
    y=A/B;
endfunction

//Initialisation des N processus 

X=zeros(1,N);
V=[];
VV=[];

//Boucle d'évolution du processus de diffusion non linéaire
// (N particules en interaction = transitions gaussiennes avec interaction au niveau des moyennes)

for t=1:n
 
    
// Calcul empirique de la dérive b(x,p)=1-int H(x-y) p_t(y) dy des particules 
// en remplaçant  p_t(y) dy par la mesure d'occupation du système de N particules
   
Y=[];  
for i=1:N
U=X(t,i)-X(t,:);
Y=[Y,1-(1/N)*sum(h(U))]; 
end

 
  
//Evolution des N particules sur un pas de temps

X=[X;X(t,:)+Y.*(T/n)+dW(t,:)];

//Calcul empirique de la solution de l'équation de Burgers

V=[V;VN(X(t,:))];

//Valeurs exactes de la solution de l'équation de Burgers sur les 
// particules simulées

Z=[];
for i=1:N
Z=[Z,VS(t/n,X(t,i))]
end
VV=[VV;Z];

end




//Tracé des 10 premières trajectoires de la diffusion non linéaire

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

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



clf(1);xset("window",1);show_window();
l=[]; 
for k=-2:0.01:4
l=[l,VS(t/n,k)];  // Valeurs exactes de la solution sur une grille
end
subplot(221)
plot(-2:0.01:4,l,'red')
 plot2d(X(n,:),V(n,:),style=0,rect=[-2,-.1, 4, 1.1]); // Approximation empirique
 subplot(222)
plot(-2:0.01:4,l,'red')
plot2d(X(n,:),VV(n,:),style=0,rect=[-2,-.1, 4, 1.1]); //Valeurs exactes sur les points simulés

 subplot(223)
ll=[]; 
for k=-2:0.01:4
ll=[ll,VMC(t/n,k,N_MC)];  // Approximation Monte Carlo de la solution 
end
plot(-2:0.01:4,l,'red');
plot2d(-2:0.01:4,ll,style=0,rect=[-2,-.1, 4, 1.1]);