suivant: Sensibilités (programmes 8 et monter: Travaux pratiques concernant la précédent: Echantillonnage préférentiel (programmes 4
Programme 6 Dans le cas de la diffusion pure, estimer la fraction des particules entrant par un côté qui sortent du carré dans un petit segment du coté opposé. On observera que l'on rencontre des difficultés de convergences lorsque le segment est de petite taille (imprécision de l'estimation malgré un nombre élevé de trajectoires simulées).
Exemple de solution (programmé en fortran) :
c---------------------------------------------------------------------
c Calcul de la fraction sortant par le côté d'en face entre y1 et y2
c---------------------------------------------------------------------
PROGRAM prog5
IMPLICIT NONE
INTEGER iseed
COMMON/comsee/iseed
DOUBLE PRECISION pi
PARAMETER(pi=3.141592654)
c---------------------------------------------------------------------
INTEGER itir,ntir
DOUBLE PRECISION w,som_w,som_w2,moy_w,moy_w2,sig_w,sig_moy_w
DOUBLE PRECISION cote,lambda,y1,y2
DOUBLE PRECISION x,y,theta,L,sortie_y1y2
DOUBLE PRECISION r,xold,yold,Lold,d,dsecx,dsecy
DOUBLE PRECISION moy_L,sig_moy_L
DOUBLE PRECISION frac_y1y2,sig_frac_y1y2
c---------------------------------------------------------------------
OPEN(10,FILE='sorties.out')
PRINT *,'germe du generateur aleatoire : '
READ *,iseed
c---------------------------------------------------------------------
cote=0.3
lambda=0.01
y1=0.1495
y2=0.1505
ntir=100000
som_w=0.
som_w2=0.
DO itir=1,ntir
L=0.
x=0.
CALL rand_uniforme(r)
y=r*cote
CALL rand_uniforme(r)
theta=asin(2.*r-1.)
1 CONTINUE
xold=x
yold=y
Lold=L
CALL rand_uniforme(r)
d=-lambda*log(1-r)
x=x+d*cos(theta)
y=y+d*sin(theta)
L=L+d
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)
L=Lold+d
IF ((cos(theta).GE.0.).AND.(dsecx.LE.dsecy)
& .AND.(y.GE.y1).AND.(y.LE.y2)) THEN
sortie_y1y2=1.
ELSE
sortie_y1y2=0.
ENDIF
ENDIF
WRITE(10,*) itir,sortie_y1y2
w=sortie_y1y2
som_w=som_w+w
som_w2=som_w2+w**2
ENDDO
moy_w=som_w/ntir
moy_w2=som_w2/ntir
sig_w=sqrt(moy_w2-moy_w**2)
sig_moy_w=sig_w/sqrt(dble(ntir))
frac_y1y2=moy_w
sig_frac_y1y2=sig_moy_w
PRINT *,frac_y1y2,sig_frac_y1y2
c---------------------------------------------------------------------
CLOSE(10)
STOP
END
c---------------------------------------------------------------------
Exemple d'extension en géométrie tridimensionnelle (sous EDStar) : en cours.
Programme 7 Reprendre l'exemple précédent en utilisant un algorithme réciproque (de façon à résoudre le problème de convergence).
Exemple de solution (programmé en fortran) :
c---------------------------------------------------------------------
c Calcul de la fraction sortant par le côté d'en face entre y1 et y2
c par un algorithme réciproque
c---------------------------------------------------------------------
PROGRAM prog6
IMPLICIT NONE
INTEGER iseed
COMMON/comsee/iseed
DOUBLE PRECISION pi
PARAMETER(pi=3.141592654)
c---------------------------------------------------------------------
INTEGER itir,ntir
DOUBLE PRECISION w,som_w,som_w2,moy_w,moy_w2,sig_w,sig_moy_w
DOUBLE PRECISION cote,lambda,y1,y2
DOUBLE PRECISION x,y,theta,L,sortie_y1y2
DOUBLE PRECISION r,xold,yold,Lold,d,dsecx,dsecy
DOUBLE PRECISION moy_L,sig_moy_L
DOUBLE PRECISION frac_y1y2,sig_frac_y1y2
c---------------------------------------------------------------------
OPEN(10,FILE='sorties.out')
PRINT *,'germe du generateur aleatoire : '
READ *,iseed
c---------------------------------------------------------------------
cote=0.3
lambda=0.01
y1=0.1495
y2=0.1505
ntir=100000
som_w=0.
som_w2=0.
DO itir=1,ntir
L=0.
x=cote
CALL rand_uniforme(r)
y=y1+r*(y2-y1)
CALL rand_uniforme(r)
theta=asin(2.*r-1.)+pi
1 CONTINUE
xold=x
yold=y
Lold=L
CALL rand_uniforme(r)
d=-lambda*log(1-r)
x=x+d*cos(theta)
y=y+d*sin(theta)
L=L+d
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)
L=Lold+d
IF ((cos(theta).LE.0.).AND.(dsecx.LE.dsecy)) THEN
sortie_y1y2=1.
ELSE
sortie_y1y2=0.
ENDIF
ENDIF
WRITE(10,*) itir,sortie_y1y2
w=(y2-y1)/cote*sortie_y1y2
som_w=som_w+w
som_w2=som_w2+w**2
ENDDO
moy_w=som_w/ntir
moy_w2=som_w2/ntir
sig_w=sqrt(moy_w2-moy_w**2)
sig_moy_w=sig_w/sqrt(dble(ntir))
frac_y1y2=moy_w
sig_frac_y1y2=sig_moy_w
PRINT *,frac_y1y2,sig_frac_y1y2
c---------------------------------------------------------------------
CLOSE(10)
STOP
END
c---------------------------------------------------------------------
Exemple d'extension en géométrie tridimensionnelle (sous EDStar) : en cours.