Diffusion d`un liquide dans un tuyau poreux

Projet de simulation MAP311
Lucas Gerin :
gerin@cmap.polytechnique.fr
Diffusion d'un liquide dans un tuyau poreux
Mots-clés : Graphes, Modélisation, Probabilités conditionnelles, ...
Le but de ce mini-projet est d'étudier dans des questions théoriques (T) et expérimentales
(E) la vitesse de diusion d'un liquide dans un tuyau poreux, que nous modéliserons par un
graphe très simple.
1 Construction et simulation du modèle
On considère le graphe inni N × {0, 1} en forme d'"échelle", c'est-à-dire que l'on met une
arête entre les points
(k, 1) ↔ (k + 1, 1) et (k, 0) ↔ (k + 1, 0) pour tout k ∈ N (les arêtes horizontales ),
(k, 0) ↔ (k, 1) pour tout k ∈ N (les arêtes verticales ).
On modélise de façon aléatoire la porosité du tuyau de la façon suivante. Soit p ∈ [0, 1], on se
donne une famille de variables aléatoires i.i.d. (Hk,δ )k∈Z,δ∈{0,1} telles que :
Hk,δ
(
2
=
1
avec probabilité p,
avec probabilité 1 − p.
La variable Hk,δ est interprétée comme le temps que met le liquide pour aller du point (k, δ) à
(k + 1, δ). On suppose pour simplier que la porosité des arêtes verticales n'est pas aléatoire : le
liquide met toujours un temps 1 pour aller de (k, δ) à (k, 1 − δ).
On suppose que l'on injecte en continu du liquide aux points (0, 0) et (0, 1), on cherche à
étudier le temps que met le liquide à atteindre le point (n, 0), pour n grand.
(0,1)
liquide
(k,1)
H(k,1)
1
(0,0)
(k,0)
H(k,0)
Pour k ∈ Z, δ ∈ {0, 1}, on note Tk,δ ∈ N le temps que met le liquide à atteindre (k, δ).
1. (T) On a par dénition T0,0 = T0,1 = 0 , justier que pour k ≥ 1
Tk,δ = min { Tk−1,δ + Hk−1,δ , Tk,1−δ + 1 } .
Solution: Lorsque le liquide arrive en (k, δ), soit il vient de (k −1, δ), soit de (k, 1−δ).
2. (E) En xant p, simuler et tracer plusieurs trajectoires de k 7→ Tk,0 .
2 Calcul de la vitesse asymptotique
On a bien sûr pour tout n, par dénition du modèle, n ≤ Tn,0 ≤ 2n, on s'attend à ce que
Tn,0 /n converge vers une constante, dans un sens à préciser.
Dans un premier temps, on se limite à la convergence de l'espérance. On cherche à calculer
et montrer l'existence de la limite
E[Tn,0 ]
µp = lim
.
n→+∞
n
Pour cela on pose, pour k ≥ 0, ∆k = Tk,0 − Tk,1 ∈ {−1, 0, 1}. Remarquons que ∆k est indépendante de Hk,0 et Hk,1 .
3. (T) Justier que pour tout k on a P(∆k = 1) = P(∆k = −1).
Solution: On a une symétrie parfaite du problème (condition initiale comprise) et du
coup :
loi
{Tk,δ }k,δ = {Tk,1−δ }k,δ .
En particulier
loi
Tk,0 − Tk,1 = Tk,1 − Tk,0 .
4. (T) Calculer P(∆k+1 = 1|∆k = 0) et P(∆k+1 = 0|∆k = 1), en déduire l'expression exacte
de P(∆k = 1). Conclure que pour tout ε ∈ {−1, 0, 1}
P(∆k = ε)
k→+∞
→
1/3.
Solution: Supposons que ∆k = 0. On xe Tk,1 = Tk,0 = t, les 4 cas sont dans le dessin
suivant :
t
t
2
2
t+2
t
t+2
t
1
2
t+1
t
t+2
t
2
1
t+2
t
t+1
t
1
1
t+1
t+1
On voit que P(∆k+1 = 1|∆k = 0) = p(1 − p) (deuxième cas ci-dessus).
Si on suppose que ∆k = 1, les cas sont :
t
t+1
2
2
t+2
t
t+3
t+1
1
2
t+1
t
t+2
t+1
2
1
t+2
t
t+2
t+1
1
1
t+1
t+2
et P(∆k+1 = 0|∆k = 1) = p(1 − p) (troisième cas).
On a ainsi
P(∆k+1 = 0) = P(∆k+1 = 1|∆k = 1)P(∆k = 1) + P(∆k+1 = 1|∆k = 0)P(∆k = 0)
xk+1
=
(1 − p(1 − p)) xk
+
p(1 − p)(1 − 2xk ),
où l'on a noté xk = P(∆k = 1) et utilisé la symétrie vue au-dessus : 2P(∆k = 1) +
P(∆k = 0) = 1. La suite (xk ) est donc une suite arithmético-géométrique, on la résout
facilement :
1 1
P(∆k = 1) = P(∆k = −1) = xk = − (1 − 3p(1 − p))k .
3 3
5. (T) Vérier que, pour tout k ,
Tk+1,0 = Tk,0 + Hk,0 − 1∆k =1,Hk,0 =2,Hk,1 =1 .
Solution: Dans les 8 cas dessinés à la question précédente, on peut vérier que l'on
a toujours Tk+1,0 = Tk,0 + Hk,0 , sauf sans le deuxième cas en bas. On a alors
Tk+1,0 = Tk,0 + Hk,0 − 1.
6. (T) En passant à l'espérance dans la formule précédente, démontrer que la limite µp existe
et que
2p + p2
µp = 1 +
.
3
Solution:
n−1
1X
1
E[Tn,0 ] =
E[Hk,0 ] − P ∆k = 1, Hk,0 = 2, Hk,1 = 1
n
n
=
1
n
k=0
n−1
X
1 + p − P(∆k = 1)p(1 − p)
k=0
n−1
= 1 + p − p(1 − p)
1X
n→+∞
P(∆k = 1) → 1 + p − p(1 − p)/3
n
k=0
en utilisant Césaro (rappelons que ∆k est indépendante de Hk,0 et Hk,1 ).
7. (E) On cherche à vérier expérimentalement la formule précédente. Pour chaque p =
1/10, 2/10, . . . , 9/10, simuler s fois la variable aléatoire Tn,0 /n (on prendra s de l'ordre
de quelques dizaines et n de l'ordre de quelques centaines). Tracer la courbe obtenue en
fonction de p et comparer avec la courbe théorique.
3 Étude expérimentale des uctuations
Il est possible de démontrer, mais cela dépasse le cadre de ce cours, que la suite Tn,0 , bien que
n'étant pas une somme de variables i.i.d., vérie une Loi des grands nombres et un Théorème
central limite. Précisément, il existe σp > 0 tel que
Tn,0 − nµp (loi)
√
→ N (0, σp2 ).
n
Nous allons vérier expérimentalement cette dernière armation.
8. (E) Fixer n (de l'ordre de la centaine) et simuler un grand nombre de fois la variable Tn,0 .
√
Tracer l'histogramme des valeurs de (Tn,0 − nµp )/ n obtenues par cette simulation et,
sur le même graphique, tracer la courbe correspondant à la gaussienne de même variance.
Interpréter graphiquement le résultat.
9. (T)-(E) Comment faire pour obtenir une estimation de σp ? (Pour cette question il n'est
pas demandé une preuve rigoureuse.) Pour p = 1/2, quelle estimation obtenez-vous ?
Solution: Notons Zn =
Tn,0 −nµp
√
.
n
On a (pour n grand)
Zn ≈ N (0, σp2 ).
On s'attend donc à ce que (c'est réellement le cas)
E[Zn2 ] → E[N (0, σp2 )2 ] = σp2 .
On simule donc pour n grand et plein de fois des variables aléatoires Zn2 et on estime σp2 par
la moyenne empirique.
(Avec
30000 simulations de T300 , je trouve σ1/2 ≈ 0.4196 puis σ1/2 ≈ 0.4171.)