Diffusion moléculaire

In [ ]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib

Préambule

Question 1

Écrire le code de la fonction binomial.

Comme générateur aléatoire, utiliser la fonction random du module python random. Cette fonction retourne un flottant aléatoire entre 0 et 1.

In [ ]:
### À COMPLÉTER

def binomial(n, p):
    """ Simule une loi binomiale de paramètres n, p.
        n est le nombre d'épreuves de Bernoulli.
        p est la probabilité de succès d'une épreuve de Bernoulli.
        
        La fonction retourne un nombre de succès compris
        entre 0 et n.
    """
    ### À COMPLÉTER

[binomial(10, .1) for _ in range(10)]

Le problème

On considère un tube de longueur 30 cm initialement rempli d’eau. On dépose alors une petite quantité de sucre au milieu du tube et on s’intéresse à l’évolution des molécules de sucre dans le tube.

Nous nous proposons de modéliser cette situation par une marche aléatoire. Pour cela on discrétise l’axe longitudinal du tube selon $n$ points séparés par une distance $a$ et on suppose que toutes les molécules de sucre sont situés initialement sur le point central. On introduit par ailleurs un temps caractéristique $\tau$. La modélisation consiste à considèrer qu'à chaque intervalle de temps $\tau$ les molécules effectuent un saut sur un des deux points voisins. La direction du saut est aléatoire et indépendante pour chacune des molécules.

Les éléments de la simulation sont :

  • nmol : le nombre de molécules. Ce nombre devra être suffisamment grand pour limiter la fluctuation aléatoire.
  • n : le nombre de points
  • s : le nombre d'instants simulés
  • i0 : la position initiale des molécules
  • U : un tableau tel que U[k, i] représente le nombre de molécules au point $i$ à l'instant $k$.

Si on réalise la simulation telle que décrite ci-dessus, il est facile de voir qu’à chaque itération les positions occupées seront soit les positions paires soit les positions impaires. Pour éviter cela, nous nous proposons d'exécuter deux étapes à chaque itération, en passant par un tableau intermédiaire V. Le principe de la simulation est expliqué sur le schéma ci-dessous :

Question 2

Coder la fonction simul ci-dessous.

In [ ]:
def simul(n, s, i0, nmol):
    """ Réalise une simulation avec
             n : Nombre de points
             s : Nombre d'instants simulés
             i0 : Position initiale des molécules
             nmol : Nombre de molécules déposées
             
        La fonction retourne un tableau numpy U tel que
        U[k, i] soit le nombre de molécules à la position i
        au bout de k intervalles de temps.
    """
    ### À COMPLÉTER

Question 3

Tester la fonction simul sur les micros simulations suivantes :

U = simul(n=5, s=2, i0=2, nmol=1000000)
print(U[1])
U = simul(n=5, s=2, i0=0, nmol=1000000)
print(U[1])
U = simul(n=5, s=100, i0=0, nmol=1000000)
print(U[99])

Les proportions observées sont elles en adéquation avec le modèle ?

In [ ]:
 

La fonction binomial que vous avez écrit existe dans le module Numpy sous le nom nb.random.binomial. Pour des raisons de performance il est préférable maintenant d'utiliser cette fonction. Pour cela exécuter la ligne suivante :

In [ ]:
binomial = np.random.binomial

Question 4

Réaliser une simulation avec $n=100$ et $s=5000$. Un million de molécules sont déposés au milieu du tube.

Représenter la distribution des molécules dans le tube aux instants 0, 100, 200, 400, 800, 1600, 3200 et $s-1$.

Utilisez pour cela la fonction plt.bar.

In [ ]:
 

Question 5

Réaliser plusieurs simulations avec des couples $(n, s)$ tels que $\frac{n^2}{s}$ reste constant. Que peut-on observer ?

Choisissez par exemple un quotient $\frac{n^2}{s}$ égal à 0.05.

In [ ]:
 

Notons $D = \frac{a^2}{2\tau}$ et appelons ce rapport le coefficient de diffusion des molécules de sucre dans l'eau.

Relire le début du notebook pour la définition de $a$ et de $\tau$.

Question 6

Le coefficient de diffusion du sucre dans l'eau est $D=0.52\times 10^{-9}$ m²/s.

Réaliser une simulation à partir des données physiques de la façon suivante :

  1. Choisir le nombre de points $n$ de discrétisation du tube.
  2. En déduire la valeur de la valeur de $a$. On rappelle que le tube mesure 30cm. Se souvenir qu'une itération effectue deux marches aléatoires.
  3. En déduire la valeur de $\tau$ associée à $a$ en utilisant $D$.
  4. Calculer le nombre d'itérations $s$ nécessaires pour atteindre une durée $t_{\text{max}}$.
  5. Effectuer la simulation.

Évaluer le temps $t_{\text{max}}$ nécessaire pour qu'environ un millième des molécules atteignent les extrémités du tube.

In [ ]:
 

Question 7

Représenter la distribution de molécules à l'instant $t_{\text{max}}$ de la simulation précédente sur l'intervalle $[0,lg]$.

L'unité des abscisses est en mètre et celle des ordonnées en molécules par mètre.

Utiliser plt.bar avec l'argument optionnel width pour contrôler la largeur des bâtons.

In [ ]:
 

Modélisation par une EDP

Ce problème peut également être modélisé par l'équation aux dérivées partielles suivante : $$ \frac{\partial u}{\partial t}(t, x) = D\frac{\partial^2 u}{\partial x^2}(t, x)$$

où $u(t, x)$ est la concentration de molécules en $x$ à l'instant $t$.

Si $n_{\text{mol}}$ molécules sont déposées en $x_0$ à l'instant $t=0$, alors la solution de cette équation est

$$u(t, x) = \frac{n_{\text{mol}}}{\sqrt{4\pi Dt}}\exp\left(-\frac{(x-x_0)^2}{4Dt}\right)$$

Question 8

Représenter les deux solutions sur un même graphique pour vérifier qu'elles sont cohérentes entre elles.

In [ ]: