Méthodes de Monte-Carlo.

Michel ROGER
Service de Physique de l’Etat Condensé
CEA Saclay

13 octobre 2008
2
 Préface

Les chapitres 1 à 8 de ce cours concernent les bases essentielles de la méthode de Monte-

Carlo. Elles sont enseignées dans le cadre commun de trois Masters de l’INSTN :
– “Physique Médicale”
– “Rayonnement & Energie”
– “Physique & Environnement”
Cet enseignement est précédé de six séances d’une heure et demie destinées à des rap-
pels fondamentaux de théorie des probabilités (Annexes A, B) et de Physique Statistique
(Annexes C, D).
Les chapitres 9 à 12 ne sont enseignés qu’aux étudiants des masters “Rayonnement &
Energie” et “Physique & Environnement”. Ils abordent des aspects plus spécifiques de le
méthode de Monte-Carlo, en particulier les problèmes d’optimisation, la simulation des
polymères et le traitement de problèmes quantiques.
4
 Table des matières

1 Introduction 11

2 Calcul d’intégrales 13
2.1 Travail dirigé introductif : l’aiguille de Buffon . . . . . . . . . . . . . . . . . 13
2.1.1 Expérience . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1.2 Montrer que M/N “tend vers” 2L/(πD) . . . . . . . . . . . . . . . . 14
2.2 Calcul d’intégrales multidimensionnelles. . . . . . . . . . . . . . . . . . . . 17

3 Techniques d’échantillonnage 19
3.1 Génération de nombres “pseudo-aléatoires” . . . . . . . . . . . . . . . . . . 19
3.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.1.2 La méthode des congruences linéaires . . . . . . . . . . . . . . . . . 20
3.1.3 Le “test spectral” . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.1.4 Suites de Fibonacci “retardées” . . . . . . . . . . . . . . . . . . . . . 26
3.1.5 Puis-je utiliser la fonction “RAND(), random(), ...” de mon ordinateur ? 27
3.2 Echantillonnage d’une loi de probabilité non uniforme . . . . . . . . . . . . 28
3.2.1 Méthode de transformation . . . . . . . . . . . . . . . . . . . . . . . 28
3.2.2 Méthode de réjection de Von Neumann . . . . . . . . . . . . . . . . 32

4 Echantillonnage suivant l’importance. 35

4.1 Calcul stochastique d’une intégrale . . . . . . . . . . . . . . . . . . . . . . . 35
4.2 Réduction de la variance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

5 Simulation de systèmes statistiques 41

6

5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.2 Chaı̂nes de Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5.2.1 Définition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5.3 Convergence vers une loi de distribution invariante imposée . . . . . . . . . 43
5.3.1 Echantillonnage de f (x) = e−βE(x) /Z . . . . . . . . . . . . . . . . . . 44
5.4 Algorithme de Metropolis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

6 Le modèle d’Ising 47
6.1 L’Hamiltonien d’Ising . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
6.2 Modélisation par la méthode de Monte-Carlo . . . . . . . . . . . . . . . . . 47
6.3 Analyse des résultats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
6.4 Evolution en fonction de la température . . . . . . . . . . . . . . . . . . . . 52

7 Utilisation de biais. Algorithme de Wolff 55

7.1 Introduction de biais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
7.2 Exemple : algorithme de Wolff pour le modèle d’Ising . . . . . . . . . . . . 56

8 Applications à la physique médicale 59

8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
8.2 Modélisation d’un rayonnement . . . . . . . . . . . . . . . . . . . . . . . . . 59
8.3 Modélisation de la source . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
8.4 Modélisation d’une trajectoire . . . . . . . . . . . . . . . . . . . . . . . . . . 60
8.5 Interaction photon-matière . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
8.5.1 Diffusion élastique de Rayleigh . . . . . . . . . . . . . . . . . . . . . 62
8.5.2 Diffusion inélastique de Compton . . . . . . . . . . . . . . . . . . . . 63
8.5.3 Absorption photoélectrique . . . . . . . . . . . . . . . . . . . . . . . 64
8.5.4 Choix d’un des événements . . . . . . . . . . . . . . . . . . . . . . . 64
8.5.5 Pour en savoir plus... . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
8.6 Propagation dans un milieu inhomogène . . . . . . . . . . . . . . . . . . . . 64

9 Méthode de Monte Carlo et Dynamique Moléculaire 67

9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
9.2 Equation stochastique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
9.3 Comparaison avec la Dynamique Moléculaire . . . . . . . . . . . . . . . . . 69
 7

10 Optimisation 73
10.1 Méthode du “Recuit Simulé” (Simulated Annealing) . . . . . . . . . . . . . 73
10.2 “Recuit Parallèle” (Parallel Tempering) . . . . . . . . . . . . . . . . . . . . 73

11 Polymères 77
11.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
11.2 Modélisation de polymères . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
11.2.1 Modèle du “collier de perles” . . . . . . . . . . . . . . . . . . . . . . 77
11.2.2 Modèle du “Chapelet de ressorts” (String beads) . . . . . . . . . . . 78
11.3 Mouvements dans un algorithme de Metropolis . . . . . . . . . . . . . . . . 79
11.3.1 Reptation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
11.3.2 Mouvements de pivot . . . . . . . . . . . . . . . . . . . . . . . . . . 79
11.4 “Reconstruction biaisée” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

12 Monte-Carlo quantique 83
12.1 Monte-Carlo Variationnel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
12.2 “Diffusion Monte-Carlo” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
12.3 “Path Integral Monte-Carlo” . . . . . . . . . . . . . . . . . . . . . . . . . . 88
12.3.1 Valeur moyenne d’une grandeur physique . . . . . . . . . . . . . . . 88
12.3.2 La Matrice Densité . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
12.3.3 Echantillonnage de ρ(R, R 0 ; β) . . . . . . . . . . . . . . . . . . . . . . 90

A Rappels de théorie des probabilités 95

A.1 Expérience, ensemble des résultats possibles. . . . . . . . . . . . . . . . . . . 95
A.2 Notion d’événement aléatoire. . . . . . . . . . . . . . . . . . . . . . . . . . . 95
A.3 Notion de probabilité. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
A.3.1 Définition heuristique . . . . . . . . . . . . . . . . . . . . . . . . . . 96
A.3.2 Définition mathématique . . . . . . . . . . . . . . . . . . . . . . . . 97
A.4 Probabilités conditionnelles. Evénements indépendants . . . . . . . . . . . . 97
A.5 Variable aléatoire. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
A.6 Loi de probabilité d’une variable aléatoire . . . . . . . . . . . . . . . . . . . 99
A.7 Fonction de répartition d’une variable aléatoire réelle . . . . . . . . . . . . . 99
A.8 Densité de probabilité. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
8

A.9 Espérance, moments d’ordre n, variance . . . . . . . . . . . . . . . . . . . . 101

A.9.1 Valeur moyenne ou “Espérance” d’une variable aléatoire réelle . . . 101
A.9.2 Moment d’ordre n. Variance d’une variable aléatoire réelle . . . . . . 102
A.10 Fonction Caracteristique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
A.10.1 Définition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
A.10.2 Exemple : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
A.11 Inégalité de Bienaymé-Tchebichef . . . . . . . . . . . . . . . . . . . . . . . . 104
A.12 Couples de variables aléatoires . . . . . . . . . . . . . . . . . . . . . . . . . 105
A.12.1 Covariance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
A.12.2 Indépendance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
A.12.3 Somme de 2 variables aléatoires indépendantes . . . . . . . . . . . . 105
A.13 Généralisation à N variables aléatoires . . . . . . . . . . . . . . . . . . . . . 107
A.14 Somme de variables aléatoires indépendantes . . . . . . . . . . . . . . . . . 107
A.14.1 Lois “faible” des grands nombres. . . . . . . . . . . . . . . . . . . . . 107
A.14.2 Théorème de la limite centrale . . . . . . . . . . . . . . . . . . . . . 109
A.15 Divers modes de convergence . . . . . . . . . . . . . . . . . . . . . . . . . . 111

B Chaı̂nes de Markov 113

B.1 Définition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
B.2 Chaı̂nes dans un ensemble discret . . . . . . . . . . . . . . . . . . . . . . . . 114
B.2.1 Propriétés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
B.2.2 Rayon spectral de la matrice p . . . . . . . . . . . . . . . . . . . . . 115
B.2.3 Convergence vers une loi de distribution stationaire . . . . . . . . . 116
B.2.4 Chaı̂ne de Markov ergodique . . . . . . . . . . . . . . . . . . . . . . 116
B.2.5 Généralisation à un ensemble d’états continus . . . . . . . . . . . . . 117

C Rappels de Physique statistique 119

C.1 Entropie statistique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
C.1.1 Définition de Shannon . . . . . . . . . . . . . . . . . . . . . . . . . . 119
C.1.2 Théorème de Khinchin (1957). . . . . . . . . . . . . . . . . . . . . . 120
C.2 Systèmes macroscopiques à l’équilibre . . . . . . . . . . . . . . . . . . . . . 120
C.2.1 Ensemble micro-Canonique . . . . . . . . . . . . . . . . . . . . . . . 120
C.2.2 Ensemble Canonique . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
 9

C.2.3 Ensemble Grand Canonique . . . . . . . . . . . . . . . . . . . . . . . 123

C.3 Evolution vers l’équilibre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
C.3.1 Probabilité de transition . . . . . . . . . . . . . . . . . . . . . . . . . 124
C.3.2 Equation maı̂tresse . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
C.3.3 Système isolé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
C.3.4 Système S en contact avec un thermostat T . . . . . . . . . . . . . . 125
C.3.5 Propriétés à l’équilibre. Relation du bilan détaillé. . . . . . . . . . . 126

D Valeurs moyennes d’observables 129

D.1 Définition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
10
 1 Introduction

La méthode de Monte-Carlo a vu son essor à partir de la fin de la seconde guerre mondiale,

essentiellement dans le cadre du projet americain “Manhattan” concernant le développement
de l’arme nucléaire. Cette époque correspond également à la construction des premiers “or-
dinateurs”. Ce projet étant été classé “secret défense”, il est difficile de savoir exactement
qui parmi ses pionniers : Von Neumann, Ulam, Metropolis a proposé le nom de “Monte-
Carlo”. Quoi qu’il en soit, ce terme fait référence aux jeux de hasard : la capitale de la
principauté de Monaco, avec ses casinos, ayant dès cette époque une solide réputation dans
la société mondaine internationale.
En l’occurrence, la méthode de Monte-Carlo consiste à simuler sur ordinateur un phénomène
stochastique par essence : la diffusion de neutrons dans des matériaux fissiles.
Un ordinateur effectue une suite d’opérations parfaitement déterministes, comment alors simuler le
hasard sur un ordinateur ? On génère des suites de nombres “pseudo aléatoires” qui ont les mêmes
propriétés statistiques que des nombres aléatoires, mais qui ne sont pas aléatoires. Les algorithmes
les plus souvent utilisés seront étudiés au Chapitre 3.
Les applications en physique médicale que nous verrons au Chapitre 8, sont du même ordre.
Par exemple, à partir d’une source de rayon X ou γ dont on connait les caractéristiques
(distribution des photons en energie, en direction) et des processus physiques bien maı̂trisés
d’interaction des photons avec la matière (tissus, milieu osseu, ...), on va par exemple simuler
un ensemble de photons arrivant sur une tumeur cancéreuse, leur diffusion sur les cellules,
et on va en déduire un certains nombres de grandeurs statistiques, telles que le nombre
moyen de photons absorbés à une profondeur x etc...
D’autre part, nous verrons des applications de la méthode de Monte Carlo à la physique
statistique (Chapitre 5). En mécanique statistique classique, la loi de Boltzman nous dit
que pour un système à grand nombre de degrés de libertés, en contact avec un thermostat
qui le maintient à la température T , la probabilité pour que le système soit dans un état de
configuration X (X représente l’ensemble de ses degrés de liberté : positions de particules
etc ...) est proportionnelle à exp[−E(X)/kT ] où E(X) est l’energie correspondant à la
configuration X. La méthode de Monte-Carlo va nous permettre d’engendrer un grand
nombre de configurations, suivant une loi de probabilité proportionnelle à exp[−E(X)/kT ],
et de faire des moyennes statistiques sur les grandeurs macroscopiques caractéristiques du
12 Introduction

système : énergie moyenne, chaleur spécifique, etc...

Mais nous allons commencer par un autre domaine très important de la méthode de Monte-
Carlo : l’utilisation de méthodes stochastiques pour une réponse approchée à un problème
mathématique exact : par exemple l’estimation d’une intégrale multidimensionnelle (qui, à
priori, n’a rien d’aléatoire !).
Dans ce domaine, il est d’usage de citer comme premier exemple historique de la méthode
(qui ne s’appellait pas encore “de Monte-Carlo” !) le jeu mathématique proposé par le
Comte de Buffon (plus célèbre pour ses travaux naturalistes) en 1777. Ce jeu consiste à
jeter N fois au hasard une aiguille sur un parquet formé de lattes parallèles et identiques
et à compter le nombre M de fois où l’aiguille intercepte la ligne de jointure entre deux
lattes. Si la longueur L de l’aiguille est inférieure à la largeur D d’une latte, alors le rapport
M/N “tend” vers 2L/(πD). Et Laplace a fait remarquer un peu plus tard (1820) que ce
jeu constituait une méthode expérimentale de détermination du nombre π.
 2 Calcul d’intégrales

2.1 Travail dirigé introductif : l’aiguille de Buffon

A travers cet exemple simple, nous allons revoir l’utilisation des notions et principaux
résultats de la théorie des probabilités : loi des grands nombres, théorème de la limite
centrale, convergences, estimation d’erreurs etc... fondamentales pour ce cours ! Un rappel
de toutes ces notions est effectué en annexe A.

2.1.1 Expérience

Chacun dispose d’une feuille de papier sur laquelle sont tracées des lignes équidistantes de
D et d’une aiguille de longueur L. Lancer l’aiguille N = 100 fois et compter le nombre M
de fois où celle-ci intercepte une ligne. Disposer dans un tableau les valeurs M i /N obtenues
par chacun des n étudiants. Calculer une estimation le la valeur moyenne globale :
i=n
1 X Mi
m=
n N
i=1
2 pour N=100 lancers
ainsi qu’une estimation de la variance σ N
i=n  2
2 1 X Mi
σN = −m
n−1 N
i=1

En déduire une valeur approchée du nombre π et en estimer la “précision” obtenue sur les
n répétitions des N lancers. .

En langage probabiliste :
– Lancer l’aiguille est une expérience aléatoire, elle est caractérisée par l’ensemble des
résultats possibles Ω constitué par toutes les positions possibles ω de l’aiguille sur la
feuille.
– Le fait de couper une ligne définit un événement aléatoire, représenté par le sous ensemble
A de résultats de Ω qui le réalisent.
14 Calcul d’intégrales

– Si on répête l’expérience un très grand nombre de fois, on observera que le rapport

M/N tend vers une limite, un nombre compris entre 0 et 1, que l’on definira comme la
probabilité de l’événement A.

2.1.2 Montrer que M/N “tend vers” 2L/(πD)

Pour chaque résultat ω de l’expérience on mesure l’angle 0 ≤ Θ(ω) < π entre la direction
de l’aiguille et celle des droites parallèles, et la distance 0 ≤ Y (ω) < D/2 du centre de
l’aiguille à la droite la plus proche. Θ et Y sont deux fonctions de Ω dans R, dont la valeur
dépend du résultat ω d’une expérience aléatoire. Ce sont donc deux variables aléatoires à
valeurs dans R.
On peut considérer que la loi de probabilité de Y est uniforme sur le segment [0, D/2]
(toute valeur de Y entre 0 et D/2 est équiprobable), c’est à dire :

2dy
dPY = P rob({y ≤ Y (ω) < y + dy}) =
D
et que la loi de distribution de Θ est elle-aussi uniforme sur [0, π] :

dθ
dPΘ = P rob({θ ≤ Θ(ω) < θ + dθ}) =
π

L’événement “couper une ligne” se produit si et seulement si

L
Y < sin Θ
2

Le décompte que nous opérons peut lui-même être représenté par une variable aléatoire
Z(ω) fonction des deux précédentes et à valeurs dans l’ensemble {0,1} :

L
Z(ω) = 1 si Y (ω) < sin Θ(ω) et 0 sinon
2

L’espérance E(Z) de la variable aléatoire Z(Y, Θ) est l’intégrale de Z(Y, Θ) par rapport
aux “mesures de probabilités” dPY , dPθ dans le pavé D à deux dimensions représenté sur
la figure 2.1 c’est à dire :
Z Z
I = E(Z) = Z(X, Θ)dPY dPΘ
D

soit :
dθ 2dy
Z Z
I = E(Z) = 1[y< L sin θ]
D 2 π D
où
1[y< L sin θ]
2
2.1 Travail dirigé introductif : l’aiguille de Buffon 15

D/2

y < (L/2) sin θ

0
0 θ π

Fig. 2.1 – L’angle θ est réparti uniformément sur le segment [0, π], la distance y à la ligne
la plus proche est répartie uniformément sur le segment [0, D/2]. La proportion de points
qui tombent sous la courbe mesure l’aire de la surface correspondante.

est la “Fonction Indicatrice” du sous ensemble de points du pavé D = [0, π][0, D] tels que
y < L2 sin θ.
On en déduit : π
dθ 2L
Z
I = E(Z) = L sin θ =
0 πD πD
Si on répête N fois l’experience et qu’on effectue la somme S N des N résultats (0 ou 1)
obtenus pour Z, alors SN définit une nouvelle variable aléatoire, somme de N variables
aléatoires indépendantes Z1 , Z2 , ..., ZN ayant la même loi de probabilité que Z :
SN = Z 1 + Z 2 + · · · + Z N
Alors la “loi forte des grands nombres” nous assure que S N /N converge “presque sûrement”
2L
vers l’espérance E(Z). C’est dans ce sens que le rapport M/N “tend” vers πD . Cette loi
ne nous dit rien sur la vitesse de convergence.
A titre anecdotique : Mario Lanzarini annonce en 1901 avoir calculé π = 3.1415929 en
lançant 3408 fois une aiguille de 2,5 cm sur un feuille de papier comportant des lignes
parallèles équidistantes de 3 cm. VRAI ou FAUX ?
Pour répondre à cette question, nous allons déterminer la “vitesse de convergence”, c’est à
dire estimer “l’erreur” :
SN
N = − E[Z]
N
16 Calcul d’intégrales

lorsque N → ∞. Cette “erreur” N définit une nouvelle variable aléatoire

Or le “Theorème de la Limite Centrale” qui nous dit que la variable aléatoire :
√
N
N
σ
“converge en loi” vers la loi normale (loi Gaussienne réduite), de densité de probabilité :
1 2
f (x) = √ e−x /2
2π
σ est “l’écart type” correspondant à la variable aléatoire Z,
q
σ = E[{Z − E[Z]}2 ] = E[Z 2 ] − (E[Z])2
p

E[Z 2 ] représente le moment d’ordre deux de la variable aléatoire Z, et E[(Z − E[Z]) 2 ] = σ 2

sa “variance” .
Calculons exactement le moment d’ordre deux de la variable aléatoire Z(Θ, Y ).
Z Z
2
E[Z ] = Z(X, Θ)2 dPY dPΘ
D

soit : Z Z n o2 dθ 2dy
2
E[Z ] = 1[y< L sin θ]
D 2 π D
qui peut s’écrire
π
dθ 2L
Z
2
E[Z ] = L sin θ =
0 πD πD
On a donc : s  
2L 2L
σ= 1−
πD πD
En appliquant le théorème
√ de la limite centrale, nous déduisons, par exemple, que la prob-
N
abilité pour que σ |N | soit inférieur à 2 est égale à la probabilité pour qu’une variable
aléatoire gaussienne réduite X de loi
1 2
f (x) = √ e−x /2
2π
soit de module inférieur à 2, soit
2
1
Z
2
P rob(|X| < 2) = √ e−x /2 dx = 0.95
−2 2π
On définit ainsi un “intervalle de confiance à 95%” pour I :
 
SN 2σ SN 2σ
P rob −√ <I< +√ = 0.95
N N N N
2.2 Calcul d’intégrales multidimensionnelles. 17

Revenons donc à l’expérience de Lanzarini. Le rapport L/D est 2.5/3, on en déduit

σ = 0.499
√
avec N = 3408, il obtient avec une probabilité de 95% une erreur inférieure à 2∗0.499/ 3408 =
0.017 sur le nombre 2L/(πD) = 0.5305 soit une erreur relative inférieure à 3%. Il ne peut
donc estimer π qu’avec un seul chiffre significatif après la virgule... c’est un imposteur !
Pour l’expérience réalisée en cours, comparer votre ecart type experimentalement estimé à la
valeur théorique exacte. Avec le nombre total n.M de lancés effectués déduire “l’incertitude”
sur votre mesure de π.

2.2 Calcul d’intégrales multidimensionnelles.

Les méthodes usuelles de quadrature (méthode de Simpson, méthodes de Gauss...) sont

efficaces pour le calcul d’intégrales à une dimension et se généralisent facilement à plusieurs
dimensions lorsque le domaine d’intégration est simple (hypercube, par exemple).

Nous allons voir, sur un exemple : la méthode de Simpson, qu’il est par contre illusoire de
généraliser ces méthodes à un grand nombre de dimensions.
La méthode de Simpson à une dimension revient à couper l’intervalle d’intégration (que
nous supposerons [0,1]) en N segments infinitésimaux de longueur h = 1/N . Nous con-
sidérons un petit segment particulier, et par une translation nous recentrons l’origine au
milieu de cet intervalle. Dans ce nouveau système d’axes, nous approximons la fonction
Φ(x) par son développement de Taylor d’ordre deux :

∂Φ ∂ 2 Φ x2
Φ(x) = Φ(0) + x+ + O(x3 )
∂x ∂x2 2

et nous intégrons cette expression entre −h/2 et h/2. les trois premiers termes sont pris en
h/2
compte exactement. Le premier terme d’erreur d’ordre x 3 donne un terme en [x4 /4]−h/2 qui
s’annulle par symétrie. L’erreur provient donc de l’intégration du terme suivant en x 4 qui
donne une erreur d’ordre h5 . Le domaine d’intégration est divisé en N intervalles égaux,
de longueur h = 1/N l’erreur totale (N fois la précédente) est d’ordre N h 5 ≈ 1/N 4 .
A deux dimensions, on divise le pavé d’intégration [0,1]x[0,1] en N pavés élémentaires de
coté h = (1/N )1/2 . Dans un carré élémentaire, on recentre les axes sur le centre de ce carré
infinitésimal [−h/2, h/2][−h/2, h/2] et on approxime la fonction Φ(x, y) par un polynome
d’ordre 2 en x, y. L’erreur provient des termes d’ordre 3 et 4 négligés. Comme précédement
le terme d’ordre 3 donne une contribution nulle par symétrie. L’erreur provenant des termes
d’ordre 4 s’écrit :
4
1 X ∂ 4 Φ(x, y) h/2
Z Z h/2
4−α
dyy dxxα ≈ O(h6 )
4! α=0 ∂y 4−α ∂xα −h/2 −h/2
18 Calcul d’intégrales

Et à d dimensions, l’erreur dans une cellule élémentaire est d’ordre h 4+d . Si le nombre
de cellules (c’est à dire le nombre de points où on doit calculer la fonction) est N , on a
h = N 11/d et, par rapport à N , l’erreur totale est
 
1
N h4+d ≈ O
N 4/d
Lorsque d est grand, la convergence, en fonction du nombre N de points où la fonction est
évaluée (ce qui détermine le temps CPU) devient
√ très lente. Pour cette règle de Simpson,
la convergence devient plus lente que 1/ N lorsque d est supérieur à 8. Nous voyons
apparaı̂tre√ l’intérêt d’une méthode stochastique, où, comme nous allons le montrer, l’erreur
est en 1/ N quelle que soit la dimension d.

Considérons une intégrale à d dimensions :

Z
I = dxΦ(x) (2.1)

où x = (α1 , α2 , · · · , αd ) est un vecteur à d dimensions et

Z Z Z Z
dx = · · · dα1 dα2 · · · dαd

Nous pouvons écrire Φ(x) sous la forme :

Φ(x) = g(x)f (x) (2.2)

où la fonction f (x) est une densité de probabilité, c’est à dire :

Z
f (x) > 0 ∀x et dxf (x) = 1

Il y a évidemment une infinité de possibilités (que nous exploiterons) pour une telle
décomposition. Nous pouvons donc réécrire :
Z
I = E[g(X)] = g(x)dP avec dP = f (x)dx (2.3)

Considérons un ensemble de N points {x 1 , x2 , · · · xN } répartis suivant la densité de proba-

bilité f (x). Alors, d’après la “loi des grands nombres”, la quantité :
N
SN 1 X
= g(xi ) (2.4)
N N
i=1

tend vers I lorsque N → ∞,√ et le théorème de la limite centrale nous assure que l’erreur
(SN /N − I) est d’ordre 1/ N
 Techniques
3 d’échantillonnage

3.1 Génération de nombres “pseudo-aléatoires”

“Anyone who consider arithmetical methods of producing random digits is, of course, in a
state of sin”.
John von Neumann (1951)

3.1.1 Introduction

Il convient de faire la distinction entre “nombres aléatoires” et “nombres pseudo-aléatoires”.

– “nombres aléatoires” :
Ils peuvent être engendrés par traitement d’une source d’entropie extérieure. Par ex-
emple, un élément radioactif constitue une excellente source d’entropie. L’intervalle de
temps entre deux désintégrations est aléatoire. Autre exemple : un site internet : (
www.random.org ) fournit gratuitement des fichiers préenregistrés (10 Moctets) ou per-
sonnalisés de nombres aléatoires. Ces nombres aléatoires sont obtenus en réglant un
récepteur radio sur une fréquence où personne ne diffuse.
– “nombres pseudo-aléatoires” :
Ce sont des nombres obtenus par un algorithme mathématique et qui ont l’apparence de
nombres aléatoires mais dont la suite est parfaitement reproductible.
A priori, le jugement moral de John von Neumann pourrait nous donner mauvaise con-
science. Sachant toutefois qu’il a lui-même “beaucoup péché”, nous le tempérerons par
cette autre citation de Donald Knuth dont l’un des ouvrages[1] reste la référence en la
matière :

“Random number generators should not be chosen at random”

Nous illustrerons ce conseil sur les deux algorithmes mathématiques les plus utilisés actuelle-
ment :
– La méthode des congruences linéaires
20 Techniques d’échantillonnage

– Les suites de Fibonacci “retardées”

3.1.2 La méthode des congruences linéaires

a) Définition

Cet algorithme proposé par Lehmer[2] génère une suite de nombres entiers compris entre
0 et (m − 1) par la relation :

Xn+1 = (aXn + c) [mod m] (3.1)

à partir d’un nombre initial X0 appelé “semence”.
On peut toujours se ramener à des nombres réels U n appartenant à [0, 1[ en considérant la
suite Un = Xn /m
Il est évident que cette suite est périodique, de période inférieure ou égale à m. Deux
exemples avec m = 16 et a = 5 sont illustrés sur la Figure 3.1 : le premier avec c = 1,
le second avec c = 0. Les nombres successivement obtenus sont représentés sur un cercle.
Quel que soit le point de départ du cycle (semence X 0 ), on y revient après avoir parcouru
le cercle dans le sens trigonométrique.
– Exemple avec c=1 :
On obtient la periode maximale m = 16 et on a engendré une répartition uniforme dans
le sous ensemble des entiers appartenant à [0,15], ce qui semble satisfaisant. Par contre,

(a) (b)

5
13
2 12
11
15
8 6

9 1 9 1

14 0

3
7
4 10
5 13

Fig. 3.1 – (a) Congruence linéaire X n+1 = 5Xn + 1 [mod 16]. (b) Congruence linéaire
Xn+1 = 5Xn [mod 16].
3.1 Génération de nombres “pseudo-aléatoires” 21

on observe une alternance parfaite de nombres pairs et impairs. Dans une représentation
binaire des nombres obtenus, le dernier digit alterne entre 0 et 1 et n’a aucun caractère
aléatoire.
Nous retiendrons de cet exemple que pour engendrer une suite aléatoire de nombres 0
ou 1, il est fortement déconseillé d’utiliser un générateur de nombres entiers pour ne
conserver que le dernier digit. De même, dans le but de gagner du temps CPU, pour
obtenir, par exemple, deux nombres aléatoires de 16 digits on peut imaginer d’engendrer
un nombre de 32 digits puis former 2 nombres avec respectivement les 16 digits les plus
bas et les 16 les plus hauts. Ceci est fortement déconseillé : il faut engendrer deux nombres
de 32 digits dont on ne conservera que les 16 digits les plus hauts (généralement les plus
“aléatoires”)
– Exemple avec c=0 :
Par rapport à l’exemple précédent, on a juste changé la constante additive c. La suite
obtenue a une très courte période 4, et elle est très fortement corrélée puisque la différence
entre deux nombres successifs est 4.
Après cette illustration triviale de l’avertissement de Knuth, cité plus haut, nous allons
donner quelques règles simples pour le choix des paramètres m, a et c. La justification
de ces règles dépasse le cadre de cette introduction, nous renvoyons le lecteur au livre de
Knuth[1] pour la démonstration des théorèmes énoncés.

b) Choix du module m
– Le module m constituant une borne supérieure à la période, on a intérêt à le prendre le
plus grand possible, par exemple on peut choisir le plus grand entier représenté sur la
machine. Pour une machine 32 bits, en langage Fortran, on peut prendre m = 2 31 (le
32ieme ) digit étant réservé au signe. En langage C, on pourra prendre 2 32 en travaillant
avec des entiers positifs (“unsigned int”).
– Avantages et inconvenients à prendre pour m une puissance de 2 :
Si m = 2α , l’opération “modulo m” sur ordinateur devient très rapide, elle correspond
à retenir les α derniers digits du nombre considéré, on la réalise avec un masque et une
opération binaire “AND”. On évite une coûteuse division. De même, il n’y a pas à se
soucier de dépassement éventuel de l’entier maximal représenté lors de la multiplication
de a par Xn , dans ce cas la machine conserve les 32 derniers digits du résultat de la
multiplication et jette le reste.

Par contre, si m = 2l est une puissance de 2, on montre[1] que les l derniers digits de
Xn ont une période inférieure ou égale à 2 l . (Dans l’exemple de la Figure 1a, les deux
derniers digits parcourent le cycle 01, 10, 11, 00 )
– Knuth montre, qu’on peut encore programmer de manière efficace l’opération “modulo
m” lorsque m = 2α − 1 ou m = 2α + 1. Par contre, dans ce cas, les digits les plus bas
sont “aussi aléatoires” que les plus hauts. Ce choix est donc préférable.
c) Choix du facteur multiplicatif a

Le module m étant fixé, il convient de choisir a de manière à obtenir la plus grande période
22 Techniques d’échantillonnage

(inférieure ou égale à m). Ceci n’est pas suffisant pour avoir un bon générateur de nombres
quasi-aléatoires. Il faudra encore passer avec succès un certain nombres de tests statistiques
dont les pricipaux ont été recommandés par Knuth.

Nous énonçons deux théorèmes dont les démonstrations figurent dans le livre de Knuth[1].

Théorème A
La période de la congruence linéaire définie par {m, a, c, X 0 } est m, si et seulement si les
trois conditions suivantes sont simultanément vérifiées :

i) c est different de zéro et est premier par rapport à m

ii) (a − 1) est multiple de p pour tout nombre premier qui divise m
iii) (a − 1) est multiple de 4, si m est multiple de 4

Le cas où c = 0 est plus complexe

Théorème B
Lorsque c = 0, la période est toujours strictement inférieure à m.
Si
m = 2α pβ1 1 · · · pβr r (3.2)
où les p sont des nombres premiers, distincts, impairs, la période maximale possible est le
plus petit multiple commun de :

λ(2α ), λ(pβ1 1 ), · · · , λ(pβr r ) (3.3)

où

λ(pβ ) = pβ−1 (p − 1) si p est impair

λ(1) = λ(2) = 1 ; λ(4) = 2
λ(2α ) = 2α−2 (α > 2)

On obtient la période maximale si les trois conditions suivantes sont simultanément vérifiées :
β β
i) an 6= 1[mod pj j ] pour 0 < n < λ(pj j )
ii) a = 1 [mod 2] si α = 1
a = 3 [mod 4] si α = 2
a = 3 ou 5 [mod 8] si α = 1
iii) X0 et m sont premiers entre eux
3.1 Génération de nombres “pseudo-aléatoires” 23

3.1.3 Le “test spectral”

Parmi les nombreux tests statistiques qui permettent de s’assurer du caractère aléatoire des
nombres générés par un algorithme numérique [1], nous avons choisi le test spectral, non
seulement parce qu’il peut être illustré graphiquement, mais parce qu’il apparait comme
l’un des test les plus forts. Un mauvais algorithme peut passer avec succès un certain
nombre de test statistiques, excepté celui-ci.
Le test spectral examine les correlations entre k nombres consécutifs {X n+1 , Xn+2 , · · · , Xn+k },
fournis par un algorithme.
Même si un algorithme de congruence linéaire fournit des nombres individuellement “quasi-
aléatoires”, les k − uples de nombres consécutifs sont loins d’être exempts de corrélations !
Si {Xn+1 , Xn+2 , · · · , Xn+k } représentent les coordonnées d’un point P n dans un espace à k
dimensions, alors les points Pn se répartissent sur des hyperplans à n − 1 dimensions et on
montre [3] qu’il y a au plus m1/k tels plans (mais il peut y en avoir beaucoup moins !...).
Reprenons le même exemple de la Figure 3.1a, et reportons sur un plan les coordonnées
des doublets (1,6), (6,15), (15,2) etc... . Les points obtenus se répartissent sur 4 lignes (Fig.
3.2). Puisque m = 16 = 24 , on en attend au plus 24/2 = 4 lignes.
Si on considère, non pas la suite des entiers X n , mais les réels Un = Xn /m appartenant
à [0, 1[, pour les k − uples consécutifs P n = (Un+1 , Un+2 , · · · , Un+k ) représentés par des
points d’un hypercube à k dimensions de coté [0,1], la distance maximale entre hyperplans

15.0

10.0

5.0

0.0
0.0 5.0 10.0 15.0

Fig. 3.2 – Corrélations entre paires pour la congruenc linéaire : X n+1 = 5Xn + 1 [mod 16].
24 Techniques d’échantillonnage

notée 1/νk par Knuth constitue une mesure des corrélations.

Nous donnons maintenant un exemple non trivial de générateur de nombres quasi-aléatoires,

qui a été utilisé pendant plusieurs décennies sur des machines 32 bits (en particulier IBM) et
qui subsistait encore comme l’un des sous programmes intrinsèques de génération de nombre
aléatoires dans la version 1998 de “DIGITAL FORTRAN90” ! Son nom est “RANDU()”. Il
utilise un algorithme de congruence linéaire avec c = 0, a = 2 16 + 3 et m = 232 . Il est rapide
puisque l’opération [mod m] s’effectue par une opération binaire AND et la multiplication
par a peut être implémentée de manière efficace puisque (a − 3) est une puissance de 2.
Par contre, les corrélations entre k-uplets (X n , Xn+1 , · · · , Xn+k ) sont catastrophiques. La
Figure 3.3 illustre les corrélations entre triplets. L’utilisation de ce générateur de nom-
bres quasi-aléatoires pour le calcul Monte-Carlo de l’aimantation dans le modéle d’Ising
ferromagnétique à deux dimensions donne une valeur 10% plus forte que la valeur exacte !

Même pour des générateurs de nombres quasi-aléatoires relativement satisfaisant, des cor-
rélations existent, mais à une échelle plus petite. L’une de bibliothèques mathématiques
les plus utilisées : IMSL, propose un algorithme de congruence linéaire avec c = 0, a = 7 5
et m = 231 − 1. Ici, a est suffisamment faible pour que la multiplication par a puisse être
implémentée de manière efficace. Les corrélations entre triplets (Fig. 3.4) n’apparaissent
qu’à une échelle 100 fois plus faible.

Knuth donne un tableau des νk , k = 2 à 6 pour une trentaine d’algorithmes dont ceux

Fig. 3.3 – Corrélations entre triplets représentées dans un cube [0,1]x[0,1]x[0,1] pour la
congruence linéaire Xn+1 = (216 + 3)Xn [mod 232 ]
3.1 Génération de nombres “pseudo-aléatoires” 25

Fig. 3.4 – Corrélations entre triplets représentées dans un cube [0,0.01]x[0,0.01]x[0,0.01]

pour la congruence linéaire Xn+1 = 75 Xn [mod 231 − 1].

Fig. 3.5 – Corrélations entre triplets représentées dans un cube [0,0.01]x[0,0.01]x[0,0.01]

pour la congruence linéaire Xn+1 = 75 Xn [mod 231 − 1] après “shuffling”.
26 Techniques d’échantillonnage

décrits ci-dessus.

Pour briser ces corrélations, une technique appelée “shuffling” (traduction : “battre les
cartes”) peut être utilisée. L’algorithme correspondant est simple [1, 3] et consomme peu
de temps CPU. Il consiste à garder constamment en mémoire un tableau des derniers p
nombres tirés et les permuter aléatoirement en fonction de nouveaux nombres tirés au
hasard par un autre algorithme.
La bibliothèque IMSL propose l’algorithme précédemment cité avec shuffling superposé. Le
coût supplémentaire en temps CPU est de l’ordre de 50%, mais les corrélations (Fig. 3.5)
disparaissent.

3.1.4 Suites de Fibonacci “retardées”

Dans les congruences considérées précédemment, X n ne dépend que du nombre précédent

Xn−1 . On peut accroı̂tre la période en choisissant une relation où X n dépend de plus d’une
valeur précédente. La relation la plus simple, à cet effet, est :

Xn+1 = (Xn + Xn−1 ) [mod m] (3.4)

En choisissant X0 = X1 = 1, la récurrence Xn+1 = (Xn−1 + Xn−1 ) définit la “suite de
Fibonacci” et on sait que lorsque n → ∞ le rapport X n+1 /Xn de deux termes successifs
tend vers le nombre d’or. Elle ne peut conduire à un bon générateur de nombres aléatoires !
Par contre on peut généraliser :

Xn+1 = (Xn−l + Xn−k ) [mod m] (3.5)

Les nombres k et l sont appelés “lags” (i.e “retards” ou “décalages”) On montre que si
m = 2α est une puissance de deux, et si l et k sont tels que le trinome :

xk + x l + 1

est un polynome premier (i.e. n’a aucun autre polynome diviseur que lui-même ou une
constante) dans le corps des entiers modulo 2, alors la période est 2 α−1 (2k − 1) (k > l) et
le nombre de cycles différents est 2 (k−1)×(α−1)
Les couples de “lags” (24,55), (38,89), (37,100), (30,127), (83,258) satisfont cette pro-
priété[1, 5].
Pour k et l suffisament grands, tous les test statistiques s’avèrent en général excellents, bien
supérieurs aux résultats des congruences linéaires. Cette méthode présente de nombreux
avantages :
– Calcul rapide, puisqu’il n’y a pas de multiplication et si m = 2 α , l’opération “modulo
m” s’effectue par un cache et une opération binaire AND.
– La période peut être très longue.
3.1 Génération de nombres “pseudo-aléatoires” 27

– Le nombre de cycles indépendants (2 (k−1)×(α−1) ) est extrêmement grand, ce qui en fait

actuellement une méthode de choix pour implantation sur ordinateurs massivement par-
allèles.
G.J Mitchell and D.P. Moore ont, les premiers proposé ce type d’algo-rithme[4, 1] avec
m pair, p = 24, q = 55 et {X0 , X1 , X2 , · · · , X54 } entiers arbitraires, non tous pairs. Pour
m = 2α sa période est 2α−1 (255 − 1). Il figure aussi parmi les générateurs de nombres
aléatoires proposé dans “Numérical Recipes”[3].

Au lieu de la relation additive précédente, on peut préférer une relation multiplicative :

Xn+1 = (Xn−l × Xn−k ) [mod m] (3.6)

Les performances statistiques sont encore meilleures, mais au prix d’un accroissement du
temps CPU.

Enfin, on peut encore utiliser des algorithmes hybrides en associant cette méthode et celle
des congruences linéaires :

Xn+1 = (a1 Xn−1 + a2 Xn−2 + · · · + ak Xn−k ) [mod m] (3.7)

Avec m = 2, l’équation précédente peut fournir un excellent générateur de bits (0 ou 1)

quasi-aléatoires. Dans ce cas, les coefficients a i sont 0 ou 1 et le polynome

1 + a 1 x + a 2 x2 + · · · + a k xk

doit être un polynome premier dans le corps des entiers modulo 2 (cf. les travaux de
Tausworthe[6])

3.1.5 Puis-je utiliser la fonction “RAND(), random(), ...” de mon ordi-

nateur ?

La réponse dans les années 80 était dans la plupart des cas négative et à cette époque, bon
nombre de calculs stochastiques publiés étaient erronnés par suite de l’usage de générateurs
de nombres peu aléatoires !
Les choses on maintenant évolué. Sur les stations 32 bits, la plupart des constructeurs
fournissent un choix de générateurs de nombre quasi-aléatoires dont certains sont satis-
faisants (mais il faut encore se méfier -cf. fonction RANDU(), qui, malgré ses performances
pitoyables n’a pas encore complètement disparu-).
Quelquefois on a accès à une notice détaillée sur le ou les générateurs de nombres aléatoires
proposés avec description de l’algorithme utilisé et résultats de test statistiques (ceux
préconisés par Knuth).
28 Techniques d’échantillonnage

Dans le cas de générateurs “boite noire”, il est prudent de se livrer soi même à quelques
tests (ex : test spectral) avant de faire un choix.
L’avantage à utiliser un générateur de nombre aléatoires fourni par le constructeur est
que l’algorithme est programmé en langage machine en tenant compte de son architecture
spécifique. Le même algorithme, programmé en langage évolué (Fortran, C...) serait plus
lent.
Par contre, si on souhaite qu’un programme soit portable et donne exactement les mêmes
résultats, quelle que soit la machine, on peut programmer soi-même l’algorithme.
Les bibliothèques IMSL et NAG contiennent actuellement de bon sous-programmes de
génération de nombres aléatoires.

Tous ces programmes fonctionnent de la manière suivante :

– Si on ne précise pas “la semence”, c’est à dire le ou les nombres initiaux pour démarrer
l’algorithme, l’ordinateur le fait lui même. Soit il utilise toujours les mêmes nombres
(dans ce cas, si on relance le programme on obtient le même résultat), soit il les calcule
par un algorithme, en prenant pour point de départ un nombre formé à partir de la date
de son horloge (dans ce cas on a un résultat différent à chaque exécution).
– Il y a toujours possibilité de choisir les paramètres initiaux et les passer à l’ordinateur
par un sous programme prévu à cet effet.
Nous terminons par quelques comparaisons de temps de calculs obtenus pour divers générateurs
de nombre aléatoires sur un processeur Pentium IV cadencé à 3 GHz

Temps CPU pour le tirage de 109 nombres pseudo-aléatoires.

(i) Bibliothèque mathématique IMSL, congruence linéaire avec c = 0, m = 2 31 −1, a = 75 :

50 s CPU.
(ii) Bibliothèque mathématique IMSL, méthode précédente + “shuffling” : 75 s CPU.
(iii) Bibliothèque mathématique IMSL, congruence linéaire avec c = 0, m = 2 31 − 1,
a = 397204094 : 100 s CPU.
(iv) Fibonacci : Xn = (Xn−24 + Xn−50 ) [mod 109 ], programme Fortran de “Numerical
Recipes”[3] : 20 s CPU.

3.2 Echantillonnage d’une loi de probabilité non uniforme

3.2.1 Méthode de transformation

Changement de variable

Soit X une variable aléatoire réelle possédant la densité de probabilité f (x). Soit Ψ(x) une
fonction réelle, monotone croissante, continue et dérivable.
3.2 Echantillonnage d’une loi de probabilité non uniforme 29

A X nous associons la variable aléatoire Y définie par :

Y = Ψ(X)

Nous allons montrer que la variable aléatoire Y possède une densité de probabilité g(y) que
nous allons expliciter.

D’après la définition de la densité de probabilité :

P ({x < X < x + dx}) = f (x)dx pour dx → 0

mais puisque Ψ(x) est une fonction monotone croissante, elle réalise une bijection et le
premier membre de cette égalité peut aussi s’ecrire :
dΨ(x)
P ({Ψ(x) < Y < Ψ(x + dx)}) = P ({Ψ(x) < Y < Ψ(x) + dx})
dx
dΨ(x)
En notant y = Ψ(x) et dy = dx dx

P ({y < Y < y + dy}) = f (x)dx

Puisque Ψ est monotone, elle admet une fonction réciproque Ψ −1 , donc

x = Ψ−1 (y)

On a
dy dΨ−1 (y)
dx = dΨ(x)
= dy
dy
dx
d’ou le résultat final
dΨ−1 (y)
P ({y < Y < y + dy}) = f (Ψ−1 (y)) dy
dy
Ce qui prouve que Y admet une densité de probabilité g

dΨ−1 (y)
g(y) = f (Ψ−1 (y))
dy

Dans le cas d’une fonction monotone décroissante, on obtient :

dΨ−1 (y)
g(y) = −f (Ψ−1 (y)) (3.8)
dy

Ce type de transformation peut se généraliser pour des variables aléatoires à valeurs dans
Rd . Elle fait intervenir un Jacobien.
30 Techniques d’échantillonnage

Cas particulier où ψ(x) = F (x), F (x) représentant la fonction de répartition

d’une variable aléatoire X de densité de probabilité f (x)

Z x
F (x) = f (u)du
−∞
Par construction, la fonction de répartition F est positive, monotone, croissante, de dérivée :
dF (x)
= f (x) (3.9)
dx
Elle varie entre F (−∞) = 0 à F (∞) = 1

La function F (x) admet donc pour 0 ≤ y ≤ 1 une fonction inverse F −1 (y) ayant pour
dérivée
dF −1 (y) 1 1
= = −1
dy f (x) f (F (y))
Et la densité de probabilité g(y) de la variable aléatoire Y = F (X) s’écrit donc :

f (F −1 (y))
g(y) = =1 0≤y≤1
f (F −1 (y))

Elle correspond à la loi uniforme sur le segment [0,1].

Echantillonnage d’une densité de probabilité quelconque f

– Etant donnée
Rx une densité de probabilité f (x), on détermine la fonction de répartition
F (x) = −∞ f (u)du.
– On tire des nombres η au hasard suivant la loi uniforme sur le segment [0,1]
– D’après ce qui précède les nombres x = F −1 (η) sont répartis suivant la densité de prob-
abilité f .

Quelques exemples

loi exponentielle
La loi exponentielle de densité de probabilité :

f (x) = λ exp(−λx) 0≤x<∞

intervient souvent en physique, elle régit le temps qui sépare deux désintégrations succes-
sives d’un élément radioactif. Sa fonction de répartition est :
Z x
F (x) = λ exp(−λu)du = − [exp(−λu)] u0 = 1 − exp(−λx)
0
3.2 Echantillonnage d’une loi de probabilité non uniforme 31

La fonction réciproque F −1 est donnée par l’inversion de la relation :

y = 1 − exp(−λx)

en
exp(−λx) = 1 − y
soit
x = −Ln(1 − y)/λ
Donc on tire un nombre η uniformément dans l’intervalle [0, 1] et les nombres x = −Ln(1 −
η)/λ sont répartis suivant une loi exponentielle.

loi gaussienne
La loi gaussienne a pour densité :

x2
 
1
f (x) = √ exp − 2
2πσ 2 2σ
sa fonction de répartition :
x
u2
 
1
Z
F (x) = √ exp − 2
2πσ 2 −∞ 2σ

s’exprime à l’aide de la fonction “erreur”, qui elle même n’admet pas de fonction inverse
analytique. On va donc devoir utiliser une astuce.

Il s’avère qu’il est plus simple de tirer un couple (X 1 , X2 ) de deux variables aléatoires
gaussiennes indépendantes.
Puisque X1 et X2 sont deux variables aléatoires indépendantes, la densité de probabilité
du couple (X1 , X2 ) est le produit f (x1 )f (x2 ) des densités de probabilité de chacune des
deux variables aléatoires. C’est à dire que :
 2
x1 + x22 dx1 dx2

P ({x1 < X1 < x1 + dx1 } ∩ {x2 < X2 < x2 + dx2 }) = exp −
2σ 2 2πσ 2
En coordonnées polaires :
x1 = ρ cos θ
x2 = ρ sin θ
dx1 dx2 = ρdρdθ
le second membre s’écrit :
ρ2
 
ρ dθ
exp − 2 2
dρ
2σ σ 2π
c’est à dire que pour échantillonner deux variables aléatoires gausiennes indépendantes, on
peut :
32 Techniques d’échantillonnage

– Tirer un angle théta uniformément sur le segment [0, 2π]

– Tirer un rayon ρ suivant la loi de probabilité

ρ2
 
ρ
h(ρ) = exp − 2
2σ σ 2

les deux nombres x1 = ρ cos θ et x2 = ρ sin θ seront répartis suivant des distributions
gaussiennes indépendantes.
Echantillonner ρ suivant la loi de probabilité h(ρ) est facile. La fonction de répartition
correspondante est :
Z ρ ρ
u2 u2 ρ2
     
u
F (ρ) = exp − 2 du = − exp − 2 = 1 − exp − 2
0 2σ σ2 2σ 0 2σ

Son inverse est obtenue en écrivant :

ρ2
 
F (ρ) = r = 1 − exp − 2
2σ

d’où p
ρ= −2σ 2 Ln(1 − r)
Donc on tire r uniformément sur le segment [0, 1] et le nombre ρ obtenu par la transforma-
tion précédente apparait avec la densité de probabilité h(ρ).
Pour les cas où la fonction de répartition n’admet pas d’inverse qui puisse s’exprimer de
manière analytique, il nous reste une méthode qui coûte plus cher en tant de calcul, mais
qui s’applique de manière tout à fait générale.

3.2.2 Méthode de réjection de Von Neumann

Cette méthode élémentaire est illustrée sur la Figure 3.6.

Soit à échantillonner la densité de probabilité proportionnelle à la fonction f (x) représentée
par le trait plein sur le segment [xmin, xmax] et égale à 0 en dehors. Soit C une constante
supérieure ou égale à F max, valeur maximale de f (x). La méthode consiste à engendrer N
couples de réels, uniformément dans le pavé [xmin, xmax] × [0, C] et à ne retenir que ceux
qui “tombent” sous la courbe f (x).
Soit M le nombre total de couples retenus et ν M (x)dx le nombre de ceux d’entre eux dont
l’abscisse est comprise entre x et x + dx (dx fixé, arbitrairement petit). Alors lorsque N
tend vers l’infini le rapport

νM (x)dx ˜ f (x)dx
tend vers : f(x)dx = R xmax
M xmin f (x)dx

Concrètement :
– La valeur une valeur de x est “tirée” suivant une loi uniforme sur le segment [xmin, xmax]
3.2 Echantillonnage d’une loi de probabilité non uniforme 33

Fmax
f(x)

xmin xmax
Fig. 3.6 – Méthode de réjection de Von Neuman

– Un nombre aléatoire η distribué uniformément dans le segment [0,1] est généré.

– La valeur x est acceptée si f (x)/C > η ; sinon, elle est rejetée.
La relation ci-dessus montre que la méthode peut permettre plus généra-lement d’engendrer
˜
une distribution de probabilité f(x), proportionnelle à une fonction positive quelconque
f (x), même lorsque l’intégrale de cette dernière n’est pas normalisée à un.
La méthode de réjection est évidemment plus coûteuse en temps de calcul que la méthode
de transformation, puisqu’une partie des tirages est simplement rejetée !
Lorsque F max est connu, on minimisera le nombre de réjections en prenant C = F max.
34 Techniques d’échantillonnage
 Echantillonnage suivant
4 l’importance.

4.1 Calcul stochastique d’une intégrale

Considérons une intégrale à d dimensions :

Z
I = dxΦ(x) (4.1)

où x = (α1 , α2 , · · · , αd ) est un vecteur à d dimensions et

Z Z Z Z
dx = · · · dα1 dα2 · · · dαd

Φ(x) est une fonction de R d dans R.

Choisissons arbitrairement une fonction f (x) ayant les propriétés d’une densité de proba-
bilité, c’est à dire positive et telle que
Z
dxf (x) = 1

En définissant :
Φ(x)
g(x) =
f (x)
L’intégrale sécrit : Z
I= g(x)dP avec dP = f (x)dx (4.2)

Elle représente la moyenne stochastique (ou “Espérance”) de la variable aléatoire Y = g(X)

transformée de la variable aléatoire X par la fonction g

Si nous “tirons” N vecteurs xi de Rd suivant la loi de probabilité P , nous pouvons associer

au résultat de cette expérience la variable aléatoire :
SN Y1 + Y 2 + · · · + Y N g(X1 ) + g(X2 ) + · · · + g(XN )
= =
N N N
36 Echantillonnage suivant l’importance.

D’après la loi forte des grands nombres (Annexe A), la “moyenne expéri-mentale” :
y1 + y 2 + · · · + y N g(x1 ) + g(x2 ) + · · · + g(xN )
= avec yi = g(xi )
N N
converge “presque sûrement” vers I lorsque N → ∞.
Si σg2 représente la variance de la variable aléatoire g(X) :
Z  2
Φ(x)
Z
2 2
σg = [g(x) − I] f (x)dx = − I f (x)dx (4.3)
f (x)
alors la variance de la variable aléatoire S N /N est σg2 /N (cf. Annexe A), et l’écart type
√ √
σg / N qui caractérise l’erreur décroı̂t en 1/ N .

4.2 Réduction de la variance

Nous allons maintenant mettre à profit le caractère arbitraire dans le choix de la densité
de probabilité f pour réduire la variance, donc l’erreur.
Il est toutefois illusoire d’essayer d’annuller la variance. Ceci ne peut se faire que si Φ(x)/f (x) =
I dans la relation précédente, mais pour cela il faut connaı̂tre I, c’est à dire avoir résolu le
problème !

Nous allons d’abord illustrer l’optimisation du choix de f sur un exemple à une dimen-
sion et nous reprendrons celui proposé par Hammersley et Handscomb dans leur ouvrage
introductif aux méthodes de Monte-Carlo[7].

Considérons la fonction dans R :

ex − 1
Φ(x) = 2
e−1
variant de 0 à 2 lorsque x varie de 0 à 1 et représentée sur la figure 4.1.
Nous nous proposons de calculer par la méthode de Monte-Carlo l’intégrale :
Z 1
I= Φ(x)dx
0

Le premier choix qui vient à l’esprit pour la densité de probabilité f u (x) est la loi uniforme
sur le segment [0, 1] :
fu (x) = 1 si x ∈ [0, 1] et f (x) = 0 sinon.
La fonction choisie est intégrable analytiquement, on a I=0.836... et l’écart-type :
s
Z 1
u
σg = (Φ(x) − I)2 dx = 0.5726
0
4.2 Réduction de la variance 37

3.0

φ(x)

2.0

fl(x)

fu(x)
1.0

0.0

−1.0
−0.5 0.0 0.5 1.0 1.5

Fig. 4.1 – Trait plein : fonction Φ(x) à intégrer entre 0 et 1. Tirets : densité de probabilité f u
uniforme. Traits mixtes : une densité de probabilité f l simple réalisant un échantillonnage
suivant l’importance.

On en d’eduit, par exemple, que pour obtenir I à 1% près, dans un “intervalle

√ de confiance”
à 95%, on doit sommer sur un nombre de points N défini par 2σ/( N I) = 0.01 soit
N = 18800 points

Ce premier choix n’est pas très astucieux, car en “tirant” des points avec une loi uniforme,
on a la même densité de points vers x = 0 où la fonction est faible que vers x = 1 où elle
est maximale. On a plutôt intérêt à choisir une loi de probabilité qui donne une densité de
points “tirés” plus importante là où la fonction est maximale, d’où l’idée “d’échantillonnage
suivant l’importance”.

Prenons la densité de probabilité linéaire (cf. Fig. 4.1) :

fl (x) = 2x si x ∈ [0, 1] et f (x) = 0 sinon.

qui est beaucoup plus proche de la fonction Φ(x) que le choix précedent. Alors l’écart type
est : s
Z 1
l Φ(x)
σg = ( − I)2 2xdx = 0.1048
0 2x
38 Echantillonnage suivant l’importance.

soit 5.49 fois plus faible. On obtiendra donc la même précision avec 5.49 2 = 30 fois
moins de points soit seulement N = 630. Cet exemple à une dimension reste toutefois
académique, dans la mesure où une méthode de Gauss à 5 points donne une précision
nettement supérieure !

L’échantillonnage suivant l’importance devient par contre spectaculaire pour des intégrales
à grand nombre de dimensions. Pour rester dans le cadre de calculs que l’on peut encore
faire analytiquement, généralisons l’exemple précédent à d dimensions en prenant sur R d
(cf. Negele et Orland[8]) :

Ψ(α1 , α2 , · · · , αd ) = Φ(α1 )Φ(α2 ) · · · Φ(αd )

que nous intégrons sur un hypercube [0, 1] × [0, 1] × · · · × [0, 1] Et prenons pour densité de
probabilité le produit :
f (α1 )f (α2 ) · · · f (αd )
Du fait de la séparabilité de l’intégrale, on obtient aisément la valeur de l’intégrale à d
dimensions : Id = I d ainsi que la variance :

d Z 
" d Z #2
Φ(αi ) 2
Y  Y
σd2 = f (αi )dαi − f (αi )dαi
f (αi )
i=1 i=1

soit :
σd2 = (σg2 + I 2 )d − I 2d

et en divisant chaque membre par I 2d

σd2 σg2
= ( + 1)d − 1
Id2 I2

avec σgu = 0.5726 pour la loi uniforme on obtient, pour d grand :

σdu
≈ 1.4691d/2
Id

avec σgl = 0.1048 pour un echantillonnage suivant l’importance (loi linéaire) on a, pour d
grand
σdl
≈ 1.0157d/2
Id
Pour d = 100, ce qui reste modeste pour un calcul Monte-Carlo, on améliore l’écart-type,
donc l’erreur statistique, d’un facteur

(1.469/1.016)50 ≈ 108

par l’échantillonnage suivant l’importance.

4.2 Réduction de la variance 39

Pour fixer les idées, pour obtenir l’intégrale à d=100 dimensions, dans un intervalle de
confiance à 95%, avec un précision relative de 1%, il nous faut, avec la loi uniforme un
nombre N tel que :
2 ∗ 1.46950
√ = 0.01
N
soit N = 2 ∗ 1021 points. Et pour la même précision, il ne nous faudra avec la loi linéaire
que N = 2 ∗ 1021 /1016 = 2 ∗ 105 points
Pour fixer les idées, en supposant qu’avec un Pentium 4, cadencé à 3 GHz, le “tirage” d’un
point x et le calcul de Ψ(x) nécessite 300 cycles. L’échantillonnage avec N = 2 × 10 21
points nécessiterait 2 × 1014 s, soit six millions d’années.
Avec un échantillonnage suivant l’importance élémentaire, on obtient le même résultat en
vingt millisecondes !...
40 Echantillonnage suivant l’importance.
 Simulation de systèmes
5 statistiques

5.1 Introduction

Dans le cadre de l’ensemble “canonique”, considérons un système S en contact avec un

réservoir d’énergie (ou thermostat) qui impose sa température T . Supposons d’abord que
le système est décrit par un ensemble d’états discrets. La probabilité pour le système S de
se trouver dans un état l d’énergie E l est :
exp[−El /(kB T )]
Pl =
Z
où X
Z= exp[−El /(kB T )]
l
est la fonction de partition du système.
La valeur moyenne de l’energie est donc :
P
El exp[−El /(kB T )]
< E >= l
Z
et plus généralement, la valeur moyenne de toute grandeur O l associée à l’état l est :
P
Ol exp[−El /(kB T )]
< O >= l
Z
Si nous passons à la limite continue : l’état du système S est caractérisé par un ensemble
de variables que nous noterons collectivement x, x ∈ R d représentant un vecteur à d
dimensions (par exemple l’ensemble des coordonnées x = (x 1 , x2 , ..., xn ) des particules qui
le constituent). L’energie correspondant à l’état x est E(x). La probabilité de trouver le
système dans un état dont les variables x sont comprises entre x et x + dx est :
exp[−E(x)/(kB T )]
dx
Z
La fonction de partition Z s’exprimant alors comme une intégrale multiple :
Z
Z = exp[−E(x)/(kB T )]dx
42 Simulation de systèmes statistiques

La valeur moyenne de toute grandeur O(x) associée à l’état x (par exemple, l’énergie
O(x) = E(x) est définie par :
R
O(x) exp[−E(x)/(kB T )]dx
< O >= (5.1)
Z
L’intérêt de la méthode de Monte-Carlo, avec “échantillonnage suivant l’importance” pour
le calcul de l’intégrale multidimensionnelle 5.1 est immédiat :
– Si on est capable d’engendrer un ensemble de configurations x correspondant à une
variable aléatoire X, à valeurs dans R d ayant pour densité de probabilité :

exp[−E(x)/(kB T )]
(5.2)
Z
– Alors la valeur moyenne < O > de la grandeur physique O(x) apparaı̂t comme l’espérance
de la variable aléatoire O(X).
Le problème crucial est que l’on ignore dans l’expression de la densité de probabilité 5.2 la
valeur de la constante de normalisation Z, c’est à dire la fonction de partition ! On pourrait
éventuellement s’en sortir en utilisant la méthode de réjection de Von Neumann, mais dans
ce problème à très grand nombre de dimensions, cette méthode s’avère tout à fait inefficace.
Une méthode élégante et efficace a été proposée par Métropolis. Elle construit un ensemble
de configurations suivant la densité de probabilité 5.2, à partir des états successifs d’une
chaı̂ne de Markov ayant atteint sa distribution d’équilibre.

5.2 Chaı̂nes de Markov

5.2.1 Définition

Une suite infinie ordonnée (X1 , X2 , · · · , Xt , · · · ) de variables aléatoires constitue une chaı̂ne
de Markov si la loi de probabilité conditionnelle de X t+1 lorsqu’on se donne les valeurs de
X1 , X2 , ...,Xt , se réduit à la loi de probabilité conditionnelle de X t+1 lorsqu’on se donne
seulement la valeur de Xt .
[t désigne ici un entier naturel, faisant référence à un temps discrétisé]

Remarque importante : pour une suite ainsi définie, la connaissance du passé influe sur
l’avenir ; mais, d’un ensemble d’événements passés, seule subsiste l’influence du plus récent.
Si l’on se donne Xt , la connaissance de Xt0 pour t0 > t n’est en rien précisée par la donnée
des valeurs de Xt−1 , Xt−2 , ...,X1 .

Nous nous restreindrons aux Chaı̂nes de Markov “homogènes” où la loi de probabilité
conditionnelle de Xt+1 lorsqu’on se donne la valeur de Xt ne dépend pas de t. Alors la
chaı̂ne de Markov est entièrement déterminée par cette loi de probabilité conditionnelle
appelée “loi de transition” et par la loi de probabilité initiale de la variable aléatoire X 1 .
5.3 Convergence vers une loi de distribution invariante imposée 43

Nous supposerons aussi que la loi de probabilité conditionnelle de X t+1 lorsqu’on se donne
la valeur de Xt possède la densité de probabilité p(y/x),

P ({y < Xt+1 < y + dy}/{Xt = x}) = p(y/x)dy (5.3)

Et nous noterons :
p(y/x) = p(x → y)
Comme toute densité de probabilité, la fonction p(y/x) = p(x → y) satisfait quel que soit
x à la condition de normalisation :
Z Z
p(y/x)dy = p(x → y)dy = 1 (5.4)

Au cours des rappels de théorie des probabilités (Annexe B), nous avons montré, dans le cas
de chaı̂nes de Markov à valeurs dans un ensemble discret, que lorque la chaı̂ne de Markov
possède la propriété d’ergodicité, c’est à dire que pour tout couple d’états discrets (a i , aj ),
il existe une probabilité non nulle de passer de l’un à l’autre après n pas, alors cette chaı̂ne
converge vers une loi de distribution invariante unique.
Nous admettrons que ce résultat reste vrai dans le cas continu : si pour tout couple d’état
(x, y) il existe une probabilité non nulle de passer de l’un à l’autre en un nombre fini
d’étapes n, alors il y a convergence vers une loi de distribution invariante unique.
Il nous reste alors à imposer une condition à la loi de transition p(y/x) = p(x → y) pour
que cette loi de distribution invariante limite soit :

exp[−E(y)/(kB T )]
f (y) =
Z

5.3 Convergence vers une loi de distribution invariante im-

posée

Théorème
– Si “la loi de transition” p(x → y) est ergodique, c’est à dire que pour tout couple d’états
(x, y) il existe une probabilité non nulle de passer de l’un à l’autre en n étapes [–donc
tous les états sont “visités”–]
– si p(x → y) satisfait à la condition “de microréversibilité”

[f (x)p(x → y) − f (y)p(y → x)] = 0 (5.5)

Cette condition de microréversibilité est encore appelée principe du bilan détaillé par
réference à l’équation maı̂tresse et au principe du bilan détaillé qui en découle en ther-
modynamique (cf. Annexe C)
Alors la loi de distribution de la chaı̂ne de Markov correspondante converge vers une dis-
tribution invariante qui est proportionnelle à f (x).
44 Simulation de systèmes statistiques

Preuve :
Puisque la chaı̂ne de Markov est ergodique, nous savons que sa loi de distribution converge
vers une distribution invariante et que cette distribution invariante est unique. Il nous suffit
donc de montrer que la distribution f (x) est invariante par rapport à la loi de transition
p(x → y).

Si “à l’instant t”, Xt est distribuée suivant la loi π t (x) = f (x), alors, en vertu de l’équation
de Chapman-Kolmogorov B.3, “à l’instant (t + 1)”, X t+1 est distribuée suivant la loi de
probabilité :
Z
π t+1 (x) = duf (u)p(u → x)

Si p(u → x) satisfait la relation du bilan détaillé 5.5, la relation précédente peut encore
s’écrire : Z Z
t+1
π (x) = duf (x)p(x → u) = f (x) dup(x → u)

et le second membre est simplement :

π t+1 (x) = f (x)

R
puisque dup(x → u) = 1 ∀x

5.3.1 Echantillonnage de f (x) = e−βE(x) /Z

Etant donnée une chaı̂ne de Markov de loi de transition p(y → x), nous allons “échantillonner”
une suite de points {x1 , x2 , · · · , xt , · · · } (xi ∈ Rd ) suivant cette chaı̂ne, i.e. : xt ayant été
obtenu, le point suivant xt+1 est “tiré” suivant la loi de probabilité p(x t → xt+1 ).

Si p(xt → xt+1 ) satisfait la relation 5.5 du bilan détaillé, avec f (x) = e −βE(x) /Z, c’est à
dire : " #
e−βE(xt ) e−βE(xt+1 )
p(xt → xt+1 ) − p(xt+1 → xt ) = 0
Z Z

alors après un certain “temps d’équilibre”, c’est à dire pour n > N (N suffisamment grand)
les xt sont distribués suivant la densité de probabilité e −βE /Z.

Remarque importante : on peut multiplier chaque membre de l’équation ci-dessus par Z :

h i
e−βE(xt ) p(xt → xt+1 ) − e−βE(xt+1 ) p(xt+1 → xt ) = 0

La connaissance de la constante de normalisation (ici fonction de partition Z) devient

absolument inutile pour l’échantillonnage suivant la loi de probabilité e −βE /Z.
5.4 Algorithme de Metropolis 45

5.4 Algorithme de Metropolis

L’algorithme le plus utilisée pour engendrer une chaı̂ne de Markov qui satisfasse le principe
du bilan détaillé a été proposé par Metropolis et al.[9].

Partant d’une configuration xt = x, d’energie E(x), on définit une configuration d’essai

y en déplaçant une ou plusieurs particules et on calcule sa nouvelle energie E(y). Cette
configuration d’essai est acceptée pour nouvelle configuration avec la probabilité p(x →
y) = min{1, e−βE(y) /e−βE(x) }, sinon, elle est rejetée et on prend pour nouvelle configuration
xt+1 = xt = x c’est à dire :
– xt+1 = y avec la probabilité p = min{1, e−βE(y) /e−βE(x) }
– xt+1 = xt = x avec la probabilité (1 − p)
En pratique :
– Si E(y) ≤ E(x), alors xt+1 = y, la configuration d’essai est acceptée avec certitude
(p = 1)
– Si E(y) > E(x), alors on “tire” un nombre η suivant la loi uniforme dans le segment [0,1]
et on accepte pour nouvelle configuration la configuration d’essai y si e −βE(y) /e−βE(x) >
η, sinon la nouvelle configuration est identique à l’ancienne : x t+1 = xt = x.
Cette technique simple, similaire à la méthode de réjection de Von Neuman décrite précédem-
ment, impose le fait que la configuration d’essai est acceptée avec la probabilité p(x → y) =
min{1, e−βE(y) /e−βE(x) }.

Il nous reste à verifier que cette probabilité p(x → y) satisfait effectivement l’équation du
bilan détaillé.

Supposons E(y) > E(x), alors :

p(y → x) = 1 et p(x → y) = e−βE(y) /e−βE(x)

la relation
p(y → x)e−βE(y) = p(x → y)e−βE(x)
est bien vérifiée.
La vérification pour le cas E(y) ≤ E(x) est analogue.
46 Simulation de systèmes statistiques
 6 Le modèle d’Ising

Nous allons “expérimenter” l’algorithme de Metropolis sur un modèle simple, dont les
propriétés thermodynamiques sont déjà d’une extrême richesse.

6.1 L’Hamiltonien d’Ising

On considère un ensemble de spins s i , porteurs d’un moment magnétique, qui ne peut

prendre que deux valeurs :
– +1 si le moment magnétique est dirigé “vers le haut”
– −1 si le moment magnétique est dirigé “vers le bas”
Deux spins premiers voisins subissent une interaction J qui tend à les aligner. A température
nulle, les spins seront donc tous alignés pour minimiser l’énergie : ceci constitue le modèle
le plus simple pour un corps ferromagnétique. En l’absence de champ magnétique, l’état
fondamental du système est doublement dégénéré car les spins peuvent être alignés tous
vers le haut ou tous vers le bas. Si on applique une induction magnétique B dirigée vers
le haut, par exemple, on lève la dégénérescence : tous les spins s’orientent vers le haut.
L’Hamiltonien s’écrit :
1erX
vois. X
H = −J sα sβ − B sα (6.1)
<α,β> α

La première somme est effectuée sur toutes les paires distinctes de premiers voisins.
Nous nous limiterons, en Travaux Pratiques au cas où l’induction magnétique est nulle
B = 0.

6.2 Modélisation par la méthode de Monte-Carlo

On considère un réseau carré de taille L × L comprenant un nombre N = L 2 de sites. Les

sites sont repérés par leurs coordonnées α = (i, j), i et j variant de 0 à L − 1. On applique
des conditions aux limites périodiques :
s(i+pL,j+qL) = s(i,j)
48 Le modèle d’Ising

On travaille dans l’ensemble canonique à une température donnée T .

On utilise l’algorithme de Metropolis en envisageant des mouvement élémentaires locaux
(retournement d’un spin).

1. On choisit une configuration initiale, par exemple :

– Tous les spins alignés
– Spins aléatoirement distribués vers le haut ou vers le bas
2. On choisit un spin, au hasard (suivant une loi de distribution uniforme)
3. On envisage comme tentative de mouvement le retournement de ce spin
4. On calcule la différence d’énergie d’interaction des spins ∆ entre la nouvelle configu-
ration avec un spin retourné et la configuration de départ.
5. Si ∆ est négatif, ou nul c’est à dire, si retourner un spin diminue l’energie, ou la laisse
invariante, on accepte la nouvelle configuration
6. Si ∆ est positif :
– On tire un nombre η, au hasard, suivant une loi uniforme sur le segment [0, 1]
– Si η < exp(−∆/kT ) on accepte la configuration avec spin retourné comme nouvelle
configuration. Sinon on la rejette et la configuration a l’etape suivante est identique
à la configuration précédente.
7. On calcule un certain nombre de variables :
P
– énergie d’interaction E = −J
P <α,β> sα sβ ,
– Aimantation par spin m = α sα /N
– etc...
8. On réitère

6.3 Analyse des résultats

Fixons la temperature à kT /J = 3 et partons d’une configuration initiale entièrement

désordonnée. On obtient, en fonction du nombre de pas Monte-Carlo l’energie E représentée
sur la figure 6.1. On peut estimer que l’équilibre est atteint lorsque l’energie ne dérive plus,
c’est à dire environ 20*N itérations (i.e après avoir tenté en moyenne 20 fois de retourner
chaque spin)
L’image instantanée du réseau au bout de 100*N itérations est représentée sur la figure 6.2.
Les spins ne sont pas répartis au hasard, on voit se développer des domaines s = +1 et des
domaines s = −1 dont la taille va croitre si on diminue la température.
L’évolution de l’aimantation par spin est représentée sur la figure 6.3.
A partir du moment où l’équilibre est atteint, on va pouvoir calculer les valeurs moyennes
des grandeurs thermodynamiques intéressantes :
– L’energie moyenne < E >
– L’aimantation moyenne par site < m >
6.3 Analyse des résultats 49

-0.5

Energie -0.6

-0.7

-0.8

-0.9
0 10 20 30 40 50 60 70 80 90 100
(nb de pas Monte-Carlo)/N

Fig. 6.1 – Evolution de l’énergie “instantanée” en fonction du nombre de pas Monte-Carlo

pour kT /J = 3. Une estimation grossière du nombre de pas nécessaires pour atteindre
l’équilibre est donné par le nombre d’iterations au bout duquel l’energie ne dérive plus,
c’est à dire environ 20*N pas (i.e on a tenté de retourner chaque spin en moyenne 20 fois)

Fig. 6.2 – Image instantanée du réseau au bout de 100*N itérations.

50 Le modèle d’Ising

0.5
aimantation par spin

-0.5

-1
0 20 40 60 80 100
(nb de pas Monte-Carlo)/N

Fig. 6.3 – Evolution de l’aimatation m pour k B T /J = 3.

– La chaleur spécifique à volume constant :

" #
∂ Ē ∂ X
Cv = = Ei exp(−Ei /kT )/Z
∂T ∂T
i

Cv = (< E 2 > − < E >2 )/kT 2

Il nous faut maintenant chiffrer l’erreur sur le calcul de ces valeurs moyennes en fonction
du nombre d’itérations.
Considérons, par exemple l’aimantation moyenne par site :
Ptmax
teq m(t)
< m >=
tmax − teq

(Chaque itération est repérée par la variable discrète t, t eq correspond au nombre d’itérations
au bout desquelles on estime avoir atteint l’équilibre et t max correspond au nombre total
d’itérations).
On peut estimer un intervalle de confiance pour cette valeur moyenne en calculant l’écart
quadratique moyen :
sP
tmax
t=teq (m(t)− < m >)2 p
σ= = < m2 > − < m > 2
tmax − teq

Mais le théorème de la limite centrale suppose que les (t max − teq ) mesures de m sont
indépendantes, ce qui n’est pas vrai. Ayant fait une mesure, il nous faut estimer le “temps”
6.3 Analyse des résultats 51

autocorrelation

0.1
0 5 10 15 20 25
(nb de pas Monte-Carlo)/N

Fig. 6.4 – Autocorrélation de l’aimantation χ(t)/χ(0)

τ (i.e. nombre d’itérations Monte-Carlo) qu’il faut attendre pour obtenir une nouvelle
mesure indépendante de la première.
Ceci peut être chiffré en calculant une “fonction d’autocorrélation”. Pour l’aimantation par
site, cette fonction est :
Z
χ(t) = dt0 [m(t0 )− < m >][m(t0 + t)− < m >]

soit Z
χ(t) = dt0 [m(t)m(t0 + t)− < m >2 ]

Pour une variable t discrète, on calculera explicitement :

tmax
X−t
1
χ(t) = m(t0 )m(t0 + t)
tmax − t
t0 =0
tmax
X−t tmax
X−t
1 1
− m(t0 ) × m(t0 + t) (6.2)
tmax − t tmax − t
t0 =0 t0 =0

qui correspond à une discrétisation de l’intégrale précédente. Dans les mêmes conditions
que pour les figures précédentes (kT /J = 3), on obtient le χ(t)/chi(0), représenté sur la
figure 6.4, en échelle semi logarithmique.
52 Le modèle d’Ising

Cette fonction suit, au départ, une loi exponentielle décroissante :

χ(t) = χ(0) exp(−t/τ )

τ représente le temps de corrélation (ici τ ≈ 15 × N )

Pour estimer la valeur moyenne de m, on pourrait n’additionner que des valeurs sélectionnées
à des intervalles de temps supérieurs à τ . On préfère additionner les valeurs obtenues ‘a
chaque temps t, mais alors le nombre de mesures indépendantes ne sera pas N = t max − teq
mais N 0 = (tmax −teq )/tau et l’intervalle de confiance pour la mesure de < m > sera obtenu
à partir de :
τ
r
σ
(tmax − teq )

6.4 Evolution en fonction de la température

On part de T = 0 et d’un état complètement ordonné. On augmente la température par pas

discrets et pour chaque température on applique l’algorithme précédent sur 200 × N pas
en effectuant les valeurs moyennes sur les 100 × N derniers. Lorqu’on varie la température,
on prend pour état initial, l’état final obtenu à la température précédente. On obtient alors

2
aimantation <m> et chaleur specifique Cv

1.5

0.5

0
0 0.5 1 1.5 2 2.5 3
Temperature: kT/J

Fig. 6.5 – Aimantation moyenne par site et chaleur specifique par site en fonction de la
température
6.4 Evolution en fonction de la température 53

pour l’aimantation moyenne et pour la chaleur spécifique les courbes représentées sur la
figure 6.5.
On observe une transition de phase vers un état complètement ordonné à kT /J ≈ 2.3. On
constate que l’erreur, en particulier sur la chaleur specifique est d’autant plus grande que
l’on s’approche de la transition. La raison est que lorsqu’on s’approche de la transitions, les
“longueurs de corrélation” qui se traduisent par la taille des domaines ordonnés à courte
distance divergent. Notre algorithme “local” n’est plus adéquat.
Ce phénomène est appelé “ralentissement critique”
54 Le modèle d’Ising
 Utilisation de biais.
7 Algorithme de Wolff

7.1 Introduction de biais

Dans l’algorithme utilisé précédemment pour la modèle d’Ising, la probabilité de choisir un

spin parmi N pour le retourner est toujours 1/N, quelle que soit la configuration de spins.
On peut envisager des algorithmes plus généraux, pour lesquels un certain mouvement de
la configuration x vers la configuration y est choisi avec une probabilité P gen (x → y) puis,
il est accepté avec une probabilité P acc (x → y). Dans ce cas, la “loi de passage” de x à y
est le produit de ces deux probabilités :
P (x → y) = Pgen (x → y)Pacc (x → y)
En travaillant dans l’ensemble canonique, la relation du bilan détaillé s’écrit :
Pgen (x → y)Pacc (x → y) exp(−E(x)/kT ) =
Pgen (y → x)Pacc (y → x) exp(−E(y)/kT ) (7.1)
Si on utilise l’algorithme de Metropolis, la règle d’acceptation sera :
 
Pgen (y → x)
Pacc (x → y) = M in 1, exp[−(E(y) − E(x))/kT ]
Pgen (x → y)
On retrouve la règle usuelle de Métropolis dans le cas où P gen est symétrique :
Pgen (x → y) = Pgen (y → x)
Par contre rien n’empêche de créer un dissymétrie telle que la probabilité de générer x à
partir de y soit supérieure à la probabilité de générer y à partir de x, d’où le nom de “biais”.
L’idéal serait de réaliser :
Pgen (x → y) exp(−E(x)/kT ) = Pgen (y → x) exp(−E(y)/kT )
de manière à avoir une probabilité d’acceptation ègale à 1. Ceci n’est possible que dans
certains cas particulier : nous allons en étudier un exemple, l’agorithme de Wolff pour le
modèle d’Ising. De manière plus générale on cherchera a augmenter la probabilité d’accep-
tation sans toujours pouvoir atteindre la limite idéale : 1.
56 Utilisation de biais. Algorithme de Wolff

7.2 Exemple : algorithme de Wolff pour le modèle d’Ising

On envisage maintenant, dans le cadre de la simulation du modèle d’Ising, le retournement

global d’une petit ensemble de spins voisins de même orientation

Définition de l’amas

1. On choisi un des N sites toujours suivant une loi uniforme comme premier site de
l’amas.
2. On regarde successivement ses 4 voisins. Tout voisin qui a la même orientation de
spin, est ajouté à l’amas avec une probabilité P (P est un nombre entre 0 et 1 fixé
une fois pour toutes).
3. On regarde les voisins des nouveaux spins ajoutés à l’amas à l’étape précédente et
qui n’appartiennent pas à l’amas. Tous ceux qui ont la même orientation que l’amas
sont ajoutés à l’amas avec la probabilité P
4. On réitère l’étape précédente
5. On s’arrête lorsqu’à une certaine étape on n’a rien rajouté à l’amas.

Retournement de l’amas

La configuration y est celle obtenue à partir de x en retournant en bloc tous les spins de
l’amas (cf. exemple de la Fig. 7.1).

Probabilité d’acceptation de la configuration retournée

1. On considère uniquement les spins de l’amas situés à la surface de l’amas. On compte

le nombre ν de leurs voisins n’appartenant pas à l’amas et qui pointent dans le même
sens que ceux de l’amas. Le nombre νx correspondant à la configuration x est en
général différent de celui νy correspondant à l’amas retouné de la configuration y. A
titre d’exemple compter les nombres ν x et νy correspondant à la Fig. 7.1.
Montrer que :
pchoix (x → y) = A(1 − P )νx
où le préfacteur A ne dépend que de la topologie de l’amas
Montrer que :
pchoix (y → x) = A(1 − P )νy
avec le même prefacteur A.
2. Montrer que :
pacc (x → y)
pacc (y → x)
 2J 
s’exprime simplement en fonction de P , exp kT et de (νx − νy ).
7.2 Exemple : algorithme de Wolff pour le modèle d’Ising 57

3. Montrer qu’un choix judicieux de P conduit à

pacc (x → y) = pacc (y → x) = 1

c’est à dire à l’acceptation certaine du retournement à chaque étape.

4. Pourquoi cet algorithme est-il nettement supérieur àu précédent ?

Fig. 7.1 – Choix d’un amas de la configuration de départ x (à gauche) et retournement
des spins de l’amas pour une nouvelle configuration d’essai y (à droite).
58 Utilisation de biais. Algorithme de Wolff
 Applications à la physique
8 médicale

8.1 Introduction

Ce chapitre est consacré à la simulation de processus stochastiques par essence : le transport

de radiations (photons) ou de particules massives (neutrons, électrons ...) dans la matière.
Il illustre l’une des applications les plus importantes de la méthode de Monte-Carlo. N’ou-
blions pas que cette méthode a vu son essor, à partir des dernières années de la seconde
guerre mondiale, à travers la simulation sur ordinateur du transport de neutrons dans la
matière.
Une particule, massive (neutron, électron...) ou non massive (photon X, photon γ) se
déplaçant dans la matière a une certaine probabilité de rencontrer un atome ou une molécule
avec lequel (ou laquelle) elle va interagir. Après interaction, sa trajectoire, son energie vont
être modifiées suivant certaines lois de probabilité. Entre deux interactions, la trajectoire
de la particule est rectiligne.
Une particule, dans la matière suit un chemin aléatoire constitué de segments de droite.
Les étapes (segments) successives de ce chemin constituent une chaı̂ne de Markov, car la
probabilité pour une particule, –ayant subi une interaction au point de coordonnée x au
temps t–, de subir une nouvelle interaction au point de coordonée y au temps t+1, après
un parcours rectiligne ne dépend que de la position x au temps t (indépendemment de tout
ce qui s’est passé aux temps antérieurs à t).
Nous nous intéresserons essentiellement ici au transport de photons (rayons X, rayons γ)

8.2 Modélisation d’un rayonnement

Un rayonnement X ou γ peut être décrit par un ensemble de photons. La simulation peut

être effectuée en suivant l’évolution (trajectoire, énergie) d’un photon. On pourra aussi
simuler en parallèlle l’évolution d’un ensemble de photons, et dans ce cas l’usage d’or-
dinateurs massivement parallèles s’avérera particulièrement efficace (bien qu’il soit aussi
possible de le faire sur un ordinateur monoprocesseur !)
60 Applications à la physique médicale

A un instant donné, un photon individuel est caractérisé par sa position r dans l’espace
(trois variables), sa direction d̂ (vecteur unitaire defini par deux variables indépendantes)
et son énergie E.
Entre deux collisions, le photon conserve sa direction et son énergie. Le parcours du photon
sera donc caractérisé par un ensemble discret d’états {r n , d̂n , En } représentant sa position,
sa direction son énergie juste après une interaction (“collision”) avec une particule du milieu
dans lequel il évolue.

8.3 Modélisation de la source

La première étape de la simulation est bien sûr de définir létat initial {r 0 , d̂0 , E0 } du photon.
Ces variables initiales seront sélectionnés suivant une densité de probabilité représentant la
source de rayonnement.
S(r0 , d̂0 , E0 )

8.4 Modélisation d’une trajectoire

Supposons que la trajectoire d’un photon ait été simulée jusquà l’état {r n , d̂n , En } juste
après la nieme collision. Le photon va continuer son chemin en ligne droite le long de la
direction donnée par le vecteur unitaire d̂n en parcourant une distance s puis subir une
(n + 1)ieme collision au point :
rn+1 = rn + sd̂n
Comment choisir cette distance aléatoire s ? Nous allons montrer que s est une variable
aléatoire de densité de probabilité :

1
p(s) = exp(−s/λ)
λ
Nous considérons un milieu homogène constitué de molécules de section efficace σ. Le
nombre de molécules par unité de volume est N .
La probabilité d’interaction, c’est à dire de rencontre d’une molécule lorsque le photon
parcourt un petit chemin ds est indeépendante de la position r et de l’orientation d̂ de
sa trajectoire. Elle est simplement proportionnelle à la longueur du chemin élémentaire
parcouru ds soit :
ds
λ
où λ est une constante —que nous préciserons par la suite— qui ne dépend que du milieu
homogène considéré.
Soit p(s) la densité de probabilité correspondant à la loi de probabilité de la longueur s du
chemin parcouru entre deux collisions.
8.4 Modélisation d’une trajectoire 61

La probabilité de parcourir un chemin de longueur s sans interagir est

Z s Z ∞
F(s) = 1 − p(u)du = p(u)du (8.1)
0 s

La probabilité de parcourir un chemin de longueur s sans interagir, puis d’interagir entre

s et s + ds est le produit des deux probabilités précédentes :

ds
p(s)ds = F(s)
λ
soit ∞
ds
Z
p(s)ds = p(u)du
λ s

p(s) obéit donc à l’équation intégrale :

∞
1
Z
p(s) = p(u)du
λ s

En dérivant chaque membre de cette équation par rapport à s, on obtient :

dp(s) 1
= − p(s)
ds λ
p(s) est donc exponentielle : p(s) = A exp(−s/λ). La constante A est déterminée par la
normalisation : Z ∞
p(s)ds = 1
0
Ce qui impose :
1
p(s) = exp(−s/λ)
λ
Il est immédiat de vérifier que l’on a :
Z ∞
s p(s)ds = λ
0

λ correspond à l’espérance ou valeur moyenne de la variable aléatoire associée au chemin

parcouru entre deux collisions. λ correspond donc au libre parcours moyen du photon. En
fonction de la densité N de molécules et de leur section efficace, il s’exprime par la relation :

λ = 1/(N σ)

La question est maintenant d’échantillonner des valeurs de s suivant la densité de probabilité

exponentielle p(s), ce qui est facile par la méthode de transformation de variables. La
fonction de partition correspondant à p(s) est :

1 s
Z
η = F (s) = exp(−u/λ)du = 1 − exp(−s/λ)
λ 0
62 Applications à la physique médicale

On inverse cette relation :

s = −λ ln(1 − η)

Donc concrètement on tire des nombres η i suivant la loi uniforme sur le segment [0,1], et
les nombres
si = −λ ln(1 − ηi )

sont alors répartis suivant la densité de probabilité exponentielle souhaitée. Remarque : il

est équivalent de choisir :
si = −λ ln(ηi )

puisque si η est réparti uniformément sur le segment [0,1], alors 1 − η l’est aussi.

8.5 Interaction photon-matière

Il nous reste à échantillonner le nouveau vecteur directeur d̂n+1 après l’interaction suivante
et la nouvelle énergie. Ceci va dépendre du type d’interaction.

8.5.1 Diffusion élastique de Rayleigh

Pour une diffusion élastique, il y a changement de la direction de propagation, mais con-

servation de l’énergie. La nouvelle direction de propagation d̂n+1 est définie par rapport
à la précédente d̂n+1 par les deux angles φ et θ indiqués sur la figure précédente. L’angle
azimuthal φ obéit à une loi de probabilité de densité uniforme sur le segment [0, 2π]. Par
contre, la loi de probabilité de l’angle θ n’est pas uniforme.
Pour un photon d’énergie E subissant une diffusion élastique sur un atome de numéro
atomique Z, la probabilité pour le photon d’être diffusé dans un petit angle solide dΩ

dn+1
φ
dn θ

Fig. 8.1 –
8.5 Interaction photon-matière 63

autour de la direction donnée par {θ, φ} est, dans le cadre de l’approximation de Born :
1 + cos2 θ 2
dσRayleigh = re2 F (q, Z)dΩ
2
q = |~q| est le module du transfert de moment :
q~ = ~kn+1 − ~kn
où ~kn et ~kn+1 sont respectivement les vecteurs d’onde incident et diffusés de directions d̂n
et d̂n+1 et de même module :
|~kn | = |~kn+1 | = E/c
(c est la vitesse de la lumière). On a donc, en fonction de E et θ :
q = 2(E/c) sin(θ/2)
re est le rayon de Bohr de l’électron et F (q, Z) est le facteur de forme de l’atome, c’est à
dire la transformée de fourier de sa densité électronique ρ(~r). Pour un atome de symétrie
sphérique, on a : Z ∞
sin(qr/~) 2
F (q, Z) = 4π ρ(r) r dr
0 qr/~
Pour échantillonner θ suivant la densité de probabilité correspondant à la “section efficace
de diffusion” dσRayleigh /dΩ exprimée ci dessus, on aura en général recours à la méthode de
réjection de Von Neumann.

8.5.2 Diffusion inélastique de Compton

Dans une diffusion de Compton, un photon d’énergie E n interagit avec un électron d’un
atome. Le photon transfère une partie de son énergie à l’électron qui est alors éjecté de
l’atome. Le photon est alors diffusé suivant les angles {θ, φ} par rapport à la direction
incidente d̂n avec une énergie En+1 < En . L’angle azimutal φ est, ici aussi, distribué
uniformément sur le segment [0, 2π].
Nous nous restreindrons ici au cas où l’electron éjecté est, au départ, au repos. La conserva-
tion de l’énergie et de l’impulsion du système total {photon+électron} impose une relation
entre l’angle θ et les energies :
me c2 me c2
− = 1 − cos θ
En+1 En
Les quantités me c2 /E = λCompton représentent les “longueurs d’onde de Compton” ( m e
est la masse de l’électron au repos).
La section efficace de diffusion, c’est à dire la densité de probabilité correspondant à l’angle
θ est donnée par la relation de Klein-Nishina :
re2 En+1 2 En+1
   
dσCompton En 2
= + − sin θ
dΩ 2 En En En+1
64 Applications à la physique médicale

8.5.3 Absorption photoélectrique

L’absorption photoélectrique est un processus où le photon est absorbé par un électron qui
effectue une transition vers un état d’énergie plus élevée. La trajectoire du photon prend
donc fin à cet endroit.

8.5.4 Choix d’un des événements

En fonction du milieu considéré, de l’énergie des photons etc... les probabilités de chacun des
trois événements décrits ci-dessus : P Rayleigh , PCompton , Pabsorption sont connues (PRayleigh +
PCompton + Pabsorption = 1). A chaque étape n il faudra tirer un de ces trois événement
suivant ces probabilité.

8.5.5 Pour en savoir plus...

Il existe des logiciels programmés pour simuler toutes les interactions possibles de partic-
ules non massives (photons) ou massives (neutrons, electrons, positrons), chargées ou non
chargées dans la matière. L’un d’entre eux, développé par un groupe de scientifique est
nommé “PENELOPE”.Le code source est payant mais les méthodes utilisées sont en libre
accès sur internet.
http ://www.nea.fr/html/dbprog/penelope-2003.pdf

8.6 Propagation dans un milieu inhomogène

En physique médicale, il est rare d’avoir un milieu homogène. Dans la plupart des cas,
un photon traverse une succession de milieux de densités et sections efficaces différentes.
Exemple : pour la simple radiographie d’une racine dentaire, les rayon X vont traverser la
gencive, puis le tissus osseux, puis la pulpe etc...
Comme indiqué sur la figure 8.2, considérons un photon, qui après avoir subi une interaction
en rt traverse successivement différents tissus T 1 , T2 , ... , Tn−1 avant de subir l’interaction
suivante en rt+1 dans le tissus Tn . Nous noterons respectivement λ1 , λ1 , ... , λn les libres
parcours moyens du photon dans les tissus T 1 , T2 , ... , Tn . Comme l’indique la figure, le
photon a parcouru une distance s1 dans le tissus T1 , une distance s2 dans le tissus T2 ...
une distance sn−1 dans le tissus Tn−1 et enfin une distance (s − sn−1 − sn−2 − · · · − s2 − s1)
dans le dernier tissus Tn .
D’après ce qui précède, la probabilité de parcourir une distance s i dans la tissus Ti sans
interagir est (cf. Equation 8.1)
∞    
1 u si
Z
F(si ) = exp − du = exp −
λi si λi λi
8.6 Propagation dans un milieu inhomogène 65

λ2
λ4
s2
s1
s4
s3

λ3
λ1

λ5

Fig. 8.2 –

Et la probabilité d’avoir interagi entre s et s + ds dans le milieu T n est :

ds
λn
La probabilité d’avoir parcouru le chemin de longueur s à travers les divers tissus sans
interagir, puis d’avoir interagit entre s et s + ds est donc :

[s − s1 − s2 − · · · − sn−1 ] ds
       
s1 s2 sn−1
p(s)ds = exp − exp − · · · exp − exp −
λ1 λ2 λn−1 λn λn

En definissant, par simple translation la variable :

λn λn λn
s̃ = s1 + s2 + · · · + sn−1 + s − s1 − s2 − · · · − sn−1
λ1 λ2 λn−1
on a  
1 s̃
p(s̃) = exp −
λn λn
On peut donc générer s̃, comme précédemment suivant la loi exponentielle :
– On tire des nombres η suivant une loi uniforme sur le segment [0,1]
66 Applications à la physique médicale

– les nombres
s̃ = −λn ln η
sont générés suivant la loi exponentielle ci dessus
– on en déduit :
 
s1 s2 sn−1
s = s1 + s2 + · · · + sn−1 − λn ln η + + +··· + (8.2)
λ1 λ2 λn−1
On procédera donc, en général de la manière suivante. La direction du photon étant donnée,
on connnait donc ses intersections avec les différentes frontières du milieu et les longueurs
s1 , s2 , ...
– On tire un nombre η uniformément distribué sur le segment [0,1]
– si − ln η < λs11 alors l’interaction suivante a lieu dans le milieu T 1 , il n’y a rien de nouveau
par rapport aux paragraphes précédents.
– si λs11 < − ln η < λs11 + λs22 alors l’interaction suivante a lieu dans le milieu T 2 et on prend :
 
s1
s = s 1 − λ2 ln η +
λ1
– etc ...
– si
s1 s2 sn−1 s1 s2 sn−1 sn
+ + ··· + < − ln η < + + ··· + +
λ1 λ2 λn−1 λ1 λ2 λn−1 λn
alors l’interaction a lieu dans le milieu T n et s est donné par la relation 8.2
 Méthode de Monte Carlo et
9 Dynamique Moléculaire

9.1 Introduction

Bien que la méthode de “Dynamique Moléculaire” soit parfaitement déterministe, donc

en dehors du cadre de ce cours, nous en donnons ici les principes, car dans les simula-
tions numériques elle est souvent employée pour les mêmes problèmes, conjointement à
la méthode de Monte-Carlo. Dans ce chapitre, nous soulignerons les avantages et les in-
convénients de chacune de ces méthodes.
La “dynamique moléculaire” simule ab initio l’évolution d’un système physique –en parti-
culier son évolution vers l’équilibre–, à partir des équations qui régissent ce système. Pour
un système classique, ces équations sont simplement celles de la mécanique newtonienne.
En préambule, et pour rendre la comparaison entre les deux méthodes plus évidente, nous
allons montrer qu’une variante de la Méthode de Monte-Carlo, appliquée à la simulation de
la situation déquilibre d’un système classique, peut se construire autour de la simulation
d’une équation différentielle stochastique.

9.2 Equation stochastique

Nous allons montrer qu’une chaı̂ne de Markov, satisfaisant le principe du bilan détaillé peut
être obtenue à partir de l’intégration d’une équation différentielle stochastique (i.e. faisant
intervenir des variables aléatoire).

Nous reprenons notre système S caractérisé par les coordonnées : x = (α 1 , α2 , · · · , α3N )

de ses N particules. Son energie est E(x), et nous introduisons “l’action” S(x) = βE(x).
Nous considérons l’équation de Langevin :
∂αi ∂S(x)
= −Γi + ξi (τ ) (9.1)
∂τ ∂αi
où les αi sont des fonctions d’une variable continue τ . Γ i est un paramètre d’echelle arbi-
traire et les ξi (τ ) représentent un ensemble de variables aléatoires indépendantes gaussi-
68 Méthode de Monte Carlo et Dynamique Moléculaire

ennes centrées de même variance 2Γ i

E[ξi (τ )ξj (τ 0 )] = 2Γi δi,j δτ −τ 0 (9.2)

Cette équation représente “un bruit blanc”. Remplaçons l’équation de Langevin continue

par une équation aux différences finies :

 
∂S(x)
αi (τn+1 ) = αi (τn ) + ∆τ −Γi + ξi (τn ) (9.3)
∂αi (τn )

de sorte que αi (τn ) correspond au nieme élément d’une chaı̂ne de Markov.

La distribution de probabilité des variables aléatoires ξ i correspondant à l’équation 9.2 est :

r
Y ∆τ − Pi Rττnn+1 4Γdt
ξi (t)2
e i
4πΓi
i

soit r
Y ∆τ − Pi 4Γ
∆τ
ξi (τn )2
e i (9.4)
4πΓi
i

La règle p(x → y) pour cette chaı̂ne de Markov est que la probabilité de passer de x =
(α1 , α2 , · · · , α3N ) à y = (β1 , β2 , · · · , β3N ) est égale à la probabilité pour la variable aléatoire
(bruit blanc) ξi dêtre égale à
βi − α i ∂S(x)
+ Γi
∆τ ∂αi

Il nous reste à vérifier que cette règle obéit à l’équation du bilan détaillé. On a
P “ βi −αi ”2
∆τ ∂S
− 4Γ +Γi ∂α
p(x → y) e i i ∆τ i
= ”2
p(y → x)
“
∆τ P αi −βi ∂S
− 4Γ i +Γi ∂β
e i ∆τ i

soit lorsque ∆τ → 0
p(x → y) P ∂S
− i (βi −αi ) ∂α
=e i
p(y → x)
ce qui peut encore s’écrire en identifiant S et son développement au premier ordre :

p(x → y) e−S(y)
→ −S(x) lorsque ∆τ → 0
p(y → x) e

CQFD
9.3 Comparaison avec la Dynamique Moléculaire 69

9.3 Comparaison avec la Dynamique Moléculaire

Nous avons jusqu’ici toujours considéré l’ensemble canonique : notre système S est en
contact avec un réservoir ou thermostat R. La température T du système est fixée par
le thermostat, mais son energie fluctue autour d’une valeur moyenne < E >. La valeur
moyenne d’une observable O(x) est :

O(x)e−βH(p,x) dpdx
R
< O >= (9.5)
Z
où Z représente la fonction de partition :
Z
Z = e−βH(p,x) dpdx (9.6)

x = (α1 , α2 , · · · , α3N ) représente les coordonnées des particules et p = (p 1 , p2 , · · · , p3N )

les variables conjuguées ou impulsions : p i /mi = dαi /dt.
X p2
i
H= + E(x) (9.7)
2mi
i

représente l’Hamiltonien du système. Le premier terme représente l’énergie cinétique et le

second terme E(x) l’energie potentielle.

Examinons la formulation de la méthode de dynamique moléculaire dans l’ensemble micro-

canonique.

On considère un système isolé, ayant donc une énergie E fixée. Les trajectoires des particules
sont déterminées par les équations de Hamilton :

dαi /dt = pi /mi

dpi /dt = −∂E(x)/αi (9.8)

Ces équation de mouvement sont intégrées numériquement. En suivant l’évolution d’une

observable O(x(t)) pendant un temps très long, on peut calculer sa valeur moyenne dans
le temps :
1 T
Z
< O >t = limT →∞ O(x(t))dt (9.9)
T 0

Si le système est ergodique, alors on montre que, dans la limite thermodynamique, la

moyenne d’ensemble < O > [relation 9.5] calculée dans l’ensemble canonique est égale
à la moyenne dans le temps < O >t calculée dans l’ensemble microcanonique, pour une
valeur fixée E de l’énergie égale à la valeur moyenne < E > de l’energie dans l’ensemble
canonique.
70 Méthode de Monte Carlo et Dynamique Moléculaire

Dans la limite thermodynamique, on passe de l’ensemble canonique à l’ensemble micro-

canonique en remarquant que, dans l’ensemble canonique, l’energie fluctue de moins en
moins autour de sa valeur moyenne < E >. On peut alors appliquer l’approximation du
col à l’intégrale représentée par la relation 9.5. Dans ce cas, la contribution essentielle à
l’intégrale provient de l’hypersurface d’energie constante définie par :

H(p, x) =< E >

En général, les méthodes de dynamique moléculaire procèdent à l’intégration des équation

du mouvement par discrétisation du temps. On suit l’évolution du système à des temps
τk = k∆τ , espacés régulièrement de ∆τ . Les équations du mouvement :

d2 αi ∂E(x)
mi =− (9.10)
dt2 ∂αi
conduisent aux équations d’évolution en temps discret :
1 ∂E(x)
αi (τn+1 ) = 2αi (τn ) − αi (τn−1 ) − (∆τ )2 (9.11)
mi ∂αi
Ces équations constituent le fondement de l’algorithme de Verlet[10] le plus courammnent
utilisé en dynamique moléculaire[11]. On peut remanier les termes de cette équation pour
les écrire de manière plus suggestive :
1 ∂E 1
αi (τn+1 ) = αi (τn ) − (∆τ )2 + [αi (τn+1 ) − αi (τn−1 )] (9.12)
2mi ∂αi 2

Il intéressant de comparer cette équation à celle obtenue après discrétisation de l’équation

de Langevin pour un échantillonnage suivant la méthode de Monte-Carlo. Dans l’équations
9.3, le paramètre Γi est arbitraire. Prenons
∆τ
Γi =
2βmi
alors, on obtient :
1 ∂E
αi (τn+1 ) = αi (τn ) − (∆τ )2 + ∆τ ξi (τn ) (9.13)
2mi ∂αi

En comparant les équations 9.12 et 9.13, on constate que les vitesses des particules :
1
vi = [αi (τn+1 ) − αi (τn−1 )]
2∆τ
qui correspondent au dernier terme de l’équation d’évolution 9.13 en dynamique moléculaire
sont remplacées, dans la méthode de Monte-Carlo par des variables aléatoires gaussiennes
indépendantes ξi (bruit stochastique)[8, 12].
9.3 Comparaison avec la Dynamique Moléculaire 71

La méthode de dynamique moléculaire a l’avantage de conduire directement dans les régions

de l’espace des phases {p, x} qui contribuent le plus à l’intégrale donnant la valeur moyenne
d’une observable, mais dans certains cas non triviaux où un temps très long est nécessaire
pour parcourir l’espace des phases, on peut “manquer” certaines régions importantes de cet
espace. La méthode de Monte-Carlo explore tout l’espace des phases, mais on peut perdre
du temps à stagner dans une région de l’espace des phases qui contribue peu à l’intégrale.

La comparaison ci-dessus des deux méthode ouvre la possibilité d’algorithmes hybrides[12,

13] où l’évolution du système aux temps discrets successifs correspond à une chaı̂ne de
Markov déterminée par l’equation :
1 ∂E
αi (τn+1 ) = αi (τn ) − (∆τ )2 + ∆τ ζi (τn ) (9.14)
2mi ∂αi

ζi est soit la vitesse (deterministe !) des particules v i , soit un bruit stochastique ξi . La

manière d’alterner dans la chaı̂ne de Markov entre ces deux possibilités est propre à chaque
problème. Cette méthode est à la fois plus rapide que la méthode de Monte-Carlo conven-
tionnelle et plus sûre que la dynamique moléculaire pure.
72 Méthode de Monte Carlo et Dynamique Moléculaire
 10 Optimisation

10.1 Méthode du “Recuit Simulé” (Simulated Annealing)

Lorsqu’on souhaite minimiser une fonction à grand nombre de variables, les méthodes
classiques n’ont aucune difficulté pour trouver des minima locaux (en général le plus proche
d’un point défini comme point de départ de la méthode). Il est par contre pratiquement
impossible de trouver le minimum absolu. Seule, la méthode de Monte-Carlo permet de
résoudre ce problème.
La fonction E(x1 , x2 , .., xN ) à minimiser est considérée comme l’énergie d’un système fictif
ayant pour espace de configuration R N . Ce système est supposé en interaction avec un
thermostat à température T . On simule l’evolution du système dans l’ensemble Canonique
par la Méthode de Monte-Carlo.
Comme en métallurgie, le “recuit” consiste a “chauffer le système” à une température
élevée puis à diminuer très doucement la température de manière à obtenir l’équilibre
thermodynamique à chaque température. Lorsque la température s’approche de zéro, le
systéme atteint son énergie minimale, c’est à dire le minimum absolu de la fonction E.

10.2 “Recuit Parallèle” (Parallel Tempering)

Lorsque la fonction E(x1 , x2 , .., xN ) est particulièrement compliquée avec un très grand
nombre de minima relatifs, la méthode simple précédente ne permet pas forcément à basse
température d’atteindre l’équilibre en un temps raisonnable, donc de converger vers le
minimum absolu à température nulle. On peut facilement resté “piégé” à basse température
dans une région de l’espace multidimensionnel, dont on ne s’échappe plus.
Dans ce cas, une méthode dite “Parallel Tempering” en anglais, —que nous traduirons
par “Recuit parallèle”— a été proposeée, il y a une vingtaine d’années par les physiciens
des “verres de spin”. Elle s’applique aussi bien à l’obtention de l’équilibre d’un système
thermodynamique comprenant un très grand nombre d’états d’énergies voisines —c’est le
cas des solides désordonnés ou “verres”, ou d’assemblée de spins avec fort désordre dans les
interactions entre spins (“verres de spins”)— qu’à un problème de minimisation de fonction
74 Optimisation

présentant les mêmes caractéristiques.

Il s’agit d’effectuer parallèlement plusieurs simulations du même système à des températures
différentes, avec des échanges de configurations, selon des règles que nous allons définir.
Considérons le cas (schématisé à une dimension sur la figure 10.1 où deux régions de l’espace
A et B sont séparées par une importante barrière d’énergie ∆E. A haute température
kTH >> ∆E, il est évident que l’algorithme de Métropolis permet d’explorer les deux
régions A et B. Par contre, à basse température kT B << ∆E, on va resté piégé dans une
des deux régions A ou B.
L’idée est de simuler en parallèle le même système à deux températures différentes T H et
TB et d’échanger “de temps en temps” les configurations des deux systèmes, avec une prob-
abilité que nous allons définir. Cet “échange” va permettre au système à basse température
de visiter les deux régions A et B.
La condition d’ergodocité (visiter tout l’espace) sera ainsi réalisée. Il nous reste à définir la
probabilité d’échange des configurations de manière à satisfaire la “microréversibilité”.
Soit EH L’énergie de la configuration du système à haute température T H et EB celle
du système à basse température T B . Nous allons montrer que si on accepte l’échange des
configurations avec la probabilité :

exp[−(βB − βH )(EH − EB )] si (EH − EB ) > 0 (10.1)

1 sinon

alors la condition de microréversibilité est satisfaite (on note β H = 1/kTH et βB = 1/kTB )

Considérons la probabilité conjointe p µ,ν pour que le système à basse température soit
dans l’état µ et que le système à haute température soit dans l’état ν. A l’équilibre les deux
systèmes doivent satisfaire la distribution canonique de Boltzmann :
1
pµ,ν = exp[−βB E(µ)] exp[−βH E(ν)] (10.2)
ZB ZH
Toutes les M itérations, nous tentons un échange des configurations µ et ν avec la probabilité
P ({µ, ν} → {ν, µ}) En appliquant la loi 10.1, on a :

P ({µ, ν} → {ν, µ})

= exp[−(βB − βH )(E(ν) − E(µ))]
P ({ν, µ} → {µ, ν})

Mais la loi de Boltzmann 10.2 conduit à :

pµ,ν exp[−βB E(µ)] exp[−βH E(ν)]
= = exp[−(βB − βH )(E(µ) − E(ν))]
pν,µ exp[−βB E(ν)] exp[−βH E(µ)]

On a donc :
P ({µ, ν} → {ν, µ}) pν,µ
=
P ({ν, µ} → {µ, ν}) pµ,ν
et la condition de microréversibilité est bien vérifiée.
10.2 “Recuit Parallèle” (Parallel Tempering) 75

kTH
E(X)

0.5

A B
kTB

0
0 0.5 1
X

Fig. 10.1 –

En ce qui concerne les M-1 autres itérations, nous tentons d’abord une transition µ → µ 0 du
système à basse température, suivant l’algorithme habituel de Métropolis, et il est facile de
vérifier à partir de 10.2 que la condition de microréversibilité est vérifiée. Puis nous tentons
une transition ν → ν 0 du système à haute température, et il en est de même.
Comment choisir M ?
Entre deux tentatives déchange, il convient que le système à haute température soit passé
à un état décorrélé du précédent. M sera donc de l’ordre du temps d’autocorrélation de
l’énergie.
Généralisation à plus de deux systèmes
Le schéma que nous venons d’établir pour deux systèmes se généralise de manière triv-
iale à n systèmes placés à n températures différentes. Il se prête particulièrement bien à
l’usage d’ordinateurs massivement parallèles, d’autant plus qu’il est absolument équivalent
d’échanger les températures plutôt que les états, et ceci sera beaucoup moins coûteux en
temps de communication entre processeurs. On peut découper un intervalle de températures
en segments relativement petits, avec une probabilité importante d’échanger les états de
deux systèmes de températures voisines.
76 Optimisation
 11 Polymères

11.1 Introduction

Dés les années 50, les pionniers de la méthode Monte-Carlo, en particulier les époux Rosen-
bluth, on simulé des systèmes de polymères.
Nous allons étudier quelques algorithmes spécifiques pour la simulation de longues chaı̂nes
de molécules.

11.2 Modélisation de polymères

11.2.1 Modèle du “collier de perles”

C’est le modéle le plus simple. Chaque monomère est modélisé par une sphère dure de
rayon a. Chaque sphère dure est accolée à ses deux voisines. Les mouvements de l’une par
rapport à l’autre de ces sphères sont libres.

Fig. 11.1 – Modélisation suivant le “collier de perles”

78 Polymères

11.2.2 Modèle du “Chapelet de ressorts” (String beads)

1. La cohésion de la chaı̂nes est modélisée par un potentiel attractif de “ressort” entre

deux molećules consécutives :

Eel = K|Ri+1 − Ri |2 (11.1)

2. On prend en compte un potentiel entre paires de monomères, fortement répulsif à

courte distance a, —de l’ordre de grandeur du diamètre de sphère dures du modèle
plus simple du collier de perles—, et faiblement attractif à longue distance (pour
tenir compte des forces de Van-der-Waals). On prendra par exemple le potentiel de
“Lennard-Jones” :
 
a a
Vi,j =  − (11.2)
|Ri − Rj |12 |Ri − Rj |6

R1
R4

R2
R3 R5

R6

R7

Fig. 11.2 – Modélisation suivant le “chapelet de ressorts”

11.3 Mouvements dans un algorithme de Metropolis 79

11.3 Mouvements dans un algorithme de Metropolis

11.3.1 Reptation

On supprime le premier monomère d’une chaı̂ne puis on le reconstruit au bout de la chaı̂ne

11.3.2 Mouvements de pivot

1. On choisit, au hasard un monomère P de la chaı̂ne
2. On choisit suivant une répartition uniforme un vecteur unitaire û dans un angle solide
centré autour du segment defini par le monomère P et le suivant (P + 1)
3. On choisit un angle φ suivant une loi uniforme sur le segment [0, 2π]
4. On effectue une rotation en bloc d’un angle φ autour de û de la partie du polymère
comprise entre le monomère (P + 1) et le dernier monomère.

11.4 “Reconstruction biaisée”

La méthode du “pivot” que nous avons décrite au paragraphe précédent est efficace pour
un ensemble de polymères “dilués i.e. pas trop proches les uns des autres. Pour un système

Fig. 11.3 – Mouvement de Pivot

80 Polymères

dense, lorsqu’on tente de bouger un grand nombre de monomères par la méthode du pivot,
par exemple, on a de grandes chances que dans la nouvelle configuration d’essai certains
“coeurs durs” se recouvrent, auquel cas le mouvement est rejeté.
En 1954, les époux Rosenbluth[14] ont proposé une technique astucieuse consistant à
supprimer une chaı̂ne de polymr̀es pour la reconstruire monomère par monomère à un
autre endroit en utilisant un “biais” pour éviter le recouvrement de couers durs. Cette
méthode a été perfectionée[15] et est maintenant couramment utilisée sous le nom anglais
de “configurational-bias Monte-Carlo algorithm”

1. Dans la configuration initiale x, on choisit une chaı̂ne Γ a comprenant N monomères

Γai
2. On supprime cette chaı̂ne pour la reconstruire ailleurs, pas à pas de la manière suiv-
ante :
– On choisit le premier monomère Γb1 de manière uniforme.
– On suppose ici, pour simplifier que les segments joignant deux monomères ne peu-
vent prendre que k orientations discrètes, ce qui est le cas si on travaille sur un
espace discrétisé, mais la méthode sétend facilement au cas continu[15].
On choisit pour le second monomère Γ b2 de la cha ı̂ne une des k orientations possibles
avec la probabilité :
exp(−ub2 /kT )
P2b =
Z2b
avec
k
X
Z2b = exp(−ub2,j /kT )
j=1

ub2,j représentant l’energie d’interaction du monomère Γ b2 avec toutes les autres

particules du système.
– On construit le second monomère de la même manière, etc.. jusqu’au N ieme monomère.
– A la fin du processus de construction, la conformation finale Γ b du polymère a été
générée avec une probabilité :

b exp(−ubi /kT ) exp[−(ub1 + ub2 + ub3 + · · · + ubN )/kT ]

Pgen = ΠN
i=2 =
Zib Wb

où
W b = exp(−ub1 /kT )ΠN b
i=2 Zi

Mais (ub1 + ub2 + ub3 + · · · + ubN ) n’est autre que la contribution totale U b de la
“nouvelle chaı̂ne”, dans sa nouvelle conformation Γ b à l’energie totale du système.
Et on peut c̀rire :
b exp(−U b /kT )
Pgen =
Wb
W b est appelè le “facteur de Rosenbluth”
11.4 “Reconstruction biaisée” 81

– Pour satisfaire la condition de microreversibilité, il faut calculer le facteur de Rosen-

bluth :
W a = exp(−ua1 /kT )ΠN a
i=2 Zi

de la chaı̂ne dans sa conformation initiale Γ a et on a :

a exp(−U a /kT )
Pgen =
Wa
U a représentant la contribution de la chaı̂ne dans son ancienne conformation Γ a à
l’energie totale du système.
3. On accepte la nouvelle conformation Γ b avec la probabilité :
 a 
Pgen
M in 1, exp[−(U b − U a )/kT ] b = M in[1, W b /W a ]
Pgen
82 Polymères
 12 Monte-Carlo quantique

12.1 Monte-Carlo Variationnel

Pour calculer l’energie fondamentale E 0 (à température nulle) d’un système quantique
décrit par un Hamiltonien H, on utilise couramment l’approximation variationnelle. Cette
approximation consiste à :
– Partir d’une fonction d’onde |ψλ,µ,··· ,ν i censée représenter un bonne approximation de
l’état fondamental du système (on n’exigera pas forcément que la norme de cette fonction
soit 1), et dépendant d’un certain nombre de param t̀res ajustables λ, µ, · · · , ν
– Calculer :
hψλ,µ,··· ,ν |H|ψλ,µ,··· ,ν i
E=
hψλ,µ,··· ,ν |ψλ,µ,··· ,ν i

On montre (théorème de Rietz ou principe variationnel) l’inégalité E 0 < E

– Le principe variationnel nous dit que la meilleure fonction d’essai est celle qui minimise
l’energie E. On ajustera donc les paramètres λ, µ, · · · , ν pour obtenir l’énergie minimale.
Que l’on parte d’une représentation de |ψ λ,µ,··· ,ν i dans l’espace réel à 3N dimensions, N
représentant le nombre de particules, ou d’uneD représentation dans E l’espace
D réciproque E
∗ ∗
(vecteurs d’ondes K), chacune des quantités ψλ,µ,··· ,ν |H|ψλ,µ,··· ,ν et ψλ,µ,··· ,ν |ψλ,µ,··· ,ν
représente une intégrale multidimensionnelle. Dans l’espace réel :

∗
R
ψλ,µ,··· ,ν (r1 , r2 , ..., rN )Hψλ,µ,··· ,ν (r1 , r2 , ..., rN )dr1 dr2 ...drN
E= R ∗
ψλ,µ,··· ,ν (r1 , r2 , ..., rN )ψλ,µ,··· ,ν (r1 , r2 , ..., rN )dr1 dr2 ...drN

D’ou l’idée de calculer directement ce rapport de deux intégrales par la méthode de Monte-
Carlo.
Il est alors naturel de choisir pour densité de probabilité :

∗
ψλ,µ,··· ,ν (r1 , r2 , ..., rN )ψλ,µ,··· ,ν (r1 , r2 , ..., rN )dr1 dr2 ...drN
p(r1 , r2 , ..., rN ) = R ∗
ψλ,µ,··· ,ν (r1 , r2 , ..., rN )ψλ,µ,··· ,ν (r1 , r2 , ..., rN )dr1 dr2 ...drN
84 Monte-Carlo quantique

On échantillonne cette densité de probabilité par la méthode de Metropolis, par exemple.

Une estimation de E est donnée par l’espérance de :
∗
ψλ,µ,··· ,ν (r1 , r2 , ..., rN )Hψλ,µ,··· ,ν (r1 , r2 , ..., rN )
∗
ψλ,µ,··· ,ν (r1 , r2 , ..., rN )ψλ,µ,··· ,ν (r1 , r2 , ..., rN )

12.2 “Diffusion Monte-Carlo”

Nous allons décrire un algorithme qui permet de trouver, sans approximation, l’énergie et
l’état fondamental de l’équation de Schrödinger :

H|Ψ >= E|Ψ >

Par souci de simplicité, nous considérons une particule de masse m, soumise à un potentiel
V (x), se déplaçant sur un axe x. La méthode est bien-sûr généralisable à un système à
plusieurs dimensions.
L’Hamiltonien s’écrit donc simplement :

~2 ∂ 2
H =− + V (x)
2m ∂x2
et l’évolution de la fonction d’onde Ψ(x, t) en fonction du temps est gouvernée par l’équation
de Schrödinger :
∂|ψ >
i~ = H|ψ >
∂t
Soit {|φ >n } une base orthonormée de fonctions propres associées aux energies E n , c’est à
dire :
H|φn >= En |φn >
avec la condition : Z ∞
< φn |φm >= φ∗n (x)φm (x)dx = δnm
−∞

Alors on peut écrire en développant |Ψ(x, t) > sur cette base :

∞
X
|Ψ(x, 0) = cn Φn (x)
n=0
avec Z ∞
cn =< φn |Ψ(x0 , 0) >= φn (x0 ))Ψ(x0 , 0)dx0
−∞

on a donc
∞
X En
|Ψ(x, t) >= cn φn (x)e−i ~
t

n=0
12.2 “Diffusion Monte-Carlo” 85

En combinant ces deux dernières équations, on obtient après permutation des signes somme
et intégrale : "∞ #
Z ∞ X En
|Ψ(x, t) >= dx0 φ∗n (x)e−i ~ t φn (x0 ) Ψ(x0 , 0)
−∞ n=0

On obtient donc |Ψ(x, t) > sous forme d’une équation intégrale :

Z ∞
|Ψ(x, t) >= dx0 K(x, t|x0 , 0)Ψ(x0 , 0)
−∞

où le noyau :
∞
X En
K(x, t|x0 , 0) = φ∗n (x)e−i ~
t
φn (x0 )
n=0

appelé opérateur d’évolution représente la probabilité pour une particule se trouvant en x 0

à l’instant t = 0 de passer en x à l’instant t.
On peut écrire formellement :
∞
X t
K(x, t|x0 , 0) = φ∗n (x)e−i ~ H φn (x0 )
n=0

Pour résoudre l’équation intégrale par Monte-Carlo, on a besoin de fonctions réelles posi-
tives, et on va se débarrasser du facteur imaginaire par le changement de variable :

τ = it

D’autre part on va effectuer une translation arbitraire de l’origine des énergies :

En → E n − E R

V (x) → V (x) − ER
L’équation de Schrödinger devient :

∂Ψ(x, τ ) ~2 ∂ 2 Ψ(x, τ )
~ = − [V (x) − ER ]Ψ(x, τ )
∂τ 2m ∂x2
Le premier terme du second membre correspond à une équation de diffusion où Ψ représente
la densité de particules qui diffusent. Une telle équation peut être simulée avec une marche
aléatoire de particules dans l’espace de configuration.
Le second terme est semblable à une équation décrivant un processus de mort et de naissance
d’individus dans une population.
L’équation globale peut être simulée par la combinaison d’un processus de diffusion et
d’un processus de branchement, dans lequel le nombre de particules diffusant augmente ou
diminue de façon à réduire la densité de probabilité dans les régions où V (x) est grand et
à l’augmenter dans les zones d’énergie potentielle favorables.
86 Monte-Carlo quantique

On a :
∞
X (En −ER )
|Ψ(x, τ ) >= cn φn (x)e− ~
τ

n=0

et le comportement asymptotique de |Ψ(x, τ ) > pour τ → ∞ ne dépend que du signe de

(E0 − ER ), E0 représentant l’energie fondamentale :
– ER > E0 : |Ψ(x, τ ) >→ 0 pour τ → ∞
– ER = E0 : |Ψ(x, τ ) >→ c0 |φ0 (x) > pour τ → ∞
– ER < E0 : |Ψ(x, τ ) >→ ∞ pour τ → ∞
Sous forme intégrale, Ψ(x, τ ) obéit à l’équation :
Z ∞
|Ψ(x, τ ) >= dx0 K(x, τ |x0 , 0)Ψ(x0 , 0) (12.1)
−∞

avec
∞
X (En −ER )
K(x, t|x0 , 0) = φ∗n (x)e− ~
τ
φn (x0 )
n=0
que l’on peut écrire formellement :
∞
X (T̂ +V̂ )
K(x, t|x0 , 0) = φ∗n (x)e− ~
τ
φn (x0 )
n=0

avec
~2 ∂ 2 Ψ(x, τ )
T̂ =
2m ∂x2
et
V̂ = −[V (x) − ER ]
Les opérateurs T̂ et V̂ ne commutent pas, donc en général :
(T̂ +V̂ ) T̂ V̂
e− ~
τ
6= e− ~ τ e− ~ τ

Cependant, pour δτ petit, on a [formule de Trotter-Suzuki] :

(T̂ +V̂ ) V̂ δτ T̂ V̂ δτ
e− ~
δτ
= e− ~ 2 e− ~ τ e− ~ 2 + O(δτ 3 ) (12.2)

L’opérateur :
∞
X V̂
Kdif f (x, t|x0 , 0) = φ∗n (x)e− ~ τ φn (x0 )
n=0

est connu exactement car il correspond à la propagation d’une particule libre (potentiel
constant).
Considérons une particule de masse m dans une boite unidimensionnelle de longueur L.
Les fonctions propres sur l’intervalle [-L/2,L/2] sont

1 2πn ~2 2
φ n = √ e kn x , avec kn = et En = k
L L 2m n
12.2 “Diffusion Monte-Carlo” 87

Ce qui donne :
∞ 2
X 1 kn
Kdif f (x, t|x0 , 0) = √ e−ikn (x−x0 ) e−it~ 2m
n=0
L
soit en passant à la limite continue, avec L → ∞ :
Z ∞ 2
1 kn
Kdif f (x, t|x0 , 0) = e−ik(x−x0 ) e−it~ 2m dk
2π −∞

soit après intégration et substitution de it par τ :

h m i1/2 m 2
Kdif f (x, τ |x0 , 0) = e− 2~τ (x−x0 ) (12.3)
2π~τ

On va donc diviser τ en N petits intervalles δτ = τ /N et pour δτ petit approximer

K(x, t|x0 , 0) en utilisant la relation de Trotter-Suzuki :
h m i1/2 δτ (V (x)−ER ) m 2 δτ (V (x0 )−ER )
K(x, δτ |x0 , 0) = e− ~ 2 e− 2~δτ (x−x0 ) e− ~ 2 (12.4)
2π~τ
Et en itérant N fois la relation intégrale 12.1 on obtient :
 (
N −1 Z ∞ N
)
Y Y
Ψ(x, τ ) =  dxj  W (xn )P (xn , xn−1 ) Ψ(x0 , 0)
j=0 −∞ n=1

avec
m  21 m(xn − xn−1 )2
  
P (xn , xn−1 ) = exp −
2π~δτ 2~δτ
et  
(V (xn ) − ER )δτ
W (xn ) = exp −
~
La fonction : P (xn , xn−1 ) représente une densité de probabilité car :
Z ∞
P (x, y)dy = 1 ∀x
−∞

ceci va nous permettre de construire une chaine de Markov ayant pour loi de transition
P (x, y).
Les poids W (xn ) dépendant uniquement du potentiel V (x) et de l’énergie de référence E R
vont être traduits par un processus de naissance-mort de particules. C’est à dire qu’au lieu
de raisonner sur une seule trajectoire de la particules, on va suivre en parallèle l’évolution
de N trajectoires.
– On part de N points distribués suivant une probabilité proportionnelle à |Ψ(x 0 , 0)|2
– A chaque étape τn = nδτ on tire N nouvelles abscisses suivant la loi de probabilité
P (xn , xn−1 )
88 Monte-Carlo quantique

– On définit mn = min[int(W (xn ) + η), ν] où η est un nombre aléatoire uniformément

distribué sur [0,1], ν est un petit nombre arbitraire (en pratique ν ≈ 3) Si m n = 0, on
arrête le processus ; si mn = 1 la particule continue normalement son chemin ; enfin si
mn > 1 on crée (mn − 1) nouvelles particules fictives au même endroit.
Au cours de l’évolution, on ajuste E R pour obtenir une “population” stable, on a alors

E0 = ER =< V >

Donc, en pratique on ajuste à chaque étape (E R )n à < V >n Lorsque le processus a

convergé, la densité de points sur l’axe x représente le module au carré de la fonction
d’onde correspondant à l’état fondamental.
Le rôle de ν est d’éviter au début un trop grand nombre de branchements. Il n’intervient
pas dans le résultat final, car lorsque l’équilibre est atteint, la probabilité pour que W (x)
soit supérieur à ν est faible.

12.3 “Path Integral Monte-Carlo”

12.3.1 Valeur moyenne d’une grandeur physique

En mécanique quantique, un système est décrit par un hamiltonien H. Soit |Φ i > un

ensemble de vecteurs propres orthonormés associés respectivement aux valeur propres E i :

H|Φi >= Ei |Φi >

δτ

2δτ

3δτ

4δτ
X

Fig. 12.1 – Processus de diffusion + naissance-mortalité

12.3 “Path Integral Monte-Carlo” 89

Une grandeur physique est représentée par un opérateur O. Dans l’ensemble canonique la
probabilité du système de se trouver dans l’état |Φ i > est proportionnelle à exp(−Ei /kT )
et la valeur moyenne de la grandeur physique correspondante est :
P
< Φi |O|Φi > exp[−βEi ]
< O >= i
Z
X
Z= exp[−βEi ]
i

est la fonction de partition et β = 1/kT

Cette valeur moyenne peut encore s’écrire sous la forme plus générale, indépendante de la
représentation (base othonormée de l’espace de Hilbert) choisie :

T r O exp[−βH]
< O >= (12.5)
Z
avec Z = T r exp[−βH]. La trace “T r” d’un opérateurs étant invariante, quelle que soit la
représentation choisie..

12.3.2 La Matrice Densité

On choisit, en général pour représentation non pas une base de vecteurs propres telle que
{|Φi >} mais une base d’états orthonormés {|R >} décrivant les positions dans l’espace
des particules. Cette base est continue et les sommes discrètes relatives à la base discrète
{|Φi >} considérée précédemment sont remplacées par des intégrales. En particulier, la
trace d’un opérateur, somme de ses valeurs propres de vient une intégrale :
Z
T r A = dR0 < R0 |A|R0 >

Alors la relation générale 12.5 s’écrit :

dR0 < R0 |O exp(−βH)|R0 >

R
< O >= (12.6)
Z
En utilisant l’orthonormalité de la base {|R >} qui se traduit par la relation “de fermeture”
Z
dR|R >< R| = 1 (12.7)

on peut rendre plus explicite le produit des deux matrices O et exp(−βH) pour arriver à
l’expression :
dRdR0 < R0 |O|R >< R| exp(−βH)|R0 >
R
< O >= (12.8)
Z
En introduisant la “Matrice Densité :

ρ(R, R0 ; β) =< R| exp(−βH)|R0 > (12.9)

90 Monte-Carlo quantique

on a
dRdR0 < R0 |O|R > ρ(R, R0 ; β)
R
< O >= (12.10)
Z
avec Z
Z= dRρ(R, R, β)

Pour calculer l’intégrale ci-dessus par la méthode de Monte-Carlo, nous allons donc générer
de configurations {R, R 0 } suivant une probabilité proportionnelle à ρ(R, R 0 ; β) et la moyenne
des élements de matrice < R 0 |O|R > correspondant à ces configurations nous donnera accès
à une approximation de < O >

12.3.3 Echantillonnage de ρ(R, R0 ; β)

Une particule à une dimension

Nous illustrons d’abord la méthode sur le cas simple utilisé précedemment pour l’algorithme
de “Diffusion Monte-Carlo” : une particule à une dimension, de masse m, soumise à un
potentiel V (x). L’hamiltonien s’écrit :

~2 ∂ 2
H =− + V (x)
2m ∂x2
soit
H = T̂ + V̂
en séparant la partie cinétique :
~2 ∂ 2
T̂ = −
2m ∂x2
de la partie potentielle
V̂ = V (x)
La matrice densité sécrit :

ρ(x, x0 ; β) =< x| exp[−β(T̂ + V̂ )]|x0 >

Cette quantité est très difficile à calculer, alors que les quantités < x| exp[−β T̂ ]|x00 > et
< x00 | exp[−β V̂ ]|x0 > sont relativement simples. Malheureusement, les opérateurs T̂ et V̂
ne commutent pas et
exp[−β(T̂ + V̂ )] 6= exp[−β T̂ ] exp[−β V̂ ]
Par contre, en divisant β en M intervalles δτ = β/M , on pourra recourir à la formule de
Trotter-Suzuki (cf. Equation 12.2)
h iM
exp[−β(T̂ + V̂ )] = lim exp(−δτ T̂ ) exp(−δτ V̂ ) (12.11)
M →∞
12.3 “Path Integral Monte-Carlo” 91

Enfin, en intercalant des relations de fermeture :

Z
dxi |xi >< xi | = 1

entre chaque produit d’opérateurs, on aboutit à la relation :

ρ(x0 , xM ; β) =
Z Z Z
dx1 dx2 · · · dxM −1 < x0 | exp(−δτ T̂ ) exp(−δτ V̂ )|x1 >< x1 | exp(−δτ T̂ ) exp(−δτ V̂ )|x2 >

· · · < xM −1 | exp(−δτ T̂ ) exp(−δτ V̂ )|xM >

Calcul de < xi−1 | exp(−δτ T̂ ) exp(−δτ V̂ )|xi >

L’opérateur V̂ est diagonal, on a donc :

< xi−1 | exp(−δτ T̂ ) exp(−δτ V̂ )|xi >=< xi−1 | exp(−δτ T̂ )|xi > exp[−δτ V̂ (xi )]

Mais l’élément de matrice :

h m i1/2 m
< xi−1 | exp(−δτ T̂ )|xi >= 2
exp[− 2 (xi − xi−1 )2 ]
2π~ δτ 2~ δτ
pour une particule libre a été calculé au paragraphe précédent (cf. Relation 12.3)
En reportant dans l’expression précédente, le calcul de ρ(x 0 , xM ; β) se présente comme une
intégrale à M − 1 dimensions :
M h
!
m
Z Z Z X i
2
ρ(x0 , xM ; β) = C dx1 dx2 · · · dxM −1 exp − (xi−1 − xi ) + δτ V (xi )
2~2 δτ
i=1
(12.12)
avec h m iM/2
C=
2π~2 δτ
La somme figurant en exposant du facteur exponentiel correspond à l’énergie d’une chaı̂ne
de particules de coordonnées x0 , x1 , ... ,xM , chaque particule étant reliée à la précédente par
un ressort d’energie proportionnelle au carré de la distance qui les sépare. Chaque particule
i étant soumise à un potentiel V (x i ). On obtient l’image d’un polymère (cf section 11.2.2)

N particules discernables à 3 dimensions

Ces calculs se généralisent facilement pour une particule à deux ou trois dimensions... et
pour un ensemble de N particules à trois dimensions soumises à un potentiel dı́nteraction
de paires
V (|Ri − Rj |
92 Monte-Carlo quantique

ne dépendant que de la distance entre paires de particules. Dans ce cas, chaque particule
p est représentée par une chaı̂ne de M monomères de coordonnées R p,i . La matrice densité
s’exprime alors comme une intégrale à 3N M dimensions :
3N M/2 R
ρ(R1,0 , R2,0 , ..., RN,0 , R1,M , R2,M , ..., RN,M ; β) = 2π~m2 δτ
 R
dR1,1 · · · dRN,M
h PN PM i h PN P PM i
exp − 2~m 2 δτ p=1 i=1 (R p,i−1 − R p,i ) 2 exp −δτ
p=1 k>p i=1 V (|R p,i − R k,i |

N Bosons

Au paragraphe précédent, nous avons raisonné sur des particules discernables. Or, en
mécanique quantique, les particules sont indiscernables.
Pour des bosons, l’expression de la matrice densité ρ B doit être complètement symétrique
par rapport à une permutation quelconque des particules.
On doit donc sommer sur toutes les permutations possibles de l’expression d’une matrice
densité de particules discernables :
1 X
ρB (R1,0 , ..., RN,0 , R1,M , ..., RN,M ; β) = ρB (R1,0 , ..., RN,0 , RP (1),M , ..., RP (N ),M ; β) =
N!
P

où P (1), P (2), ... ,P (N ), représente une permutation quelconque des N particules.
Les permutations seront elles-mêmes échantillonnées par un algorithme de Monte-Carlo.
Nous n’aborderons pas le problème des Fermions pour lesquels l’antisymétrisation conduit
à une alternance de signe qui constitue un des problèmes les plus difficiles pour l’application
de la méthode de Monte-Carlo.

Pour en savoir plus...

Le livre de Feynman [16] est essentiel pour approfondir la formulation de la mécanique

quantique sous forme “d’intégrales de chemin”, que nous avons brièvement introduite ci-
dessus.
Pour une approche plus complète des méthodes de Monte-Carlo quantique, nous conseillons
les articles de revue et cours de Ceperley [17]
Annexes : Probabilités et
Physique Statistique
94 Monte-Carlo quantique
 Rappels de théorie des
A probabilités

Nous nous contenterons de rappeler brièvement les notions indispensables à la justification

des méthodes de Monte-Carlo. Pour un exposé complet et didactique sur la théorie des
probabilités, nous conseillons le livre (niveau second cycle) de J. Bass[18] ou le cours de
l’Ecole Polytechnique de J. Neveu[19].

A.1 Expérience, ensemble des résultats possibles.

Une expérience aléatoire est une expérience dont le résultat est soumis au hasard. Elle se
décrit par la donnée de l’ensemble des résultats possibles. Nous noterons ω un tel résultat
et Ω l’espace formé par tous les résultats possibles.

Exemple : l’espace Ω associé au jet aléatoire de deux dés est composé des 36 couples d’entiers
ω = (x, y) tels que 1 ≤ x, y ≤ 6

A.2 Notion d’événement aléatoire.

Un événement aléatoire est représenté par l’ensemble des résultats ω de l’expérience qui le
réalisent. Dans l’exemple précédent du jet de deux dés, on peut définir comme événement :
“obtenir un total de points supérieur à 8”. Cet événement est représenté par le sous ensemble
de résultats A ⊂ Ω :

A = {(3, 6), (4, 5), (5, 4), (6, 3), (4, 6), (5, 5), (6, 4), (5, 6), (6, 5), (6, 6)}

Nous devons maintenant envisager des opérations logiques sur ces événements.
– A tout événement A est associé son contraire A c représenté par l’ensemble complémentaire
de A dans Ω. (Le contraire de A dans l’exemple précédent est “obtenir un total de points
inférieur ou égal à 8”).
– Pour tout couple dévénements (A 1 , A2 ), l’événement “A1 et A2 ” est par définition celui
qui est réalisé si les événements A 1 et A2 sont réalisés à la fois. En langage de théorie des
96 Rappels de théorie des probabilités

ensembles, il correspond à l’intersection A 1 ∩ A2 . En prenant toujours le même exemple

de jet de deux dés, considérons, A 1 : “obtenir un total de points supérieur à 8” et
A2 : “chacun des dés est tombé sur un nombre pair”, l’événement “A 1 et A2 ” est donc
représenté par :
A3 = A1 ∩ A2 = {(4, 6), (6, 4), (6, 6)}
– L’événement impossible sera représenté par l’ensemble vide ∅. L’équation A 1 ∩ A2 = ∅
signifie que les événements A1 et A2 sont incompatibles, c’est à dire que les parties de Ω
qui les représentent sont disjointes
– Pour tout couple dévénements (A 1 , A2 ), l’événement “A1 ou A2 ” est par définition celui
qui est réalisé si l’un au moins des deux événements A 1 ou A2 est réalisé (“ou” non
exclusif). En langage de théorie des ensemble, il correspond à l’union A 1 ∪ A2 .
– L’événement certain sera noté Ω puisqu’il est réalisé quel que soit le résultat ω de l’-
expérience.
En conclusion, une expérience aléatoire est décrite mathématiquement par la donnée d’un
espace Ω et d’une classe A de parties de Ω. Chaque élément A de A représente un événement
que l’on considère relativement à l’expérience. Nous exigerons de la classe A les propriétés
suivantes :
– Elle contient l’événement certain : l’ensemble Ω lui-même.
– Si A appartient à A, l’ensemble Ac complémentaire de A par rapport à Ω (i.e. “le contraire
de A”) est aussi un élément de la classe A
– Si A1 et A2 appartiennent à A, la réunion et l’intersection de A 1 et A2 appartiennent à
A.
– Lorsque l’espace Ω n’est pas fini, il est encore indispensable d’exiger que la réunion d’une
infinité dénombrable d’ensembles de A appartienne à A.
On dit que A est une “tribu de parties de Ω”, ou encore une “σ-algèbre construite sur Ω.”

Par exemple, dans R, la tribu engendrée par les segments ouverts [i.e la plus petite tribu
contenant les segments ouverts] est nommée “tribu Borélienne”. D’après les propriétés
énoncées ci-dessus, elle contient aussi les segments fermés (complémentaires d’ouverts),
l’intersection et la réunion de nombres infinis dénombrables d’entre eux.

A.3 Notion de probabilité.

A.3.1 Définition heuristique

Une propriété essentielle d’une expérience est de pouvoir être répétée indéfiniment. Si
nous comptons, au cours de N répétitions d’une expérience, le nombre de fois N A où
un événement A est réalisé, nous observons que la fréquence statistique N A /N tend vers
une limite p ≤ 1 lorsque N → ∞. D’où l’idée d’associer un nombre positif ou nul P (A) = p,
inférieur à 1, à tout événement A.
Si deux événements sont incompatibles (A 1 ∩A2 = ∅, la fréquence statistique de l’événement
A.4 Probabilités conditionnelles. Evénements indépendants 97

A1 ∪ A2 au cours de N répétitions d’un expérience aléatoire est égale à la somme des

fréquences statistiques des événements A 1 et A2 puisque un seul des événement, au plus,
peut être réalisé.
NA1 ∪A2 N A1 N A2
= +
N N N
donc
P (A1 ∪ A2 ) = P (A1 ) + P (A2 ) si A 1 ∩ A2 = ∅

A.3.2 Définition mathématique

– A chaque événement A ∈ A est associé un nombre réel P (A) ∈ [0, 1] appelé probabilité
de A
– La probabilité 1 est attribuée à l’événement certain : P (Ω) = 1
– Nous imposons la propriété d’additivité pour toute suite (finie ou infinie-dénombrable)
d’événements deux à deux indépendants :

P (A1 ∪ A2 ∪ · · · ∪ An ) = P (A1 ) + P (A2 ) + · · · + P (An ) (A.1)

si
Ai ∩ A j = ∅ ∀(i, j), i 6= j
En langage mathématique, P est une “mesure positive, de masse totale P (Ω) = 1, définie
sur la tribu A de parties de Ω”.

On note “espace probabilisé” {Ω, A, P } le triplet formé par l’ensemble des résultats possi-
bles de l’expérience Ω, la tribu de parties de cet ensemble A et la mesure de probabilités
P.

A.4 Probabilités conditionnelles. Evénements indépendants

Considérons une expérience qui se réalise en deux temps. En reprenant toujours l’exemple
du jet de deux dés, on peut lancer successivement un premier dé puis le second. On peut
définir des événements relatifs aux résultats du premier temps (exemple B : “le résultat
du premier dé est pair”) et des événement relatifs au résultat global après le second temps
(exemple A : “le total des points est 8”). Ces deux événements, bien sûr, appartiennent à
une même tribu A de parties de Ω
Nous devons alors distinguer heuristiquement entre la probabilité, dans l’absolu, P (A)
de l’événement A, et la probabilité de l’événement A, sachant que B s’est produit, que
nous appellerons “probabilité de l’événement A conditionnelle en l’événement B”. Nous la
noterons P (A/B).
De même que la fréquence NA /N des réalisations de l’événement A au cours de N répétitions
de l’expérience aléatoire, donne une estimation statistique de la probabilité P (A) de l’événe-
ment A, nous pouvons penser que le rapport N A∩B /NB de réalisations de l’événement A
98 Rappels de théorie des probabilités

parmi les NB répétitions de l’expérience qui ont donné B fournisse une estimation statis-
tique de la probabilité de l’événement A conditionnelle en l’événement B.
Après avoir remarqué que :
NA∩B
NA∩B N
= NB
NB N
il devient naturel de poser comme définition :

P (A ∩ B)
P (A/B) = (A.2)
P (B)

Deux événement A et B seront dits indépendants si la probabilité de l’événement A condi-

tionnelle en l’événement B ne dépend pas de l’événement B, c’est à dire si P (A/B) = P (A),
ce qui d’après la relation précédente correspond à l’égalité

P (A ∩ B) = P (A)P (B) (A.3)

A.5 Variable aléatoire.

Une variable aléatoire est définie par référence à une expérience aléatoire comme un nombre
réel, ou plus généralement un vecteur X = φ(ω) dont la valeur dépend du résultat ω de
l’expérience. Par exemple, dans l’expérience du jet de deux dés, la somme S des points
amenés par les deux dés est une variable aléatoire réelle. Le vecteur V = (m, n) où m est le
nombre de points amené par le premier dé et n le nombre de points amenés par le second
dé est un vecteur aléatoire qui prend ses valeurs dans R 2 .
Une variable aléatoire est donc représentée par une application φ de Ω dans R d . Elle est
“aléatoire” du fait du caractère aléatoire de l’expérience, dont chaque résultat possible est
un point ω de l’ensemble Ω.
Nous exigerons une condition essentielle sur cette application :
L’ensemble des résultats ω de l’expérience, tels que :

a < X = φ(ω) < b (A.4)

appartient à la tribu A

Si X est une variable aléatoire réelle, a et b représentent deux réels quelconques tels que
a < b.
Si X est une variable aléatoire vectorielle à valeurs dans R d , a = (a1 , a2 , · · · , ad ) et b =
(b1 , b2 , · · · , bd ) représentent deux vecteurs de R d et les inégalités sont vérifiées composante
par composante. Dans ce cas, l’inégalité

a < x < b i.e. a i < αi < bi ∀i

A.6 Loi de probabilité d’une variable aléatoire 99

définit un pavé B de Rd .

Donc, l’image réciproque par l’application X = φ(ω) d’un segment B de R, pour une
variable aléatoire réelle, (ou plus généralement d’un “pavé” B de R d , pour une variable
aléatoire vectorielle) est un élément (i.e “un événement”) de la tribu A de sous ensembles
de Ω.
Cette condition est indispensable pour que nous puissions parler de la probabilité de
“l’événement : {a < X = φ(ω) < b}. Elle correspond à exiger que l’application X = φ(ω)
soit une “application mesurable”.

Remarque importante : la classe B de parties de R d formée par les pavés de R d , leurs

intersections et la réunion d’un nombre fini ou infini-dénombrable d’entre eux constitue
une tribu (cf. ci-dessus) appelée “tribu Borélienne”

A.6 Loi de probabilité d’une variable aléatoire

Soit AB l’image réciproque dans Ω d’un pavé B de R d (–B est un élément de la tribu
Borélienne B–) par l’application X = φ(ω), c’est à dire :
AB = {ω tels que X = φ(ω) ∈ B}
alors les sous ensembles “événements” A B forment une sous tribu de la tribu A appelée “la
tribu engendrée par la variable aléatoire X.

La loi :
P(X) (B) = P (AB )
définit une probabilité P(X) sur Rd par rapport à la tribu B “Borélienne”) formée par les
pavés ouverts ou fermés de R d et la réunion d’un nombre fini, ou infini-dénombrable d’entre
eux.

Cette probabilité P(X) sur Rd , image de la probabilité P par l’application X = φ(ω) de Ω

dans Rd est appelée “la loi de probabilité de la variable aléatoire X”

L’intérêt de l’introduction de ces “lois de probabilité” tient à ce qu’il est plus aisé de
manipuler des “mesures” sur R d que sur l’espace Ω. Pour les applications en physique,
nous utiliserons le plus souvent la mesure de Lebesgue ou la mesure de Rieman sur R d .

A.7 Fonction de répartition d’une variable aléatoire réelle

Considérons le segment semi infini B =] − ∞, z] alors

F (z) = PX (B) = P (AB ) = P ({−∞ < X(ω) ≤ z})
100 Rappels de théorie des probabilités

est la fonction de partition relative à la loi de probabilité de la variable aléatoire X. C’est

une fonction positive, croissante qui atteint sa valeur maximale pour z → ∞.

lim F [z] = P (Ω) = 1

z→∞

C’est une fonction à valeurs dans le segment [0, 1]

A.8 Densité de probabilité.

Si la loi de probabilité P(X) de la variable aléatoire X peut s’écrire :

Z
P(X) (B) = P ({X = φ(ω) ∈ B} = f (x)dx (A.5)
B

où dx désigne la mesure de Lebesgue ou la mesure de Rieman.

Alors nous dirons que la variable aléatoire X possède la densité de probabilité f (x)

(f est la densité de la “mesure de probabilité” P (X) )

f est une fonction positive, et il résulte de l’égalité P (X) (Rd ) = 1 que :

Z
f (x)dx = 1 (A.6)
Rd

Si une variable aléatoire X possède une densité de probabilité f (x), on peut écrire si X est
définie dans R (d = 1) :

P ({x < X = φ(ω) < x + dx}) = f (x)dx lorsque dx → 0 (A.7)

et plus généralement pour une variable aléatoire définie dans R d , en notant x = (α1 , α2 , · · · , αd )
un élément de Rd :

P ({αk < (X)k = φ(ω)k < αk + dαk }) = f (x)dα1 dα2 · · · dαd (A.8)

Etant donné que, dans les applications physiques que nous allons étudier, nous nous restrein-
drons au cas où les variables aléatoires considérées possèdent une densité de probabilité f ,
le lecteur, quelque peu dérouté par le formalisme général peut se “raccrocher” aux relations
ci-dessus en tant que définition d’une densité de probabilité.

la fonction de répartition d’une variable aléatoire réelle admettant une densité de probabilité
f (x) est Z z
F (z) = f (x)dx
−∞
A.9 Espérance, moments d’ordre n, variance 101

Quelques exemples de densités de probabilité dans R

– La loi uniforme sur un segment :
f (x) = 1/(b − a) si a<b<x et 0 sinon
– La loi exponentielle
f (x) = λe−λx x≥0 et 0 sinon
qui décrit un certain nombre de phénomènes physiques, dont le temps de désintégration
d’un élément radioactif.
– La loi de Gauss (ou loi “normale”)
1 (x−m)2
f (x) = √ e− 2σ2 (A.9)
2πσ
très importante en physique et statistique

A.9 Espérance, moments d’ordre n, variance

A.9.1 Valeur moyenne ou “Espérance” d’une variable aléatoire réelle

Considérons tout d’abord une variable aléatoire discrète X, c’est à dire à valeurs dans
l’ensemble des entiers relatifs Z ⊂ R. Au cours de N répétitions de l’expérience aléatoire,
chaque valeur {X = xi } “sort” ki fois. Si nous voulons évaluer empiriquement une valeur
moyenne m de la variable aléatoire X, nous prendrons :
X X  ki 
m= ki xi /N = xi
N
i i

Lorsque N → ∞, on associe la probabilité p i à la fréquence ki /N , d’où la définition de la

valeur moyenne ou espérance dans le cas discret :
X
E[X] = m = xi pi (A.10)
Z

Dans le cas continu, la somme discrète est remplacée par une intégrale :
Z Z
E[X] = XdP = xdP(X) (A.11)
Ω R
où les intégrales sont prises au sens général de la théorie de l’intégration par rapport aux
mesures abstraites dP , dP(X) .

Dans le cas où la variable aléatoire X possède la densité de probabilité f , alors, on peut
écrire plus simplement : Z
E(X) = m = xf (x)dx (A.12)
R
où dx représente la mesure de Lebesgue ou la mesure de Rieman sur R
102 Rappels de théorie des probabilités

A.9.2 Moment d’ordre n. Variance d’une variable aléatoire réelle

Toute fonction Y = g(X) = g[φ(ω)] = g ◦ φ(ω) définit une nouvelle variable aléatoire. Son
espérance est : Z
E[Y ] = E[g(X)] = g(x)f (x)dx (A.13)
R

Considérons en particulier g(X) = X n , alors

Z
µn = E[X n ] = xn f (x)dx (A.14)
R

est le moment d’ordre n de la variable aléatoire réelle X

Le moment d’ordre deux : Z

2
E[X ] = x2 f (x)dx (A.15)
R

joue un rôle particulièrement important.

On appelle variance le moment d’ordre deux de (X − m) = (X − E[X]) et écart type σ la
racine carrée de cette quantité
p
σ= E[(X − E[X])2 ] (A.16)

La variable aléatoire (X − E[X])/σ a pour espérance 0 et pour écart type 1. On dit qu’elle
est réduite.

A.10 Fonction Caracteristique

A.10.1 Définition

La fonction caractéristique d’une variable aléatoire réelle X est la transformee de Fourier

de la loi de probabilité P(X) de la variable aléatoire X :

F(t) = E[exp (itX)]

Son développement de Taylor à l’origine s’exprime en fonction des moments d’ordre n

successifs de la variable aléatoire X

F(t) = 1 + itE(X) − t2 E[X 2 ]/2 − it3 E[X 3 ]/3!...

Si X admet une densité de probabilité f (x), alors sa fonction caractéristique est la trans-
formée de Fourier de la fonction f (x).
A.10 Fonction Caracteristique 103

A.10.2 Exemple :

Loi gaussienne :
Nous choisissons pour exemple la loi gaussienne de densité de probabilité :

x2
 
1
fG (x) = √ exp − 2
2πσ 2 2σ

car sa fonction caractéristique est très simple, et le résultat nous servira à la démonstration
du théorème de la limite centrale.
Z ∞
x2
 
1
FG (t) = √ exp (itx) exp − 2 dx
2πσ 2 −∞ 2σ

soit
∞  2 2
(x − σ 2 it)2
 
1 σ t
Z
FG (t) = √ exp − exp − dx
2πσ 2 −∞ 2 2σ 2
Nous allons montrer que l’intégrale :
∞
(x − σ 2 it)2
 
1
Z
I=√ exp − dx
2πσ 2 −∞ 2σ 2

est simplement 1. Pour cela nous considérons dans le plan complexe la fonction générale :

Z2
 
1
G(Z) = √ exp − 2
2πσ 2 2σ

et nous l’intégrons sur le contour fermé représenté sur la figure A.1, Puisqu’il n’y a aucun
pôle à l’intérieur de ce contour, le théorème des résidus nous dit que l’intégrale sur le
contour est nulle :
Z Z Z Z
G(Z)dZ + G(Z)dZ + G(Z)dZ + G(Z)dZ = 0
C1 C2 C3 C4

Lorsque A → ∞ les deux derniers termes tendent vers 0 car |Z| 2 → ∞. Le second terme
n’est autre que I et le premier terme représente l’intégrale sur l’axe réel :
−∞
x2
 
1
Z
√ exp − 2 dx = −1
2πσ 2 ∞ 2σ

Donc I = 1 et la fonction caractéristique de la loi gaussienne est :

σ 2 t2
 
FG (t) = exp − (A.17)
2
104 Rappels de théorie des probabilités

-A 0 A
C1
C3 C4
C2

−σ t i
2

Fig. A.1 –

A.11 Inégalité de Bienaymé-Tchebichef

Nous allons démontrer une inégalité triviale, mais néanmoins importante introduite par
Bienaymé et Tchebichef.

Pour toute variable aléatoire réelle X, de carré intégrable (i.e µ 2 = E[X 2 ] est fini), et pour
tout nombre réel a strictement positif :
µ2
P (|X| ≥ a) ≤ 2 (A.18)
a

Démonstration : Désignons par

1{|x|≥a}
la fonction qui vaut 1 si |x| ≥ a et 0 sinon, alors on a
1 2
1{|x|≥a} ≤ x
a2
et on a aussi, entre les variables aléatoires définies par les fonctions correspondant à chacun
des deux membres :
1
1{|X|≥a} ≤ 2 X 2
a
En prenant l’intégrale par rapport à la mesure de probabilité P de chacun des deux mem-
bres, on obtient :  2
X
P (|X| ≥ a) ≤ E
a2
A.12 Couples de variables aléatoires 105

Et comme
X2
 
E = E[X 2 ]/a2
a2
l’inégalité de Bienaymé-Tchebichev est ainsi démontrée

A.12 Couples de variables aléatoires

Soient deux variables aléatoires X 1 et X2 de lois de probabilité PX1 et PX2 . On peut définir
la loi de probabilité conjointe du couple (ou vecteur) U = (X 1 , X2 ) :

P(U ) (A × B) = P ({X1 ∈ A} ∩ {X2 ∈ B})

A.12.1 Covariance

On appelle covariance des deux variables aléatoires :

cov(X1 , X2 ) = E[(X1 − E[X1 ])(X2 − E[X2 ])] = E[X1 X2 ] − E[X1 ]E[X2 ]

A.12.2 Indépendance

Deux variables aléatoires sont indépendantes si quels que soient A et B :

P(U ) (A × B) = P(X1 ) (A)P(X2 ) (B)

soit
P ({X1 ∈ A} ∩ {X2 ∈ B}) = P ({X1 ∈ A})P ({X2 ∈ B})

On en déduit :

E(X1 X2 ) = E(X1 )E(X2 ) donc σ(X1 , X2 ) = 0

Si les variables alátoires X1 et X2 admettent des densités de probabilité f 1 (x1 ), f2 (x2 ) alors
la loi de probabilité conjointe du couple admet pour densité de probabilité :

h(x1 , x2 ) = f1 (x1 )f2 (x2 )

A.12.3 Somme de 2 variables aléatoires indépendantes

La fonction caractéristique de la somme S = X 1 + X2 de deux variables aléatoires X1 et

X2 s’écrit :
F(t) = E[exp it(X1 + X2 )] = E[exp(itX1 ) exp(itX2 )]
106 Rappels de théorie des probabilités

Si deux variables aléatoires X1 et X2 sont indépendantes, les variables aléatoires exp(itX 1 )

et exp(itX2 ) sont indépendantes et l’espérance du produit de ces dernières est égale au
produit des espérances. Donc

F(t) = E[exp(itX1 )]E[exp(itX2 )]

si X1 et X2 admettent pour densités de probabilité f 1 (x1 ) et f1 (x1 ) alors F(t), transformée

de fourier de la densité de probabilité de la somme X 1 + X2 est le produit des transformées
de Fourier des densités de probabilité f 1 et f2 . On en déduit que la densité de probabilité
de X1 + X2 est le produit de convolution :
Z
f1 ∗ f2 (x) = f1 (y)f2 (x − y)dy

Ce résultat peut aussi est obtenu de manière élémentaire : L’ensemble des valeurs du couple

X2
s+
ds
s

s-x +ds
s-x

x X1

Fig. A.2 –

U = (X1 , X2 ) étant représenté par le plan (x1, x2) de la figure A.2 on a :

P ({x1 < X1 < x1 + dx1 } ∩ {x2 < X2 < x2 + dx1 }) = h(x1 , x2 )dx1 dx2

Puisque les variables sont indépendantes, la densité de probabilité h(x 1 , x2 ) du couple

(X1 , X2 ) est le produit des densités de probabilité des variables X 1 et X2 soit :

P ({x1 < X1 < x1 + dx1 } ∩ {x2 < X2 < x2 + dx1 }) = f1 (x1 )f2 (x2 )dx1 dx2
A.13 Généralisation à N variables aléatoires 107

On a alors : Z
P ({s < S < s + ds}) = f1 (x1 )f2 (x2 )dx1 dx2
[aire en gris]

soit Z ∞ Z s−x+ds
P ({s < S < s + ds}) = f1 (x)dx f2 (x2 )dx2
−∞ s−x
Et dans la limite où ds → 0 :
Z ∞ 
P ({s < S < s + ds}) = f1 (x)f2 (s − x))dx ds
−∞

A.13 Généralisation à N variables aléatoires

Tous les raisonnements précédents pour un couple de deux variables aléatoires se généralisent
à un n-uplet U = (X1 , X2 , · · · , Xn ) (ou vecteur à n dimensions) de variables aléatoires.
Nous dirons que les n variables X1 , X2 , ..., Xn sont indépendantes si pour tout ensemble
A1 , A2 , ... An de segments ouverts de R, la loi conjointe :

P(U ) (A1 × A2 × · · · × An ) = P(X1 ) (A1 )P(X1 ) (A1 ) · · · P(Xn ) (An )

soit
P ({X1 ∈ A1 } ∩ {X2 ∈ A2 } ∩ · · · ∩ {X2 ∈ A1 }) =
P ({X1 ∈ A1 })P ({X2 ∈ A2 }) · · · P ({X2 ∈ A1 })

On a donc :
E[X1 X2 · · · Xn ] = E[X1 ]E[X2 ] · · · E[Xn ]

Et si chacune de ces lois a une densité de probabilité f i (xi ), alors le “n-uplet”

(X1 , X2 , · · · , Xn )(x1 , x2 , · · · , xn )

admet pour densité de probabilité le produit :

h(x1 , x2 , · · · , xn ) = f (x1 )f (x2 ) · · · f (xn )

A.14 Somme de variables aléatoires indépendantes

A.14.1 Lois “faible” des grands nombres.

Soit X1 , X2 , ..., XN , ... une suite de variables aléatoires réelles deux à deux indépendantes
et ayant toutes la même moyenne E[X i ] = m et le même écart type fini : σ = E[X i2 ]. Alors
la moyenne arithmétique :
X1 + X 2 + · · · + X N
N
108 Rappels de théorie des probabilités

converge “en probabilité” vers la moyenne stochastique m lorsque n → ∞.

C’est à dire, quel que soit le nombre positif  donné (si petit soit-il) :
  
 X1 + X 2 + · · · + X N 
limN →∞ P   − m >  = 0

N

Démonstration :

Considérons la variable aléatoire :

ZN = X 1 + X 2 + · · · + X N

On a :  
ZN
E[ZN ] = N m, ou E =m
N
On peut se ramener à une somme de variables aléatoires centrées :

ZN − E[ZN ] = (X1 − m) + (X2 − m) + · · · + (XN − m)

En utilisant l’indépendance deux à deux des X i , i.e

E[Xi Xj ] = E[Xi ]E[Xj ] ∀i 6= j

On peut écrire :
N
X
E[(ZN − E[ZN ])2 ] = E[(Xi − m)2 ] = N σ 2
i=1

Soit, en divisant par N2 : " 2 #

ZN σ2
E −m = (A.19)
N N

En appliquant à la variable aléatoire (Z N /N − m) l’inégalité de Bienaymé-Tchebichef,

démontrée au paragraphe précédent :

σ2
  
 ZN 
P   − m >  <
 (A.20)
N N 2

Quel que soit  fixé, si petit soit-il, le second membre tend vers 0 lorsque N → ∞ C.Q.F.D.

Remarque importante :
√
La relation (A.18) montre que l’ecart type sur la variable aléatoire Z N /N est N fois plus
faible que l’écart type σ sur la variable aléatoire X
A.14 Somme de variables aléatoires indépendantes 109

Plus généralement soit une suite de variables aléatoires X 1 , X2 , ..., XN , ... deux à deux
indépendantes, à valeurs dans R d , ayant toutes la même loi de probabilité qu’une variable
aléatoire donnée X. Pour toute fonction g de R d dans R, la somme

SN g(X1 ) + g(X2 ) + · · · + g(XN )

=
N N
converge “en probabilité” vers la moyenne stochastique
√ ḡ = E[g(X)] lorsque N → ∞
et l’écart
p type sur la variable aléatoire S N /N est N fois plus faible que l’ecart type :
σg = E[(g(X) − ḡ)2 ] sur la variable aléatoire Y = g(X)

Cette généralisation est évidente :

Considérons les variables aléatoires Y n = g(Xn ), si les Xn sont deux à deux indépendantes,
alors les Yn le sont aussi. Nous pouvons donc appliquer la loi forte des grands nombres à
la suite de variables aléatoires réelles Y n et le théorème énoncé ci-dessus est démontré.
Le type de convergence (convergence en probabilité) que nous venons de démontrer, à partir
d’une inégalité triviale (de Bienaymé-Tchebichef) est “faible”.
On peut prouver une loi beaucoup plus forte, dite “loi forte des grands nombres” que nous
ne redémontrerons pas dans ce cours :
Pour toute suite de variables aléatoires indépendantes équidistribuées, la somme :

(X1 + X2 + · · · + Xn )/N

converge presque sûrement vers la moyenne m = E[X 1 ].

La convergence presque sûre de la variable aléatoire S n vers S a lieu si et seulement si
Si = φi (ω) tend vers S = φ(ω) pour tout résultat ω ∈ Ω, sauf peut-être pour certains
résultat ω particuliers appartenant à un ensemble de mesure nulle (i.e. de probabilité nulle).
[En théorie générale de la mesure on emploie le terme “convergence presque partout”]

A.14.2 Théorème de la limite centrale

La “loi des grands nombres” démontrée au paragraphe précédent nous dit que la vari-
able aléatoire ZN /N , somme de N variables aléatoires indépendantes tend vers la moyenne
stochastique E[ZN /N ]. Par contre, l’inégalité triviale de Bienaymé-Tchebichef utilisée pour
la démonstration ne donne pas une estimation précise de “l’erreur” :  N = (ZN /N −
E[ZN /N ]). Nous pouvons le faire dans le cas où les variables aléatoires X n sont indépendantes
dans leur ensemble (la loi de probabilité de (X 1 , X2 , · · · , XN ) est le produit des lois de prob-
abilités de X1 , X2 ,..., XN ).

Soit une suite de variables aléatoire indépendantes dans leur ensemble :

X1 , X 2 , · · · , X N
110 Rappels de théorie des probabilités

ayant toutes la même moyenne m et un même écart type σ fini.

Considérons la variable aléatoire :
√  
N ZN ZN (X1 − m) + (X2 − m) + · · · + (XN − m)
N = −E = √
σ N N σ N
Cette variable aléatoire a une moyenne nulle son écart type est l’unité. Elle est donc réduite.
√
Nous allons montrer que lorsque N → ∞ la loi de probabilité de σN N converge “en loi”
vers la “loi normale” (loi Gaussienne) réduite ayant pour densité de probabilité :
1 2
f (x) = √ e−x /2
2π
√
Nous pouvons réécrire σN N sous la forme :
√
N X1 − m X2 − m XN − m
SN = N = √ + √ + ··· + √ = Y1 + Y2 + · · · + Y N
σ σ N σ N σ N
Les variables :
Xi − m
Yi = √
σ N
√
sont indépendantes, de moyenne nulle et de même écart type 1/ N . A l’ordre deux leur
fonctions caractéristiques ont le même développement en série de Taylor :
t2 t2 1 √
FYi (t) = E[exp(itYi )] = 1 + E(Yi )it − E[Yi2 ] + ... = 1 − + O(1/ N )
2 2N N
la fonction caractéristique de la somme sécrit :
FSN (t) = E[exp it(Y1 + Y2 + ... + YN )] = E[exp(itY1 ) exp(itY2 ) · · · exp(itYN )]
et puisque les varaibles aléatoires exp(itY i ) sont indépendantes, l’espérance du produit est
le produit des espérances :
FSN (t) = E[exp(itY1 )]E[exp(itY2 )] · · · E[exp(itYN )] =
N
t2 √

1
1− + O(1/ N )
2N N
En prenant le logarithme :
t2 √ t2 √
 
1
LnFSN (t) = N Ln 1 − + O(1/ N) = − + O(1/ N )
2N N 2
Donc LnFSN (t) tend vers −t2 /2 lorsque N → ∞ et :
t2
lim FSN (t) = exp(− )
t→∞ 2
2
Mais exp(− t2 )
est la fonction caractéristique de la loi normale réduite. On en déduit que
SN converge en loi vers la loi normale réduite.
La convergence “en loi” est équivalente à la convergence des fonctions caractéristiques.
A.15 Divers modes de convergence 111

A.15 Divers modes de convergence

A travers les divers exemples que nous avons rencontrés, nous avons vu qu’en théorie des
probabilités, plusieurs types de convergence, plus ou moins fortes, pour des suites X n de
variables aléatoires peuvent avoir lieu :
– Nous venons de voir la convergence “en loi” qui a lieu si les fonctions caractéristiques
convergent simplement pour tout réel t.
– Nous avons vu lors de la démonstration de la loi faible des grands nombres “la convergence
en probabilités” :
P ({|Xn − X| > }) → 0
– Nous avons cité lors de l’énoncé de la loi forte des grands nombres la convergence “presque
sûre” (ou “presque partout”) :
Xn (ω) → X(ω)
pour tout événement ω de Ω sauf peut-être pour ω appartenant à un sous ensemble de
Ω de probabilité (mesure) nulle
– Citons encore la convergence en moyenne d’ordre p qui correspond à
Z
lim |Xn − X|p dP = 0
n→∞

En particulier pour p = 2 la convergence en moyenne quadratique qui correspond à :

lim E[(Xn − X)2 ] = 0

n→∞

La convergence presque sûre entraı̂ne la convergence en probabilités, qui entraı̂ne la conver-

gence en loi. La convergence en moyenne d’ordre p entraine la convergence en probabilités.
112 Rappels de théorie des probabilités
 B Chaı̂nes de Markov

B.1 Définition

Une suite infinie ordonnée (X1 , X2 , · · · , Xt , · · · ) de variables aléatoires constitue une chaı̂ne
de Markov si la loi de probabilité conditionnelle de X t+1 lorsqu’on se donne les valeurs de
X1 , X2 , ...,Xt , se réduit à la loi de probabilité conditionnelle de X t+1 lorsqu’on se donne
seulement la valeur de Xt .

Remarque importante : pour une suite ainsi définie, la connaissance du passé influe sur
l’avenir ; mais, d’un ensemble d’événements passés, seule subsiste l’influence du plus récent.
Si l’on se donne Xt , la connaissance de Xt0 pour t0 > t n’est en rien précisée par la donnée
des valeurs de Xt−1 , Xt−2 , ...,X1 .
Ici t designe un entier, nous avons choisi cette notation (plutôt que n) parce qu’en physique
cette variable représente souvent un temps discret. Nous verrons, en fin de paragraphe la
généralisation à une variable t continue.

Nous nous restreindrons aux Chaı̂nes de Markov “homogènes” où la loi de probabilité
conditionnelle de Xt+1 lorsqu’on se donne la valeur de Xt ne dépend pas de t. Alors la
chaı̂ne de Markov est entiérement déterminée par cette loi de probabilité conditionnelle
appelée “loi de transition” et par la loi de probabilité initiale de la variable aléatoire X 1 .

Nous supposerons aussi que la loi de probabilité conditionnelle de X t+1 lorsqu’on se donne
la valeur de Xt possède une densité de probabilité que nous noterons p(y/x) ou p(x → y),
c’est à dire :

P ({y < Xt+1 < y + dy}/{Xt = x}) = p(y/x)dy = p(x → y)dy (B.1)

Nous nous contenterons d’énoncer quelques résultats sur les chaı̂nes de Markov qui sont
indispensables à la compréhension de l’algorithme de Métropolis de la Méthode de Monte-
Carlo. Nous démontrerons ces résultats dans la cas discret et nous admettrons leur généralisation
pour des variables continues.
114 Chaı̂nes de Markov

B.2 Chaı̂nes dans un ensemble discret

B.2.1 Propriétés

Nous nous restreignons ici au cas où les variables aléatoires X t ne peuvent prendre que r
valeurs discrètes :
Xt (ω) ∈ {a1 , a2 , · · · , ar }
A “l’instant” t (discret, lui-aussi), la loi de probabilité de la variable aléatoire X t est définie
par un vecteur ligne π t dont le k ieme élément est P ({Xt = ak })

π t = (P ({Xt = a1 }), P ({Xt = a2 }), P ({Xt = a3 }), ....., P ({Xt = ak }), ...)

Puisque π t définit une loi de probabilité, la somme de ses composantes est 1 :

X
(π t )k = 1
k

A la loi de probabilité conditionnelle (probabilité de transition) P (X t+1 /Xt ) on peut as-

socier une matrice : p définie par :

(p)i,j = P ({Xt+1 = aj }/{Xt = ai })

Ou chaque élément (p)i,j représente la “probabilité de transition de l’état i à létat j”.

Cette matrice est “non négative” (tous ses élements sont positifs ou nuls) et elle satisfait
la condition : X
(p)i,j = 1 ∀i
j

C’est à dire que la somme de chaque ligne est 1. On peut traduire cette propriété par la
relation :
p.1 = 1
où 1 est le vecteur colonne dont tous les éléments sont égaux à 1. Cette relation signifie
que le nombre λ = 1 est valeur propre de la matrice p avec pour vecteur propre 1.

Connaissant la loi de probabilité π t de la variable Xt à l’instant t, la loi de passage, donc

ici la matrice p permet de déterminer la loi de probabilité π t+1 de Xt+1 à l’instant t + 1 :
L’axiome des probabilités conditionnelles permet décrire :

P ({Xt+1 = ak } ∩ {Xt = aj }) = P ({Xt = aj })P ({Xt+1 = ak }/{Xt = aj })

En sommant chacun des deux membres sur tous les j, on obtient :

X X
P ({Xt+1 = ak } ∩ {Xt = aj }) = P ({Xt = aj })P ({Xt+1 = ak }/{Xt = aj })
j j
B.2 Chaı̂nes dans un ensemble discret 115

Puisque pour i 6= j les évémements X t = ai et Xt = aj sont incompatibles, le premier

membre peut s’écrire, en vertu de la propriété d’additivité A.1

P ({Xt+1 = ak } ∩ [{Xt = a1 } ∪ {Xt = a2 } ∪ · · · ∪ {Xt = an } ∪ · · · ])

soit encore :
P ({Xt+1 = ak } ∩ Ω) = P ({Xt+1 = ak })
On obtient donc
X
(π t+1 )k = P ({Xt+1 = ak }) = P ({Xt = aj })P ({Xt+1 = ak }/{Xt = aj }) (B.2)
j

soit X
(π t+1 )k = (π t )j (p)jk
j

c’est à dire au sens de la multiplication d’un vecteur ligne par une matrice :

π t+1 = π t p

On en déduit
π t+α = π t pα
C’est l’equation de Chapman-Kolmogorov

B.2.2 Rayon spectral de la matrice p

.
Nous savons que λ = 1 est valeur propre de la matrice de passage p. Nous allons montrer
que c’est la plus grande valeur propre en module de cette matrice, c’est à dire son “rayon
spectral”.
Si λ est valeur propre de la matrice p, alors il existe un vecteur colonne v tel que :

p.v = λv

et il existe un vecteur ligne u associé à cette même valeur propre :

λu = u.p

soit : X
λuj = ui .pij ∀j
i

En prenant la valeur absolue de chaque membre, on a :

X
|λ||uj | = | ui .pij | ∀j
i
116 Chaı̂nes de Markov

et puisque la valeur absolue d’une somme est inférieure ou égale à la somme des valeurs
absolues : X
|λ||uj | ≤ |ui |.|pij | ∀j
i
Si maintenat nous sommons chaque membre sur j, nous obtenons
X XX XX
|λ| |uj | ≤ |ui |.|pij | = |ui |.|pij |
j j i i j

Mais puisque la matrice p satisfait la condition :

X
(p)i,j = 1 ∀i
j

(cf paragraphe précédent), on a alors :

X X
|λ| |uj | ≤ ui
j i

donc
λ≤1
CQFD

B.2.3 Convergence vers une loi de distribution stationaire

Si pour un temps “suffisament long” la loi de distribution π t de la chaı̂ne de Markov tend

vers une distribution π L qui devient invariante, c’est à dire indépendante du temps, on dit
que la distribution π L est stationaire et on a :

π L = π L .p

C’est à dire que π L est un vecteur propre (vecteur ligne) assocé à la valeur propre λ = 1

B.2.4 Chaı̂ne de Markov ergodique

D’après les paragraphes précédents, il existe au moins une distribution invariante π L , mais
elle n’est pas forcément unique. Pour cela on doit imposer une condition supplémentaire à
la matrice p.

Matrice régulière

Une matrice est non régulière s’il existe une permutation des indices qui permette de la
mettre sous la forme :  
A1 B
0 A2
B.2 Chaı̂nes dans un ensemble discret 117

dans ce cas :  
A1 B
(0, U ) = (0, U A2 )
0 A2
et en partant d’un état correspondant à la partie droite U du vecteur ligne (0, U ), on reste
dans ce même sous espace, sans jamais atteindre le sous espace correspondant à la partie
gauche.
Une matrice est régulière si elle n’est pas décomposable sous la forme ci-dessus.

Ergodicité

Si la matrice de passage p d’une chaı̂ne de Markov est régulière alors quel que soit le couple
d’états {ai , aj } il existe un entier n tel que la probabilité de passage de a i à aj au bout d’un
temps n soit strictement positive (i.e non nulle). Dans ce cas, tous les états sont visités au
cours du temps, ce qui correspond à la propriété physique “d’ergodicité”.
Le théorème mathématique de Perron-Fröbenius nous dit que pour une matrice non-
négative, régulière la valeur propre correspondant à son rayon spectral (i.e la valeur propre
de plus grand module) est non dégénérée, c’est à dire que le vecteur propre correspondant
est unique.
Pour une chaı̂ne de Markov ergodique, la distribution stationaire π L (i.e vecteur propre
correspondant à la valeur propre maximale λ = 1) est unique.

B.2.5 Généralisation à un ensemble d’états continus

Au lieu de raisonner sur une matrice de passage p à deux dimensions, on considère la

densité de probabilité p(x → y) de la loi de transition qui est une fonction de deux variables
continues, et dans les relations précédentes, démontrées dans le cas discret, on remplace les
sommes discrètes par des intégrales.
On peut, en particulier faire un passage à l’intégrale pour l’équation de Chapman-Kolmogorov
B.2 démontée pour un ensemble discret.
Dans le cas continu, si π t (x) représente la densité de probabilité correspondant à la loi de
probabilité de la variable Xt , au temps t + 1, la densité de probabilité associée à la loi de
probabilité de la variable aléatoire X t+1 est :
Z
t+1
π (x) = π t (u)p(u → x)du (B.3)

Cette relation représente l’équation de Chapman-Kolmogorov, pour un ensemble d’états

continu.
118 Chaı̂nes de Markov
 Rappels de Physique
C statistique

Les rappels de ce chapitre sont tirés de l’excellent livre de Diu et al. [20]

C.1 Entropie statistique

C.1.1 Définition de Shannon

Soit X une variable aléatoire discrète, sur un ensemble fini, avec pour loi de probabilité :

p(xi ) = P ({X = xi })

l’entropie est définie par :

M
X
S[X] = −E[log P (X)] = − p(xi ) log[p(xi )]
i=1

avec la convention : 0 log(0) = 0. En théorie de l’information, on utilise en général le

logarithme en base 2.
Elle mesure le “manque d’information” et possède les propriétés suivantes :
1. Elle est positive, complètement symétrique
2. Les événements de probabilité nulle ne contribuent pas
3. Elle atteint son minimum S = 0 si l’un des événements est réalisé avec une probabilité
1.
4. Elle atteint son maximum lorsque tous les événements sont èquiprobables
5. Additivité : considérons deux variables aléatoire indépendantes X, Y et la loi de
probabilité du couple :
pxm ,yn = P (X = xm , Y = yn )
Alors l’entropie du couple est :
X
S=− pxm ,yn log[pxm ,yn ]
m,n
120 Rappels de Physique statistique

Si X et Y sont indépendantes, alors p xm ,yn = p(xn )p(yn ) et il est facile de montrer

que : X X
S = −{ p(xm ) log[p(xm )] + p(yn ) log[p(yn )]}
m n

c’est à dire :
S[(X, Y )] = S[X] + S[Y ]
Le manque d’information sur le couple est la somme des manques d’information sur
chaque variable.
Plus généralement, si X et Y ne sont pas indépendantes, on peut montrer la propriété

S[(X, Y )] = S[X] + S[{Y /X}]

en faisant intervenir la probabilité conditionnelle P ({Y /X})

C.1.2 Théorème de Khinchin (1957).

Khinchin a montré que la seule fonction qui satisfasse aux propriétés 1-5 est :
M
X
S = −k p(xi ) ln[p(xi )]
i=1

On a le choix de la constante k : en thermodynamique, k est la constante de Boltzman et

le logarithme est en base e (logarithme Népérien)

C.2 Systèmes macroscopiques à l’équilibre

C.2.1 Ensemble micro-Canonique

On considère un système isolé.

Son état est déterminé par un ensemble de paramètres extérieurs :
– Son énergie E définie à δE près
– Son volume V
– Son nombre de particules N
– etc ...
Le système est caractérisé par ses états microscopiques. A chaque état l correspond une
probabilité Pl . Pour un système isolé à l’équilibre, puisqu’on ne connait aucune information
plus précise sur ce système, on postule que tous les états microscopiques accessibles sont
équiprobables (—entropie maximale—) :

Pl = 1/Ω si E < El < E + δE (C.1)

Pl = 0 sinon
C.2 Systèmes macroscopiques à l’équilibre 121

Ω représente le nombre d’états du système et on a :

X
Pl = 1
l

Son entropie est donc maximale :

S µC (E, V, N, ...) = k ln(Ω) (C.2)

Cette relation est l’effigie inscrite sur la tombe de Boltzman Cette hypothèse d’équipartition
correspond donc au manque d’information maximal sur le système !
On peut considérer diverses dérivées par rapport aux paramètres extérieurs :
– La température microcanonique T µC :
1
= ∂S µC (E, N, V, ...)/∂E (C.3)
T µC
– Le potential chimique microcanonique :
µµC
= −∂S µC (E, N, V, ...)/∂N (C.4)
T µc
– La pression microcanonique
P µC
= ∂S µC (E, N, V, ...)/∂V (C.5)
T µC

C.2.2 Ensemble Canonique

On considère un système S en contact avec un système beaucoup plus grand : réservoir de

chaleur ou thermostat T .
L’union de ces deux systèmes constitue un système isolé et son énergie totale :

ES + ET = Etot

est fixée à δE près

T étant beaucoup plus grand que S on a toujours E S  ET , et on peut considérer la
température microcanonique du thermostat T TµC comme indépendante de l’énergie du petit
système S et donc constante.
La température canonique T c du petit système S sera par définition la température micro-
canonique du thermostat.

T c = TTµC
Le systm̀e SU T étant isolé, on peut lui appliquer les postulats de l’ensemble microcanon-
ique. L’union des deux systèmes est caractérisée par l’ensemble des couples états (l, L)
d’energie
El,L = El + EL
122 Rappels de Physique statistique

tels que :
Etot ≤ El,L ≤ Etot + δE
Tous les états microcanoniques sont èquiprobable, et la probabilité P lC de trouver le système
dans létat d’energie l est égale au rapport du nombre détats Ω T (ET = Etotal − El ) du
thermostat
P d’energie ET = Etot − El sur le nombre total d’états de SU T : Ω SU T (Etot ) =
Ω
l T (E T = Etotal − El )
X
PlC = ΩT (Etot − El )/ ΩT (Etot − El )
i

ΩT (Etot − El ) s’exprime en fonction de l’entropie microcanonique du Thermostat :

ΩT (Etot − El ) = exp[STµC (Etot − El )/k]

Mais l’energie El est toujours faible devant l’énergie du grand système T , et un développement
au premier ordre donne :
∂S
STµC (Etotal − El ) = STµC (Etotal ) − El (Etot )
∂E
soit
1
STµC (Etotal − El ) = STµC (Etotal ) − El
kT
Ce qui conduit à
PlC = exp[−El /kT ]/Z (C.6)
avec X
Z(T, N, V, ...) = exp[−El /kT ] (C.7)
i

On définit le potentiel thermodynamique, ou énergie libre :

F = −kT ln Z(T, N, V, ...) (C.8)

– La valeur moyenne de l’énergie est :

X ∂
Ē = El exp[−El /kT ]/Z = − ln Z(T, N, V, ...) (C.9)
∂β
i

avec β = 1/T
– La chaleur spécifique à volume constant est :

∂ Ē
Cv = (C.10)
∂T
– L’entropie Canonique est :
X
S C (T, V, N, ...) = −k Pl ln[Pl ] = −k(ln Z − Ē/kT ) (C.11)
i
C.2 Systèmes macroscopiques à l’équilibre 123

Soit
X exp(−El /kT  El

c
S (T, V, N, ...) = −k − − ln(Z)
Z kT
l
døù
T S c (T, V, N, ...) = Ē + kT ln(Z)
ce qui permet d’écrire l’énergie libre sous la forme habituelle :

F = Ē − T S C (C.12)

– Le potentiel chimique canonique est :

∂F
µC = (C.13)
∂N
– La pression canonique est
∂F
PC = (C.14)
∂V

C.2.3 Ensemble Grand Canonique

On considère un petit système S en contact avec un réservoir R de chaleur et de particules.

L’union des deux systèmes S et R constitue un système isolé auquel on peut appliquer les
principes microcanoniques.
On procède comme pour l’ensemble Canonique, mais ici un état du système S est caractérisé
par son energie El et son nombre de particules Nl .
La probabilité de trouver le système S dans un état d’énergie E l , avec un nombre de
particules Nl est proportionnelle au nombre d’état du réservoir d’énergie E R = Etot − El
et de nombre de particules NR = Ntot − Nl soit :

PlGC = c Ω(ER = Etot − El , NR = Ntot − Nl )

µC
PlGC = c Ω exp[SR (ER = Etot − El , NR = Ntot − Nl )/k]
µC
Et en développant SR au premier ordre en El et Nl :

µC µC ∂S ∂S
SR = SR (ER , NR ) − El − Nl
∂E ∂N
µC µC 1 µ
SR = SR (ER , NR ) − El + Nl
T T
On deduit :
El − µNl
PlGC = exp[− ]/Z (C.15)
kT
avec :
X El − µNl
Z= exp[− ] (C.16)
kT
l
124 Rappels de Physique statistique

On définit le Grand Potentiel thermodynamique :

G = −kT ln Z (C.17)

En regroupant les termes correspondant à un même nombre de particules N l = n, on peut

écrire :
N {Nl =n}
X X
Z= [exp(µ/kT )]n exp(−El /kT )
n=1 l

où la seconde somme porte sur le sous ensemble d’etats d’energie E l correspondant à um
mm̂e nombre de particules Nl = n. Ce n’est autre que la fonction de partition canonique
Zn d’un système à n particules.
N
X
Z= [exp(µ/kT )]n Zn
n=1

La valeur moyenne du nombre de particules est

X El − µNl
N̄ = Nl exp[− ]/Z
kT
l

soit
∂ ln(Z)
N̄ = kT = −∂G/∂µ
∂µ

C.3 Evolution vers l’équilibre

C.3.1 Probabilité de transition

On postule que la probabilité pour qu’un système qui se trouve à l’état l à l’instant t 0 se
retouve à l’état m à l’instant t0 + dt (dt long devant le temps caractéristique des transitions
microscopiques, mais court à l’échelle macroscopiques est :
– Proportionnel à dt
– indépendant de t0
Pml (t0 , t0 + dt) = aml dt (C.18)

C.3.2 Equation maı̂tresse

A un instant t, état macroscopique d’un système est caractérisé par l’ensemble des proba-
bilités Pl (t) pour qu’il se trouve dans chacun de ses états microscopiques l.
Soit dPl l’évolution de Pl entre les instants t et t + dt. On peut écrire intuitivement :
C.3 Evolution vers l’équilibre 125

X X
dPl = alm Pm (t)dt − aml Pl (t)dt (C.19)
m m

On postule “l’équation maı̂tresse suivante :

dPl X
= [alm Pm (t) − am,l Pl (t)] (C.20)
dt m

C.3.3 Système isolé

– Les alm ne relient que des états de même energie :

aµC
lm = 0 si |El − Em | > δE

– Les alm sont symétriques

aµC µC
lm = aml

C.3.4 Système S en contact avec un thermostat T

On considère l’union SUT comme système isolé.

La probabilité Pl,L (t) pour que S soit dans l’état l et T dans létat L obéit à l’équation
maı̂tresse :
dPl,L X µC
= [a(l,L)(m,M ) Pm,M (t) − aµC
(m,M )(l,L) Pl,L (t)]
dt
m,M

La probabilité Pl pour que le système S soit dans l’état l est

X
Pl = Pl,L
L

On a
Pl,L = Pl PLµC (Etot − El )
En substituant dans la relation :
dPl X dPl,L
=
dt dt
L

on obtient l’équation maı̂tresse :

dPl X
== [aC C
lm Pm (t) − am,l Pl (t)] (C.21)
dt m

avec
aµC µC
X
aC
lm = (l,L)(m,M ) PM (Etot − Em )
M,L
126 Rappels de Physique statistique

aµC µC
X
aC
ml = (m,M )(l,L) PM (Etot − El )
M,L
on a
µC 1
PM (Etot − El ) = = exp[−S µC (Etot − El )]
Ω(Etot − El )
et en développant S µC au premier ordre :
∂S µC
S µC (Etot − El ) = S µC (Etot − El
∂E
on en déduit :
aC C
lm exp(−Em /kT ) = aml exp(−El /kT ) (C.22)

C.3.5 Propriétés à l’équilibre. Relation du bilan détaillé.

A l’équilibre :

dPl /dt = 0 ∀l
L’équation maı̂tresse sécrit alors :
X
e
[alm Pm − aml Ple ] = 0
m

L’équation ci dessus, avec les contraintes :

Ple ≥ 0 ∀l
et X
Ple = 1
l
admet pour seule solution :
– La distribution microcanonique pour un système isolé
– La distribution canonique pour un système en contact avec un thermostat

1. Pour un système isolé :

P µC (E) = 1/Ω(E) E ≤ El ≤ E + δE
µC
P (E) = 0 sinon
On a donc :
PlµC (E) = Pm
µC
(E)
Puisque aµC µC
lm = aml , on a une relation terme à terme beaucoup plus forte, dite rela-
tion du bilan détaillé :
[alm Pm − aml Pl ] = 0 ∀l, ∀m (C.23)
On montre que c’est la seule solution d’équilibre.
C.3 Evolution vers l’équilibre 127

2. Pour un système en contact avec un thermostat :

PlC = exp(−El /kT )

et
aC C
lm exp(−Em /kT ) = aml exp(−El /kT )

Et on a aussi la relation du bilan détaillé :

[alm Pm − aml Pl ] = 0 ∀l, ∀m

128 Rappels de Physique statistique
 Valeurs moyennes
D d’observables

D.1 Définition

Système quantique

Un système quantique est décrit par son Hamiltonien H. Ses différents états sont les vecteurs
propres |n > de H, avec pour energies E n les valeurs propres correspondantes :

H|n >= En |n >

Une grandeur physique, ou “observable” est caractérisée par un opérateur O et la mesure

de cette grandeur, lorque le système est dans un état |n > est

< n|O|n >

Dans l’ensemble Canonique, la fonction de partition est, d’après les rappels précédents : :
X
Z= exp(−βEn )
n

avec β = 1/kT .
On peut écrire en considérant une base d’états propres |n >
X
Z= < n| exp(−βH)|n >
n

Z apparaı̂t donc comme la trace de l’opérateur exp(−βH) et on sait que la trace d’un
opérateur ne dépend pas de la base choisie. Donc on ècrit de manière plus générale :

Z = T r[exp(−βH)] (D.1)

La valeur moyenne ou “espérance d’un opérateur O sécrit, dans l’ensemble canonique (cf.
rappels précédents) : : P
< n|O|n > exp(−βEn )
< O >= n
Z
130 Valeurs moyennes d’observables

soit P
n < n|O exp(−βEn )|n >
< O >=
Z
ou encore : P
< n|O exp(−βH)|n >
n
< O >=
Z
On voit apparaı̂tre la trace de l’opérateur O exp(−βH), indépendante du choix de la base.
On écrira :
T r[O exp(−βH)]
< O >= (D.2)
Z

Système classique

Pour un système classique, décrit par un Hamiltonien classique : H(p, x), la valeur moyenne
d’une grandeur physique, ou “observable” O (par exemple : énergie, aimantation pour un
système magnétique, etc...) est, dans l’Ensemble Canonique[21, 22, 23, 20] :

O(x)e−βH(p,x) dpdx
R
< O >= (D.3)
Z
où Z représente la fonction de partition :
Z
Z = e−βH(p,x) dpdx (D.4)

On note β = 1/kB T , x = (α1 , α2 , · · · , α3N ) représente les coordonnées des particules et

p = (p1 , p2 , · · · , p3N ) les variables conjuguées ou impulsions : p i = dαi /dt. L’Hamiltonien
H du système s’écrit :
X p2
i
H= + E(x) (D.5)
2mi
i

Le premier terme représente l’energie cinétique du système et le second terme E(x) l’energie
potentielle.

Dans ce cas, les intégrales sur les p i se séparent et se simplifient entre numérateur et
dénominateur et on peut écrire plus simplement :

O(x)e−βE(x) dx
R
< O >= (D.6)
Z
avec : Z
Z= e−βE(x) dx (D.7)

La fonction :
f (x) = e−βE(x) /Z (D.8)
D.1 Définition 131

R
est positive et a pour intégrale : f (x)dx = 1, elle peut donc représenter la densité de
probabilité d’une variable aléatoire. Elle est grande dans les portions de l’espace qui nous
intéressent pour la calcul de la valeur moyenne de toute observable O et décroı̂t rapide-
ment en dehors. Elle constitue donc le meilleur choix pour un “échantillonnage suivant
l’importance”.

Il subsiste toutefois un problème important. Dans la définition de f , nous ne connaissons pas

la constante de normalisation Z, c’est à dire la fonction de partition. Il est donc impossible
d’échantillonner directement f .

Nous devons donc avoir recours à un autre procédé d’échantillonnage. Nous utiliserons pour
outil des chaı̂nes de Markov.
132 Valeurs moyennes d’observables
 Bibliographie

[1] D.E. Knuth, Seminumerical Algorithms, 2nd. ed. vol. 2 of The Art of Computer Pro-
gramming (Reading, Mass. : Addison-Wesley), (1981).
[2] D.H. Lehmer, Proc. 2nd Symp, Large-Scale Digital Calculating machinery (Harvard
Univ. Press, Cambridge, MA, 1951) pages : 141-146
[3] W. H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes (Cam-
bridge University Press, New-York) (1986).
[4] G.J. Mitchell and D.P. Moore (1958) unpublished.
[5] Heringa, Blöte, and Compagner, Int. J. Mod. Phys. C3, 561 (1992) et références citées.
[6] R. C. Tauthworthe, Math. Comput. 19, 201 (1965).
[7] J.M. Hammersley et D.C Handscomb, Les Méthodes de Monte-Carlo (Ed. Dunod,
Paris), (1967).
[8] J.W Negele et H. Orland, Quantum Many-Particle System, Chapitre 8 : “Stochastic
methods”, (Ed. Addison-Wesley, New-York) (1988).
[9] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller et E. Teller, J. Chem. Phys.
21, 1087 (1953).
[10] L. Verlet, Phys. Rev. 159, 98 (1967).
[11] M.P. Allen et D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford
(1987).
[12] J.B. Kogut, J. of Stat. Phys. 43, 771 (1986).
[13] A.D. Kennedy, Nuclear Physics B (Proc. Suppl.) 4, 576 (1988).
[14] M. Rosenbluth, A Rosenbluth, J. Chem. Phys. 23, 356 (1954).
[15] D. Frenkel, G. Mooij, B. Smith, J. Phys. Cond. Mat. 3, 3053 (1991).
[16] R.P. Feynman and A.R. Hibbs, Quantum Mechanics and Path Integrals, Ed. Mc Graw-
Hill (1965).
[17] D.M. Ceperley, Rev. Mod. Phys. 67, 279 (1995).
Voir aussi les cours de Ceperley, accesssibles gratuitement sur internet :
http ://people.physics.uiuc.edu/Ceperley/papers/163.pdf
http ://www.phys.uri.edu/ nigh/QMC-NATO/webpage/abstracts/ceperley.ps
134 BIBLIOGRAPHIE

[18] J. Bass, Eléments de calcul des probabilités ; théorique et appliqué. Ed. Masson [Paris]
(1967).
[19] J. Neveu, Cours de probabilités : Ecole Polytechnique (1974) ; même auteur : Proba-
bilités, Edition 1995, Ecole Polytechnique [Palaiseau, France].
[20] B. Diu, C. Guthmann, D. Lederer, B. Roulet, Physique Statistique, Hermann Paris
(1989).
[21] R. Balian, Du microscopique au macroscopique, Cours de Physique statistique de
l’Ecole polytechnique, tomes 1 et 2. Ellipses (1982).
[22] L. Landau et E. Lifschitz, Physique Statistique, Editions Mir, Moscou (1967).
[23] F. Reif, Fundamentals of Statistical and Thermal Physics, McGraw-Hill, New York
(1965).