"Résolution de l''équation de la chaleur par des méthodes probabilistes"

"1-Saisie des paramètres"
N = 100; //Nombre de trajectoires
b = 10**-3; //Paramètre du drift

T = 60; //Durée de la simulation (en s)
Jt = 1500; //Nombre de subdivisions en temps
L = 3; //Demi-longueur de l'intervalle de résolution (en m): [-L,L]
Jx = 100; //Nombre de subdivisions en espace

"2-Calcul d''autres paramètres"
ht = T/Jt; //Pas de subdivision en espace
TT = 0:ht:T; //Intervalle de résolution en temps
hx = L/Jx; //Pas de subdivision en temps
XX = -L:hx:L; //Intervalle de résolution en espace

"3-Condition initiale"

function y = u0(x)
    y = exp(-x.^2);
endfunction

"4-Calcul des mouvements browniens"

B = zeros(N,Jt+1);

for jt=1:Jt
    B(:,jt+1) = B(:,jt) + grand(N,1,'nor',0,ht) + b*ones(N,1); //Accroissements des mouvements browniens + drift
end

"5-Calcul de la solution approchée"

function y = sol(jx,jt)
    //jx est un entier compris entre -Jx et Jx
    //jt est un entier compris entre 1 et Jt+1
    u = jx*hx*ones(N,1) + B(:,jt);
    y = (1/N)*sum(exp(-u.^2));
endfunction

U = zeros(2*Jx+1,Jt+1)

for jx = -Jx:Jx
    for jt = 1:Jt+1
        U(jx+Jx+1,jt) = sol(jx,jt);
    end    
end

"6-Tracé de la solution"

f0 = scf(1);
plot3d1(TT,XX,U',270,0,flag=[-6,2,4]);
f0.color_map = jetcolormap(512);
title('Equation de la chaleur- du transport avec '+string(N)+' particules')
xlabel('t [s]');
ylabel('x [m]');
colorbar(min(U),max(U));

xclick;

f1 = scf(2);
a=get("current_axes")
a.data_bounds=[-L min(U) ; L max(U)];
set(gca(),"grid",[1 1]);
title('Evolution de la densité de particules');
for n=1:Jt+1
    drawlater();
    clf();
    a=get("current_axes")
    a.data_bounds=[-L min(U) ; L max(U)];
    title('Evolution de la densité de partiules');
    set(gca(),"grid",[1 1]);
    plot(XX',U(:,n),'color','blue','thickness',2);
    drawnow();
end
