Générations aléatoires (programme 1)

Programme 1 Le premier programme à pour objectif de mettre en oeuvre quelques exemples de réalisations de variables aléatoires typiques de celles rencontrées dans la Physique du transport. On cherche à simuler et visualiser une marche aléatoire de diffusion multiple sur un plan à l'intérieur d'un carré de côté $a$. Les particules partent uniformément et de façon isotrope depuis l'un des côtés, vers l'intérieur du carré. Les particules sont suivies jusqu'à leur premier retour à la frontière. Le libre parcours de diffusion $\lambda_d$ est uniforme et la fonction de phase de diffusion simple est isotrope. On fera attention à la représentation de la distribution isotrope des particules entrantes.

Remarque: Les variables aléatoires rencontrées dans ce premier exemple sont décrites ci-après:

• Supposons que dans un repère orthonormé $(Oxy)$ le carré considéré soit celui correspondant à l'ensemble des points de coordonnées $(x,y)$ telles que $0<x<a$ et $0<y<a$. Choisissons alors le côté correspondant à $x=0$ et $0<y<a$. Chaque particule entre alors en $x=0$ mais sa coordonnée $y$ prend une valeur aléatoire comprise entre $0$ et $a$. Sa distribution est uniforme. Notons $Y$ la variable aléatoire correspondante. Sa densité de probabilité est $p_Y(y)=\frac{1}{a}$.
• La direction d'entrée de chaque particule est repérée à l'aide d'un angle $\theta$ compris entre $-\frac{\pi}{2}$ et $+\frac{\pi}{2}$ (angle trigonométrique usuel). L'hypothèse d'une distribution isotrope des particules entrantes se traduit, lorsqu'on se place à la frontière, par une distribution proportionnelle à $cos(\theta)$. Si on note $\Theta_e$ la variable aléatoire correspondante, sa densité de probabilité est donc (après normalisation) $p_{\Theta_e}(\theta)=\frac{cos(\theta)}{2}$.
• La distance $d$ parcourue par une particule avant la première diffusion est distribuée selon $exp\left(-\frac{d}{\lambda_d}\right)$. Si on note $D$ la variable aléatoire correspondante, sa densité de probabilité est donc (après normalisation) $p_D(d)=\frac{1}{\lambda_d}exp\left(-\frac{d}{\lambda_d}\right)$.
• Lors de chaque diffusion, la direction de propagation après diffusion est repérée à l'aide d'un angle $\theta$ compris entre $0$ et $2\pi$ (angle trigonométrique usuel). Comme la fonction de phase de diffusion simple est isotrope, cet angle est distribué uniformément sur $[0,2\pi]$. Si on note $\Theta_d$ la variable aléatoire correspondante, sa densité de probabilité est donc $p_{\Theta_e}(\theta)=\frac{1}{2\pi}$.

Exemple de solution (programmé en fortran) :

c---------------------------------------------------------------------
c Une trajectoire de multi-diffusion depuis le bord gauche d'un carre
c---------------------------------------------------------------------
      PROGRAM prog1
      IMPLICIT NONE
      INTEGER iseed
      COMMON/comsee/iseed
      DOUBLE PRECISION pi
      PARAMETER(pi=3.141592654)
c---------------------------------------------------------------------


      DOUBLE PRECISION cote,lambda
      DOUBLE PRECISION x,y,theta
      DOUBLE PRECISION r,xold,yold,d,dsecx,dsecy


c---------------------------------------------------------------------
      OPEN(10,FILE='sorties.out')
      PRINT *,'germe du generateur aleatoire : '
      READ *,iseed
c---------------------------------------------------------------------


      cote=0.3
      lambda=0.1

      x=0.
      CALL rand_uniforme(r)
      y=r*cote
      CALL rand_uniforme(r)
      theta=asin(2.*r-1.)
      
 1    CONTINUE

      xold=x
      yold=y
      WRITE(10,*) x,y

      CALL rand_uniforme(r)
      d=-lambda*log(1-r)
      x=x+d*cos(theta)
      y=y+d*sin(theta)

      IF ((x.GE.0).AND.(x.LE.cote).AND.(y.GE.0).AND.(y.LE.cote)) THEN
         CALL rand_uniforme(r)
         theta=r*2*pi
         GOTO 1
      ELSE
         IF (cos(theta).GE.0.) THEN
            dsecx=(cote-xold)/cos(theta)
         ELSE
            dsecx=(0.-xold)/cos(theta)
         ENDIF
         IF (sin(theta).GE.0.) THEN
            dsecy=(cote-yold)/sin(theta)
         ELSE
            dsecy=(0.-yold)/sin(theta)
         ENDIF
         IF (dsecx.LE.dsecy) THEN
            d=dsecx
         ELSE
            d=dsecy
         ENDIF
         x=xold+d*cos(theta)
         y=yold+d*sin(theta)
         WRITE(10,*) x,y
      ENDIF


c---------------------------------------------------------------------
      CLOSE(10)
      STOP
      END
c---------------------------------------------------------------------