Monte-Carlo

1. Introduction

Ce document constitue une introduction aux méthodes de Monte-Carlo appliquées à la physique statistique. Les méthodes sont exposées à partir d'un exemple, le modèle d'Ising du ferromagnétisme.

2. Modèle d'Ising

2.a. Définition

Les matériaux ferromagnétiques (comme le fer ou le nickel) tiennent leurs propriétés magnétiques du spin des électrons itinérants (responsables aussi de la conduction électrique). Le spin est un moment cinétique intrinsèque de l'électron, auquel est associé un moment magnétique. Le spin est quantifié. La projection du moment cinétique de spin sur un axe Z n'a que deux valeurs possibles :

S_{z} = \pm \frac{ℏ}{2} (1)

où $ℏ$ est la constante de Planck divisée par 2π. Le moment magnétique associé au spin est proportionnel au moment cinétique :

m_{z} = - g μ_{B} \frac{S_{z}}{ℏ} (2)

où g est le facteur de Landé (égal environ à 2) et μ_B le magnéton de Bohr. En présence d'un champ magnétique, l'énergie de l'électron est :

E = - m_{z} B_{z} (3)

Les moments magnétiques ont donc tendance à s'orienter dans le sens du champ.

Dans un matériau ferromagnétique, les spins ont la faculté de s'orienter dans le même sens, même en l'absence de champ magnétique extérieur. C'est pourquoi un morceau de fer peut présenter une aimantation permanente. Pour expliquer ce phénomène, il faut une interaction entre les spins.

Deux électrons peuvent interagir par l'interaction dipolaire : le champ magnétique créé par le moment magnétique d'un électron exerce une influence sur le moment magnétique de l'autre électron. Cependant, l'interaction dipolaire est beaucoup trop faible pour expliquer le ferromagnétisme.

C'est un phénomène quantique qui explique l'interaction entre deux spins responsable du ferromagnétisme. Le principe d'exclusion de Pauli et la répulsion électrique entre les deux électrons conduisent à une interaction appelée interaction d'échange. Pour deux électrons voisins i et j, cette interaction peut être représentée par une énergie de la forme :

E = - S_{z}^{i} S_{z}^{j} (4)

Deux électrons voisins ont ainsi une énergie plus basse lorsque leurs spins sont de même sens. Ils ont donc tendance à s'orienter dans le même sens.

Le modèle d'Ising est une représentation simplifiée du matériau ferromagnétique, dans lequel la direction des spins (axe Z) est la même pour tous les électrons. Les spins sont positionnés sur un réseau. On se limite ici au modèle d'Ising à une dimension. La figure suivante montre un exemple de configuration pour un système à N spins. Les deux valeurs possibles pour chaque spin sont -1 et 1.

Figure pleine page

Le spin d'indice i interagit avec ses deux voisins i-1 et i+1. On introduit une condition limite périodique, consistant à faire interagir les spins 0 et N-1. Avec ces hypothèses, l'énergie du système s'écrit :

E = - \sum_{i = 0}^{N - 1} B S_{i} + S_{i} S_{i + 1} (5)

où B est une constante représentant le champ magnétique extérieur. Par convention, le moment magnétique est dans le même sens que le spin (physiquement, ils sont de sens opposé car l'électron a une charge négative).

2.b. Distribution de Boltzmann

En raison de l'agitation thermique, les électrons itinérants échangent en permanence de l'énergie avec les atomes du métal et les autres électrons. L'agitation thermique a tendance à maintenir un état désordonné des spins. Un résultat important de la théorie statistique est la distribution de Boltzmann, valable si le système est à l'équilibre thermique à la température T. Soit μ une configuration particulière du système de spins et E_μ l'énergie de cette configuration. La propabilité de cette configuration est proportionnelle à un facteur exponentiel (facteur de Boltzmann) :

p_{μ} = C exp (- \frac{E_{μ}}{k T}) (6)

où k est la constante de Boltzmann. Lorsque la température est faible, les configurations de faible énergie sont beaucoup plus probables que les configurations de plus haute énergie. Lorsque la température augmente, l'énergie a de moins en moins d'influence sur les rapports de probabilités. À température infinie, toutes les configurations sont équiprobables. Pour simplifier les notations, on pose :

β = \frac{1}{k T} (7)

Pour déterminer la constante de normalisation C, il faut écrire la condition de normalisation : la somme des probabilités de toutes les configurations doit être égale à 1. On obtient ainsi :

p_{μ} = \frac{e^{- β E_{μ}}}{\sum_{μ} e^{- β E_{μ}}} (8)

Cette distribution permet en principe de calculer la moyenne de différentes grandeurs physiques. Soit X une grandeur physique et X_μ sa valeur pour la configuration μ. La valeur moyenne de X, c'est-à-dire l'espérance de la variable aléatoire X, est :

\bar{X} = E (X) = \frac{\sum_{μ} X_{μ} e^{- β E_{μ}}}{\sum_{μ} e^{- β E_{μ}}} (9)

Dans le cas du modèle d'Ising, une grandeur intéressante est le moment magnétique total, défini par :

m_{μ} = \sum_{i} S_{i} (10)

Le moment magnétique est une grandeur macroscopique mesurable. La valeur mesurée est précisément la moyenne définie ci-dessus.

3. Méthodes de Monte-Carlo

3.a. Principe

Il s'agit d'évaluer l'espérance $(9)$ . Même lorsque le nombre de configurations est fini, cette somme ne peut être calculée directement car ce nombre est beaucoup trop grand. Par exemple pour le modèle d'Ising, le nombre de configurations est 2^N. Même pour un petit système de N=100 spins, le nombre de configurations à explorer dépasse de très loin les capacités d'un ordinateur.

Il faut donc s'orienter vers une méthode de Monte-Carlo, qui consiste à opérer sur un échantillon de configurations choisies aléatoirement.

3.b. Échantillonnage direct

Une première idée serait de choisir des configurations aléatoirement, toutes ayant la même probabilité. Dans le cas du modèle d'Ising, ce tirage ne pose pas de difficultés. Il suffit de choisir aléatoirement chaque spin (-1 ou 1 équiprobables) pour obtenir une configuration. Si M est le nombre de configurations tirées, l'espérance $(9)$ est évaluée par la moyenne empirique :

m o y (X, M) = \frac{1}{M} \sum_{k = 1}^{M} p_{k} X_{k} (11)

où p_k est la probabilité de la configuration k donnée par la loi de Boltzmann.

Cette méthode est cependant impossible à réaliser, car la somme figurant au dénominateur de l'expression de p_μ n'est pas connue a priori. C'est justement ce type de somme que l'on cherche à calculer.

La solution consiste à tirer aléatoirement les échantillons avec la probabilité p_μ et non pas avec une probabilité uniforme. De cette manière, l'évaluation de l'espérance devient :

m o y (X, M) = \frac{1}{M} \sum_{k = 1}^{M} X_{k} (12)

La difficulté est d'avoir une méthode pour échantillonner les configurations avec la probabilité p_μ. Les différentes méthodes d'échantillonnage de nombres aléatoires sont exposées dans Échantillonnage des distributions de probabilité discrètes.

La méthode du rejet ne nécessite pas la connaissance de la constante de normalisation. Elle consiste à tirer aléatoirement une configuration μ avec un tirage uniforme, puis à calculer sa probabilité non normalisée :

q_{μ} = e^{- β E_{μ}} (13)

La seconde étape consiste à tirer un nombre réel aléatoire x avec une densité uniforme sur l'intervalle [0,q_max]. Dans le cas présent, la valeur maximale de la probabilité non normalisée est q_max=1. Ce nombre x est comparé à q_μ, ce qui permet de conserver ou de rejeter la configuration. Si x<q_μ, la configuration est retenue, sinon elle est rejetée.

Pour le modèle d'Ising, la méthode du rejet est donc constituée des étapes suivantes :

Choisir aléatoirement les N spins, les valeurs -1 et +1 étant équiprobables.
Calculer l'énergie E_μ de la configuration et la probabilité q_μ non normalisée correspondante.
Tirer aléatoirement un réel x dans l'intervalle [0,1] avec une densité uniforme.
Si x<q_μ, la configuration est retenue, sinon elle est rejetée.

L'évaluation de la moyenne se fait avec les configurations retenues. On évalue aussi la variance comme expliqué dans méthode de Monte-Carlo. En divisant cette variance par le nombre M de tirages, on obtient la variance de la moyenne.

Dans le cas de la distribution exponentielle de Boltzmann, la méthode du rejet souffre d'un grave inconvénient : la très grande majorité des configurations a une probabilité non normalisée très faible, ce qui conduit à un taux de rejet très important, surtout lorsque la température est faible.

3.c. Méthode de Metropolis

La méthode de Metropolis ([1], [2]) a servi de point de départ pour le développement de méthodes d'échantillonnage par chaines de Markov, qui sont incomparablement plus efficaces que les méthodes d'échantillonnage direct. L'article original exposait son application au modèle des sphères dures, mais la méthode de Metropolis se généralise a tout problème de calcul statistique.

On se contente ici d'un exposé descriptif de la méthode appliquée au modèle d'Ising.

Au lieu de tirer directement chaque configuration, la méthode de Metropolis tire une nouvelle configuration en modifiant légèrement la configuration précédente. Il s'agit de générer une suite de configurations avec des règles de transition, de telle sorte que les configurations obtenues obéissent à la loi de probabilité voulue, en l'occurence la loi de Boltzmann.

L'obtention d'une nouvelle configuration à partir de la précédente se fait en deux étapes : la sélection et l'acceptation. En cas de refus de la nouvelle configuration, la configuration est conservée.

Figure pleine page

La sélection est une modification de la configuration simple à réaliser. Toute règle de sélection est possible, à condition qu'elle permette l'exploration de toutes les configurations. Dans le modèle d'Ising, on se contente de choisir un spin aléatoirement et de changer son signe. Notons μ la configuration et μ' la nouvelle configuration. Le rapport de leur probabilité est :

\frac{p_{μ'}}{p_{μ}} = e^{- β (E_{μ'} - E_{μ})} = e^{- β Δ E} (14)

Lorsqu'on modifie un seul spin, il est très facile (et peu coûteux) de calculer la variation d'énergie

Δ E = E_{μ'} - E_{μ} (15)

Lorsque cette différence est négative, la nouvelle configuration est plus probable que la précédente. Il est donc logique de l'accepter. Si la différence d'énergie est positive, il ne faut pas refuser systématiquement la nouvelle configuration, car cela pourrait empêcher d'explorer certaines configurations importantes. Il faut donc l'accepter ou la refuser en utilisant une méthode similaire à la méthode du rejet exposée plus haut. On tire un nombre réel aléatoire avec une densité de probabilité uniforme sur l'intervalle [0,1]. La nouvelle configuration est acceptée seulement si :

x < e^{- β Δ E} (16)

Dans le cas contraire, la nouvelle configuration est refusée et on garde la précédente.

Voici pour résumer les opérations à effectuer pour obtenir une nouvelle configuration, dans le cas du modèle d'Ising :

Choix d'un spin aléatoirement.
Calcul de la variation d'énergie ΔE résultant de l'inversion de ce spin.
Si ΔE<=0, le spin est inversé.
Sinon, tirage d'un réel x aléatoirement avec une densité uniforme sur l'intervalle [0,1]. Si $x < e^{- β Δ E}$ , le spin est inversé, sinon il est inchangé.

La suite de configurations ainsi générée converge vers la distribution de Boltzmann ([2]). La configuration initiale est quelconque; on peut la choisir aléatoirement. Il faut calculer un certain nombre Q de configurations (qui peut être très grand) avant que la suite de configurations obéisse effectivement à la loi de Boltzmann. À partir de cet instant, le système est à l'équilibre et on peut commencer à calculer des grandeurs moyennes. Si M est le nombre total de configurations d'équilibre générées, l'évaluation de l'espérance de X est donnée par la somme suivante :

m o y (X, M) = \frac{1}{M} \sum_{k = 1}^{M} X_{k} (17)

On calcule aussi la variance expérimentale.

Q est le nombre de configurations qu'il faut générer pour atteindre l'équilibre. Les calculs de moyenne ne peuvent commencer qu'à partir de ce rang. La principale difficulté de la méthode de Metropolis est l'évaluation de Q, qui n'est pas connu a priori.

Comparée à la méthode d'échantillonnage directe (méthode du rejet), la méthode de Metropolis est bien meilleure pour deux raisons :

Le calcul de ΔE est beaucoup moins coûteux que le calcul de l'énergie E d'une configuration.
La proportion de configurations rejetées est beaucoup plus faible que dans la méthode du rejet.

Contrairement à la méthode directe, la méthode de Metropolis nécessite un certain nombre de tirages de mise à l'équilibre (Q), avant que l'on puisse commencer à faire les calculs de moyenne.

Pour améliorer une méthode de Metropolis, il faut trouver des règles de sélection qui donne des taux de rejet faibles dans l'étape d'acceptation. Une méthode de Metropolis idéale aurait un taux de rejet nul.

4. Travaux pratiques : implémentation de la méthode de Metropolis

L'objectif est d'implémenter la méthode de Metropolis pour le modèle d'Ising à une dimension. Les deux variables aléatoires dont on cherche à calculer la moyenne sont le moment magnétique, c'est-à-dire la somme des spins :

m = \sum_{i = 0}^{N - 1} S_{i} (18)

et l'énergie :

E = - \sum_{i = 0}^{N - 1} B S_{i} + S_{i} S_{i + 1} (19)

Compte tenu du très grand nombre de configurations qu'il faudra générer, il faut prendre soin de minimiser les calculs. En particulier, le moment magnétique et l'énergie devront être stockés dans une variable, qui sera actualisée à chaque basculement de spin.

Lorsque l'inversion d'un spin est testé, il faut calculer la variation d'énergie ΔE résultant de l'inversion de ce spin. Notons S l'état du spin avant l'inversion. Le terme d'interaction avec le champ extérieur est changé de +2BS. Pour le terme d'interaction entre spins voisins, il y a trois cas suivant l'état des spins avant le changement :

Les deux spins voisins S_i-1 et S_i+1 sont de même signe que S_i. La variation d'énergie est +4.
Les deux spins voisins sont de même signe mais de signe opposé à S_i. La variation d'énergie est -4.
Les deux spins voisins sont de signes opposés. La variation d'énergie est nulle.

Finalement, il y a 6 variations d'énergie totale possibles :

2 B + 4, 2 B - 4, 2 B, - 2 B + 4, - 2 B - 4, - 2 B (20)

Le facteur e^-βΔE de ces 6 termes devra être calculé au préalable et stocké dans une liste, car le calcul d'une exponentielle est une opération très coûteuse qu'il faut éviter de faire à chaque pas de calcul.

import numpy
import random
import math
from matplotlib.pyplot import *

Pour mémoriser toutes les données relatives à un réseau de spins (tableau de spins, énergie, moment magnétique), il est commode d'utiliser une classe. On définit donc une classe dont voici le prototype :

class Ising1D:
    def __init__(self,N):
    def constantes(self,B,T):
    def init_random(self):
    def voisins(self,i):
    def deltaE(self,i):
    def metropolis(self,i,x):
    def boucle(self,n):

Le constructeur __init__(self,N) se charge de créer un tableau de N spins, initialisés à 1. On rappelle que toute variable interne à la classe doit avoir son nom précédé de self..

La fonction constantes(self,B,T) définit le champ magnétique et la température. Dans cette fonction, on calcule les termes exponentiels correspondant aux différentes variations d'énergie possibles et on les mémorise dans un tableau.

La fonction init_random(self) initialise les spins avec des valeurs aléatoires (-1 ou 1). Cela correspond à une configuration initiale avec une température infinie. La fonction doit aussi initialiser l'énergie et le moment magnétique (self.E et self.M).

La fonction voisins(self,i) renvoie l'état des deux voisins du spin d'indice i, en tenant compte des conditions limites périodiques.

La fonction deltaE(self,i) renvoie la variation d'énergie ΔE et le facteur de Boltzmann e^-βΔE pour le basculement du spin d'indice i.

La fonction metropolis(self,i,x) teste le basculement du spin d'indice i et réalise ou pas ce basculement en fonction de ΔE et de x, qui est un nombre aléatoire flottant compris entre 0 et 1.

La fonction boucle(self,n) réalise une boucle de calcul comportant n*N tentatives (réussies ou pas) de basculement de spin. Cela correspond, en moyenne, à n basculements pour chaque spin du réseau. La fonction renvoie la liste des valeurs de moment magnétique (une pour chaque itération) et la liste des valeurs de l'énergie.

class Ising1D:
    def __init__(self,N):
        self.N = N
        self.spin = numpy.ones(self.N)
        self.E = 0 # à calculer
        self.M = self.spin.sum()

    def constantes(self,B,T):
        beta = 1/T # T : k*T
        self.T = T
        self.B = B
        self.dE = numpy.array([2*B+4,2*B-4,2*B,-2*B+4,-2*B-4,-2*B])
        self.exp = numpy.exp(-beta*self.dE)

    def init_random(self):
        for i in range(self.N):
            self.spin[i] = random.randint(0,1)*2-1
        self.E = 0.0
        self.M = 0.0
        for i in range(self.N-1):
            self.E += -self.B*self.spin[i]+self.spin[i]*self.spin[i+1]
            self.M += self.spin[i]
        i = self.N-1
        self.E += -self.B*self.spin[i]+self.spin[i]*self.spin[0]
        self.M += self.spin[i]

    def voisins(self,i):
        # s1,s2,s3
        if i==0:
            s1 = self.spin[self.N-1]
            s3 = self.spin[i+1]
        elif i==self.N-1:
            s1 = self.spin[i-1]
            s3 = self.spin[0]
        else:
            s1 = self.spin[i-1]
            s3 = self.spin[i+1]
        return (s1,s3)

    def deltaE(self,i):
        s2 = self.spin[i]
        (s1,s3) = self.voisins(i)
        indice = 0
        if s2==-1:
            indice = 3
        if (s1,s2,s3)==(1,1,1) or (s1,s2,s3)==(-1,-1,-1):
            pass
        elif (s1,s2,s3)==(1,-1,1) or (s1,s2,s3)==(-1,1,-1):
            indice += 1
        else:
            indice += 2
        return (self.dE[indice],self.exp[indice])

    def metropolis(self,i,x):
        (dE,exp) = self.deltaE(i)
        if dE<=0:
            self.spin[i] = - self.spin[i]
            self.E += dE
            self.M += 2*self.spin[i]
        else:
            if x<exp:
                self.spin[i] = - self.spin[i]
                self.E += dE
                self.M += 2*self.spin[i]

    def boucle(self,n):
        M = numpy.zeros(n)
        E = numpy.zeros(n)
        for k in range(n):
            for l in range(N):
                i = random.randint(0,self.N-1)
                x = random.random()
                self.metropolis(i,x)
            M[k] = self.M
            E[k] = self.E
        return (M,E)

Voici un exemple :

N = 1000
ising = Ising1D(N)
B = 1.0
T = 0.1
ising.constantes(B,T)
ising.init_random()
(M,E) = ising.boucle(100)

figure()
subplot(211)
plot(M/N)
ylim(-2,2)
ylabel("M/N")
subplot(212)
plot(E/N)
ylim(-2,2)
ylabel("E/N")

figA.pdf

On voit sur ces courbes le moment magnétique (par spin) et l'énergie en fonction du nombre d'itérations (une itération = N tentatives de basculement). Ce type de tracé donne une estimation du nombre d'itérations qu'il faut pour atteindre l'équilibre. Pour confirmer, on refait une simulation avec une condition initiale aléatoire différente :

ising.init_random()
(M,E) = ising.boucle(100)

figure()
subplot(211)
plot(M/N)
ylim(-2,2)
ylabel("M/N")
subplot(212)
plot(E/N)
ylim(-2,2)
ylabel("E/N")

figB.pdf

Cela confirme que l'équilibre est largement atteint au bout de 100 itérations (mais on ne peut en être certain). On peut donc commencer les calculs de moyennes et d'écart-types. Voici par exemple un calcul avec 100 itérations :

n = 100
(M,E) = ising.boucle(n)
M /= N
E /= N
M_moy = M.mean()
M_std = M.std()
E_moy = E.mean()
E_std = E.std()

print((M_moy,M_std))
--> (1.0, 0.0)

print((E_moy,E_std))
--> (-2.0960000000000001, 0.0)

L'écart-type de l'estimation de la moyenne est obtenu en divisant par la racine carré du nombre de configurations utilisées pour estimer la moyenne :

print(M_std/math.sqrt(n))
--> 0.0

print(E_std/math.sqrt(n))
--> 0.0

Lorsque le programme fonctionnera correctement, on testera différentes valeurs de la température pour ce champ magnétique, puis différents champs magnétiques. Pour B=0, on constatera que le modèle d'Ising à une dimension ne peut donner de moment magnétique permanent. Pour obtenir une aimantation en l'absence de champ B (ferromagnétisme), il faut un modèle à deux dimensions.

Méthodes de Monte-Carlo en physique statistique