TangenteX.com

Quelques rappels de physique statistique

Nous avons déjà rencontré l'algorithme de Metropolis dans plusieurs pages de TangenteX: pour simuler la transition de phase ferromagnétique-paramagnétique autour de la température de Curie d'un matériau avec le modèle d'Ising, ou encore pour résoudre le problème du voyageur de commerce à l'aide de l'algorithme de recuit simulé.

Aujourd'hui, je vais revenir sur cet algorithme dans son contexte natif, c'est à dire la physique statistique. Cet algorithme est tellement utile et efficace dans l'emploi des méthodes de Monte-Carlo qu'on le retrouve dans des domaines très variés de simulation comme la finance ou la dynamique des populations. C'est un bon exemple d'essaimage des méthodes de la physique numérique aux autres sciences.

Faisons d'abord quelques rappels de physique statistiques, indispensables pour aborder l'algorithme de Metropolis.

Les différents niveaux d'analyse d'un système réel

Les physiciens définissent trois échelles pour étudier un système :

Pourquoi une "physique statistique"

Nous savons que la température d'un corps macroscopique, disons l'eau dans un verre pour fixer les idées, dépend de l'agitation thermique des molécules d'eau dans le verre.

Cette agitation thermique se modélise par le déplacement des molécules, c'est à dire par la connaissance, à l'échelle microscopique, de la position et de la vitesse de chacune de ces molécules. Restons dans l'approximation non quantique et demandons nous s'il serait possible, à partir d'informations locales microscopiques, de déterminer la température de l'eau dans notre verre.

Considérons donc notre verre d'eau, qui contient 10 cl d'eau, soit 100 grammes d'eau. 100 grammes, c'est à peu près 6 moles d'eau soit 6 fois 6.02 1023 molécules d'eau, environ 4.1024 molécules d'eau. En mécanique classique, le calcul d'une trajectoire dans l'espace de notre verre nécessite 6 informations pour chaque molécule (3 coordonnées, 3 vitesses). Pour calculer la trajectoire de l'ensemble des molécules du verre, il faudrait résoudre plus de 2.1025 équations différentielles de deuxième ordre !

Et encore, ce serait en supposant que nous puissions déterminer une trajectoire ! Car même sans entrer dans des considérations quantiques, les molécules se heurtent un nombre de fois tellement grand (de l'ordre de 108 fois par seconde, selon la température) que le calcul n'aurait pas vraiment de sens.

Tout ça pour dire qu'il est absolument inenvisageable de calculer la température de l'eau de notre verre à partir des comportements microscopiques, à supposer d'ailleurs que nous sachions les mesurer.

C'est là qu'intervient la physique statistique. Comme son nom le laisse deviner, elle consiste à déterminer des valeurs moyennes macroscopiques comme la température, à partir de l'étude des comportements microscopiques.

Passer du microscopique au macroscopique

Pour illustrer cette notion, je vais prendre l'exemple de la température, grandeur observable connue de tous. La température est une grandeur macroscopique : définir la température d'une molécule d'eau n'aurait aucun sens. Oui, mais pourquoi ?

Revenons à notre verre d'eau, posé sur mon bureau, à une température de 293 K. Je suppose que l'eau contenue dans le verre est à l'équilibre thermodynamique. Je vais m'attarder un instant sur cette notion d'équilibre thermodynamique. Lorsque je sors l'eau du frigo et que je la verse dans mon verre posé sur mon bureau, les échanges de chaleur entre l'eau à 280 K et le milieu ambiant à 293 K sont importants. Le milieu ambiant se comporte comme un thermostat qui va imposer sa température à l'eau. La température de l'eau augmente parce qu'elle reçoit de l'énergie du milieu alors que la température du milieu, elle, ne bouge pratiquement pas. La température de l'eau se rapproche de la température du milieu : on dit que l'eau se "thermalise". Toutefois, même lorsque les deux températures sont pratiquement identiques, les échanges d'énergie entre les deux milieux se poursuivent de façon aléatoire et tout à fait négligeables sur le plan macroscopique (de l'ordre de \( \dfrac{1}{\sqrt{N}}\)).

Du point de vue macroscopique, la vitesse moyenne <v> de l'ensemble de ses 1024 molécules s'exprime par \( <v> = \displaystyle  \dfrac{1}{N} \sum_{i=1}^N v_i \), où vi est la vitesse individuelle de chaque molécule d'eau, N étant bien sur égal à 1024.

J'en déduit l'énergie cinétique moyenne de l'ensemble des molécules d'eau par la formule classique  \( <E_c> = \dfrac{1}{2}m<v>^2 \) soit \( <E_c> = \displaystyle  \dfrac{1}{N} \sum_{i=1}^N \dfrac{1}{2}m v_i^2 \).

D'autre part, la thermodynamique nous fournit la relation entre <Ec> et la température T  \( <E_c> = \displaystyle  \dfrac{3}{2} k_B T \), valable à l'équilibre thermique avec kB la constante de Boltzmann (kB = 1,38.10-23 m2 kg s-2 K-1).

Nous savons donc statistiquement établir une relation entre une grandeur macroscopique observable T et des grandeurs microscopiques vi :

\( \displaystyle  \dfrac{3}{2} k_B T  = \dfrac{1}{N} \sum_{i=1}^N \dfrac{1}{2}m v_i^2 \)

Cette relation fait intervenir la valeur moyenne des vitesses des molécules d'eau. Définir la température d'une molécule n'aurait donc pas de sens avec cette définition. Le passage du microscopique au macroscopique implique la mise en jeu d'un grand nombre d'objets microscopiques.

Micro-état et distribution de Boltzmann

Considérons maintenant un système de N particules identiques et discernables à spectre d'énergie discret, en équilibre thermodynamique à la température T. Que signifie "à spectre d'énergie discret" ? C'est l'expression consacrée pour dire que chaque particule peut prendre une valeur entière d'énergie comprise entre 0 et p, p pouvant être infini. C'est le genre de système que l'on rencontre en physique quantique mais aussi lorsqu'on veut simuler un système en physique statistique.

On appelle 'micro-état' l'état de chacune des particules, défini ici par son énergie.

Comme notre système est "thermalisé", la distribution du niveau d'énergie de chaque particule est réparti selon la loi de distribution de Boltzmann. On exprime ce fait en définissant la probabilité qu'une particule occupe un niveau d'énergie Ei donné, c'est à dire la probabilité de son microétat, par :

\(  \displaystyle p(E_i) = \dfrac{1}{Z} e^{-\dfrac{E_i}{k_B T}} \)

Z est la fonction de partition de la particule et vaut \( Z = \displaystyle \sum_{i=1}^N e^{-\dfrac{E_i}{k_B T}} \). Cette fonction est définie pour normaliser à 1 la distribution de probabilités.

On remarque que la probabilité d'un micro-état est d'autant plus élevée que l'énergie de ce micro-état est plus faible.

SI j'appelle <E> l'énergie moyenne d'une particule de mon système, alors on démontre que :

\( \displaystyle \dfrac{<E_{totale}>}{E_{totale}} = \dfrac{1}{\sqrt{N}} \dfrac{<E>}{E} \)

La variation d'énergie totale du système Etotale est proportionnelle à la variation d'énergie moyenne d'une particule, le facteur de proportionnalité étant \(\dfrac{1}{\sqrt{N}} \), soit quelque chose de très très petit. Quand N devient très grand, de l'ordre du nombre d'Avogadro, c'est à dire lorsqu'on passe à l'échelle macroscopique, la variation d'énergie du système tend vers 0. Au niveau macroscopique, on peut définir une énergie interne bien déterminée, caractéristique d'un système thermalisé.

Le principe de l'algorithme de Metropolis

Le contexte

Imaginez que l'on vous demande de simuler la thermalisation d'un système de 100 particules identiques et discernables, dont chacune peut prendre 10 niveaux d'énergie différents. Cela représente 10100 micro-états différents. Pour simuler la thermalisation, vous allez être obligé de calculer l'évolution de la probabilité de chaque micro-état avec la petite formule que nous avons vu ci-dessus. Une sommation sur 10100 micro-états ... Rien que pour calculer la fonction de partition ! Dans l'état actuel de nos outils informatiques, c'est inimaginable, et pour longtemps encore j'imagine. Et pourtant un système de 100 particules c'est minuscule !

Devant une intégration (ou une sommation) plutôt "tordue", nous avons déjà abordé une méthode : la méthode de Monte-Carlo. Elle est basée sur un mécanisme assez simple dans le principe : on tire une très grande quantité de nombres aléatoires qui représenteraient des micro-états puis on calcule la valeur moyenne de la grandeur recherchée, ici l'énergie. L'idée est bonne sauf qu'il y a un petit problème : la loi de distribution utilisée est la loi uniforme. Tous les micro-états auraient la même probabilité, ce qui est contraire à la physique du système comme on l'a constaté plus haut. Il faudrait adapter la loi de distribution pour coller à la loi de distribution de Boltzmann (ou loi canonique ou Maxwell-Boltzmann).

A cette contrainte, il faut ajouter deux autres contraintes de nature à la fois mathématique et physique :

La démarche

Il s'agit tout d'abord de construire une chaîne de Markov ergodique qui parcourt l'ensemble de l'espace des états de notre système. Notre chaîne de Markov est définie par la probabilité de transition du micro-état \( \mu_i \) au micro-état \( \mu_{i+1} \). Cela nous donne une distribution de probabilités des micro-états générés par la chaîne de Markov, que je vais appeler \( S(\mu , i) \). Ce que nous attendons de cette distribution, c'est qu'elle converge vers la distribution de Boltzmann pour simuler la thermalisation de notre système.

L'algorithme

Je donne ci-dessous le principe général de l'algorithme de Métropolis pour la distribution de Boltzmann, qui devra être adapté pour chaque cas, comme je l'ai fait dans la simulation du modèle d'Ising :

Dans le chapitre suivant, nous allons implémenter cet algorithme en Python à travers un exemple simple de simulation en physique statistique.

De l'intérêt d'un bon générateur de nombres aléatoires

L'algorithme de Metropolis, comme d'ailleurs toutes les variantes de la méthode de Monte-Carlo, utilise un générateur de nombres aléatoires pour générer des distributions de variables aléatoires.

Nous savons que les générateurs utilisés en informatique sont nécessairement périodiques (notre espace d'adresses n'est pas infini), ce qui peut poser problème. L'algorithme de Metropolis appliqué à un réseau monodimensionnel utilise 3 tirages aléatoires par itération. Nous verrons dans la suite que le nombre d'itérations peut atteindre et même dépasser le million. Il faut donc s'assurer que la période du générateur soit très supérieure à ce nombre de tirages.

D'après la doc Python, le générateur de nombres aléatoires de Scipy présente une période de 219937-1. Cela devrait donc nous aller !

Simuler la thermalisation d'une population de particules

Le modèle

Imaginons une chaîne de N particules. Ces particules peuvent prendre un niveau d'énergie entier (spectre d'énergie discret) compris entre 0 et Q, Q grand mais fini. Je noterai ei l'énergie d'une particule. L'énergie totale du système ET est égale à \( \displaystyle E_T = \sum_{i=0}^{N-1} e_i \).

La simulation que je vous propose consiste à examiner la thermalisation de notre système unidimensionnel de particules lorsqu'on le met au contact d'un thermostat qui lui impose une température T.

Au temps t=0, chaque particule possède un état aléatoire, dont l'énergie variera entre 0 et 2, 2 est purement arbitraire, et son énergie sera donc non nulle.

Pendant la thermalisation, c'est à dire l'évolution vers l'équilibre avec le thermostat, les particules échangent de l'énergie entre elles et avec le thermostat par de multiples collisions. A l'équilibre, la distribution de l'énergie dans les particules aura considérablement évolué pour atteindre une distribution de Boltzmann.

Je vais appliquer l'algorithme de Metropolis à ce système pour suivre sa thermalisation, c'est à dire la convergence de ET vers une valeur d'équilibre.

Le code

Commençons par définir les paramètres de la simulation, la taille de notre système et le nombre d'itérations, ce qui simule le temps dans le processus de thermalisation :
N = 1000                    # taille du système
Iterations = 800000         # nombre d'itérations

Au cours de la manip, je ferai varier le nombre N de particules du système pour observer l'influence de ce paramètre sur la thermalisation. Le nombre d'itérations doit être suffisamment grand pour observer la convergence vers l'équilibre. Il dépend de la température du thermostat et de N. Je l'ai choisi par essais successifs.

Définissons maintenant les caractéristiques du réseau :

T = 10.0                              # température initiale du réseau
beta = 1.0/T                          # coefficient de Boltzmann simplifié
Reseau, E0 = InitReseauAleatoire(N)   # initialisation du réseau
Energie = [E0]                        # énergie initiale du réseau

La température du thermostat est assez arbitraire. Je la fixe à 10.0 (sans unité) de façon assez arbitraire. Vous pourrez essayer d'autres valeurs.

La fonction InitReseauAleatoire() permet de créer une matrice aléatoire dont l'énergie de chaque particule varie de façon entière entre 0 et 2 (là aussi, 2 est arbitraire). Le code de cette fonction est très simple :

def InitReseauAleatoire(n):
    # construction du réseau avec des états aléatoires variant entre 0 et 2   
    res = random_integers(0, 2, n)
    # calcul de l'énergie initiale
    Ei = 0   
    for i in xrange(n):
        Ei += res[i]
    return res,Ei

La fonction retroune le réseau initialisation et son énergie initiale.

Enfin, la principale pièce de ce code, l'implémentation de l'algorithme de Metropolis, adapté au besoin de notre simulation :

def Metropolis(res,beta,En):
    n = len(res)   
    # choix aléatoire de la particule à traiter
    j = randint(0,n-1)
    # je tire au hasard une variation d'énergie discrète +1 ou -1
    if res[j] == 0:
        dE = 1
    else:
        if random() < 0.5:
            dE = 1
        else:
            dE = -1   
    # si la variation d'énergie est négative, je l'accepte
    if dE == -1:
        res[j] += dE
        En += dE
   # sinon je l'accepte si la proba est inférieur à exp(-beta*dE)   
    else:
        if random() < exp(-dE*beta):
            res[j] += dE
            En += dE           
    return[res, En]

Les commentaires dans le code sont assez parlant je crois. La fonction Metropolis() retourne le nouvel état du réseau et son énergie actuelle.

Il ne reste plus qu'à définir la boucle qui parcoure la chaîne de Markov et qui constitue le squelette de la méthode de Monte-Carlo :

t0 = time.time()
for i in xrange(Iterations):
    En = Energie[i]
    Reseau,En = Metropolis(Reseau,beta,En)
    Energie.append(En)
print "temps de calcul : %3.4f secondes\n" % (time.time() - t0)

Pour avoir une idée du temps de traitement, j'ai ajouté deux instructions qui permettent de mesurer et d'afficher celui-ci (entre 7 et 9 secondes sur mon Mac).

Le code source de ce script AlgoMetropolis.py est disponible ici.

Les résultats

Observation de l'influence de N sur la simulation

La méthode de Monte-Carlo et l'algorithme de Metropolis sont des méthodes statistiques. La taille du système, fixée par le paramètre N dans mon code, doit donc influer sur la qualité des résultats.

Je vais donc faire plusieurs essais avec N= 10, N = 100 et N = 1000 sur 8.105 itérations, en traçant la courbe d'évolution de l'énergie totale du système en fonction du nombre d'itérations, qui simule le processus de thermalisation en fonction du temps.

La température du thermostat est fixée à 10.0.

Pour  N = 10, j'obtiens :

Thermalisation avec 10 particules

Les oscillations de l'énergie du système simulé sont considérables. Quelque soit le nombre d'itérations, aucune convergence vers une valeur d'équilibre n'est constatée. La population de mon système est bien trop petite pour que les lois statistiques soient applicables et efficaces.

Pour N = 100, j'obtiens :

Thermalisation avec 100 particules

On observe une tendance à converger vers une valeur moyenne de l'énergie du système simulé. Toutefois, les oscillations autour de cette valeur moyenne sont élevées et ne permettent pas de calculer cette valeur moyenne. La population de mon système est encore trop petite pour que l'algorithme de Metropolis soit efficace.

Pour N = 1000, j'obtiens :

Thermalisation avec 1000 particules

L'amplitude des oscillations a très nettement diminué. L'effectif de la population du système devient compatible avec une bonne efficacité de l'algorithme de Metropolis. D'ailleurs, vous pourrez essayer de faire tourner le script avec une population plus grande et une nombre d'itérations également plus grand pour observer la qualité de la convergence.

Lorsque l'équilibre vous semble atteint, vous pourrez également calculer la valeur moyenne de l'énergie du système et vérifier la loi de Boltzmann.

Observer l'état du système avant et après la thermalisation

J'ai initialisé le réseau avec des particules dont l'énergie était fixée aléatoirement avec une distribution uniforme entre 0 et 2. On retrouve cette distribution sur la courbe ci-dessous qui illustre la répartition de l'énergie du système avant sa thermalisation. Le nombre de particules est N = 500. La température est T = 10.0.

Etat du systeme avant thermalisation

Après la thermalisation, la répartition de l'énergie apparait ainsi :

Etat du systeme après thermalisation

On observe que la majorité des particules se concentre dans les niveaux d'énergie les plus bas. Il serait intéressant de tracer la courbe de distirbution de l'énergie pour vérifier qu'elle répond bien à la distribution de Boltzmann...

Influence de la température du thermostat

Conservons les mêmes paramètres N = 500 particules et la répartition aléatoire de l'énergie sur le réseau entre 0 et 2. Nous allons augmenter la température du thermostat en la portant de 10.0 à 30.0 et en augmentant le nombre d'itérations à 106 itérations.

L'évolution vers l'équilibre est similaire mais un peu plus lente (d'où l'augmentation du nombre d'itérations) :

Evolution vers l'état d'équilibre à 30

Vous noterez que l'énergie totale à l'équilibre est supérieure à celle observée pour T = 10.0, ce qui est conforme aux lois de la thermodynamique.

Voyons maintenant la répartition de l'énergie pour chaque particule :

Etat du systeme après thermalisation

Comme on pouvait s'y attendre, l'énergie moyenne d'une particule a augmenté et les valeurs d'énergie le plus hautes ont aussi augmenté : la distribution de Boltzmann dépend de la température...

Contenu et design par Dominique Lefebvre - www.tangenteX.com février 2019 Licence Creative Commons Contact :