{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Température\n",
    "\n",
    "Marc Lorenzi\n",
    "\n",
    "19 avril 2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (8, 8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Un jeu de hasard"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 À quoi allons-nous jouer ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $S$ une liste de taille $N$. Soit $q\\in\\mathbb N$. On initialise tous les éléments de $S$ à la valeur $q$, puis on itère l'opération suivante :\n",
    "\n",
    "On choisit au hasard un entier $i$ entre 0 et $N-1$. Si $S[i]>0$,\n",
    "- On choisit au hasard $j$ entre 0 et $N-1$.\n",
    "- On augmente $S[j]$ de 1.\n",
    "- On diminue $S[i]$ de 1.\n",
    "\n",
    "On fait cela un \"grand\" nombre de fois puis on regarde l'état de la liste $S$. Combien de cases à 0 ? Combien à 1 ? etc.\n",
    "\n",
    "L'idée est de modéliser (très naïvement) un système composé de $N$ particules. Au départ, toutes ces particules ont une énergie $qE_0$ où $E_0>0$ est l'énergie d'un \"quantum d'énergie\". Chaque particule peut interagir avec une autre particule en perdant un quantum d'énergie et en donnant ce quantum à l'autre particule."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Jouons !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `iterer` prend en paramètres \n",
    "- Un entier $N$, le nombre d'éléments de notre liste.\n",
    "- Un entier $q$ pour initialiser la liste.\n",
    "- Un entier `niter`, le nombre d'itérations à effectuer.\n",
    "\n",
    "La fonction renvoie la liste $S$ après itérations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def iterer(N, q, niter):\n",
    "    S = [q for i in range(N)]\n",
    "    for k in range(niter):\n",
    "            i = random.randint(0, N - 1)\n",
    "            if S[i] != 0:\n",
    "                j = random.randint(0, N - 1)\n",
    "                S[j] += 1\n",
    "                S[i] -= 1\n",
    "    return S"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici un exemple."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 10000\n",
    "S = iterer(N, 2, 100000)\n",
    "plt.plot(S, 'k')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ce que nous voyons ci-dessus c'est la liste $S$ après un certain nombre d'itérations. C'est assez confus, alors intéressons-nous à __l'histogramme__ $S$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 L'histogramme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `histogramme` renvoie ... __l'histogramme__ de la liste $S$. En clair, elle renvoie une liste $h$ telle que $h[0]$ est la probabilité qu'un élément de $S$ ait la valeur 0, $h[1]$ est la probabilité qu'un élément de $S$ ait la valeur 1, etc. En fait de probabilités, je devrais parler de __fréquences__. ne voulant pas entrer dans les détails pour ne pas lasser le lecteur, je serai assez léger sur la rigueur du vocabulaire tout au long du notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def histogramme(S):\n",
    "    N = len(S)\n",
    "    M = max(S)\n",
    "    h = (M + 1) * [0]\n",
    "    for i in range(N):\n",
    "            h[S[i]] += 1 / N\n",
    "    return h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = iterer(10000, 3, 100000)\n",
    "h = histogramme(S)\n",
    "plt.plot(h, 'k')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On obtient pour l'histogramme une fonction décroissante. Est-ce toujours le cas ? Faites des tests, changez les valeurs des paramètres, vous verrez qu'effectivement, si `niter` est \"assez grand\" on retrouve des courbes similaires."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quelques statistiques ? La fonction `moment` prend un histogramme $h$ et un entier $k$ en paramètres. Elle renvoie le __moment d'ordre $k$__ de $h$, c'est à dire le réel\n",
    "\n",
    "$$\\mu(h,k)=\\sum_{i=0}^{n-1} i^k h_i$$\n",
    "\n",
    "où $n$ est la taille de la liste $h$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def moment(h, k):\n",
    "    s = 0\n",
    "    n = len(h)\n",
    "    for i in range(n):\n",
    "        s += i ** k * h[i]\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour $k=0$, cela doit faire 1 à tout prix. Pour $k=1$, on obtient la __moyenne__ (ou, plus pompeusement, __l'espérance mathématique__) de $h$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il est aussi facile de définir variance et écart-type de $h$ :\n",
    "\n",
    "$$V(h)=\\mu(h,2)-\\mu(h, 1)^2$$\n",
    "\n",
    "$$\\sigma(h)=\\sqrt{V(h)}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def moyenne(h): return moment(h, 1)\n",
    "\n",
    "def variance(h):\n",
    "    return moment(h, 2) - moyenne(h) ** 2\n",
    "\n",
    "def ecart_type(h):\n",
    "    return math.sqrt(variance(h))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici ce que cela donne pour notre jeu."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = iterer(10000, 2, 100000)\n",
    "h = histogramme(S)\n",
    "print(moyenne(h))\n",
    "print(ecart_type(h))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Une moyenne de 2 ne devrait pas vous surprendre ...\n",
    "\n",
    "Tentons une petite expérience. Traçons les histogrammes obtenus pour différentes valeurs de `niter` sur un même graphique. Les courbes claires correspondent à de petites valeurs de `niter`, les courbes sombres à de grandes valeurs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ks = [2 ** k for k in range(12, 20)]\n",
    "n = len(ks)\n",
    "plt.xlim(0, 10)\n",
    "plt.ylim(0, 0.4)\n",
    "for k in range(n):\n",
    "    S = iterer(10000, 3, ks[k])\n",
    "    h = histogramme(S)\n",
    "    c = 1. - k / n\n",
    "    plt.plot(h, color=(c, c, c))\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Les histogrammes ont l'air de \"tendre\" vers une sorte d'histogramme limite : le système $\\mathcal S$ que nous essayons de modéliser tend-il vers un état limite ? Mais qu'est-ce qu'un __état__ ? Nous allons voir qu'il existe deux sortes d'états : les micro-états et les macro-états."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Un peu de théorie"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Micro-états, macro-états"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1.1 Micro-états"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reprenons notre jeu favori. La liste $S$, de taille $N$, est un système possédant une énergie égale à $E=qNE_0$, où $E_0>0$ est l'énergie d'un \"quantum\". Cette énergie est répartie entre les différentes cases de la liste. Combien y a-t-il de configurations possibles de la liste $S$ ? Un nombre gigantesque. Vraiment énorme. Nous appellerons chacune de ces configurations un __micro-état__ du système. Une façon fun de compter les micro-états est la suivante :\n",
    "\n",
    "On dispose de $Nq+N$ objets, dont $N$ sont des bâtons et $Nq$ sont des étoiles. Les étoiles représentent les quanta d'énergie, les bâtons sont les bords gauches des cases de la liste $S$. Histoire de me faire bien comprendre, regardons l'exemple ci-dessous. Prenons $N=4$ et $q=3$. Le schéma\n",
    "\n",
    "$$|*|****|***|****$$\n",
    "\n",
    "représente la liste $S=[1,4,3,4]$.\n",
    "\n",
    "Combien y a-t-il de façon de placer les bâtons parmi les étoiles ? Le premier bâton se trouve évidemment toujours au début, c'est le bord gauche de la zéroième case de $S$. Il nous reste à placer $N-1$ bâtons, c'est à dire à choisir $N-1$ objets parmi les $Nq+N-1$ restants et les décréter \"bâtons\". Les autres sont des étoiles. Au total, nous avons donc\n",
    "\n",
    "$$\\Omega=\\binom{Nq+N-1}{N-1}=\\frac{(Nq+N-1)!}{(N-1)!(Nq)!}\\quad{\\rm micro-états}$$\n",
    "\n",
    "Faisons intervenir $E=NqE_0$ :\n",
    "\n",
    "$$\\Omega=\\frac{(E/E_0+N-1)!}{(N-1)!(E/E_0)!}=\\frac 1{(N-1)!}(\\frac E {E_0}+1)(\\frac E {E_0}+2)\\ldots(\\frac E {E_0}+N-1)$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histoire d'avoir une idée de l'ordre de grandeur du nombre de micro-états de notre jeu, voici une fonction `log_Omega` qui renvoie $\\log\\Omega$ (en base $10$). En effet, il est inutile de vouloir calculer $\\Omega$, c'est un nombre bien trop gigantesque pour Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_Omega(N, q):\n",
    "    s = 0\n",
    "    for k in range(1, N):\n",
    "        s += math.log(N * q + k, 10) - math.log(k, 10)\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un petit exemple ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "log_Omega(10000, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour $N=10000$ et $q=3$ nous avons environ $10^{9800}$ micro-états, un 1 avec 9800 zéros derrière ... c'est beaucoup, le nombre d'atomes dans l'univers est, en étant TRÈS généreux, $10^{100}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1.2 Macro-états"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lorsque nous observons l'histogramme $h$, nous ne retenons que certaines informations associées à la liste. Combien de cases de valeur nulle ? De valeur 1 ? etc. Deux listes peuvent avoir le même histogramme et être totalement différentes. Nous appellerons l'histogramme un __macro-état__. Un macro-état pour notre jeu est donc une suite $(h_j)_{j\\ge 0}$ de multiples entiers de $\\frac 1 N$ telle que \n",
    "\n",
    "$$\\sum_{j=0}^{\\infty}jh_j=q \\quad{\\rm et}\\quad \\sum_{j=0}^{\\infty}h_j=1$$\n",
    "\n",
    "Combien notre jeu possède-t-il de macro-états ? Moins que de micro-états, évidemment, puisque chaque macro-état correspond à plusieurs (beaucoup de) micro-états. Le dénombrement des macro-états de notre système est une question ouverte, si quelqu'un a une idée je suis preneur."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Hypothèses de travail"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $\\mathcal S$ un système. Nous supposons que le système peut être décrit par un très grand nombre de micro-états ayant tous la même probabilité. Nous pouvons par exemple penser à un système formé de particules ayant chacune une position, une vitesse, etc. Un micro-état serait la donnée de tous les paramètres associés à une particule, ceci pour toutes les particules.\n",
    "\n",
    "Ce que le physicien mesure, en revanche, ce sont des quantités issues des macro-états du système (la pression d'un gaz par exemple). Tous les macro-états n'ont pas la même probabilité, mais nous allons faire l'hypothèse suivante :\n",
    "\n",
    "__Hypothèse__ : Parmi tous les macro-états possibles, un système \"choisit\" celui qui correspond au plus grand nombre possible de micro-états.\n",
    "\n",
    "On peut montrer que cette hypothèse est une conséquence des trois hypothèses ci-dessous :\n",
    "\n",
    "1. Tous les micro-états du système ont la même probabilité.\n",
    "2. La dynamique interne du système est telle que les micro-états du système changent continuellement.\n",
    "3. Si on attend suffisamment longtemps (cela peut être très, très, très longtemps :-), le système explorera tous les micro-états possibles et passera un temps égal dans chacun d'entre-eux.\n",
    "\n",
    "Munis de cette magnifique hypothèse, nous allons définir le mot __température__."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Température"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considérons deux systèmes qui peuvent échanger de l'énergie entre-eux, mais avec rien d'autre. Le système $\\mathcal S_1$ a une énergie $E_1$ et le système $\\mathcal S_2$ a une énergie $E_2$. L'énergie totale \n",
    "\n",
    "$$E=E_1+E_2$$\n",
    "\n",
    "est donc, avec nos hypothèses, constante. Les deux systèmes $\\mathcal S_1$ et $\\mathcal S_2$ procèdent à des échanges d'énergie, et au bout d'un certain temps un équilibre s'établit : __l'équilibre thermique__. Cet équilibre est celui qui maximise le nombre de micro-états.\n",
    "\n",
    "\n",
    "Le système $\\mathcal S_i$ ($i=1,2$) peut être dans n'importe lequel de $\\Omega_i(E_i)$ micro-états. Un micro-état de notre système global est donc un couple formé d'un micro-état de $\\mathcal S_1$ et d'un micro-état de $\\mathcal S_2$. Le nombre total de micro-états du système est ainsi\n",
    "\n",
    "$$\\Omega_1(E_1)\\Omega_2(E_2)$$\n",
    "\n",
    "Nous cherchons donc à maximiser cette quantité. Nous allons travailler dans le cadre où $\\Omega_i$ est une fonction continue de $E_i$. Mais comment une fonction à valeurs entières peut-elle être continue ? Si cela vous chagrine, rappelez-vous le nombre de micro-états de notre jeu. Et travaillez sur la fonction $10^{-9800}\\Omega_i$ qui est un nombre entre 0 et 1, avec 9800 chiffres après la virgule. Du coup, l'hypothèse de quantités continues n'est plus du tout absurde :-). Dérivons par rapport à $E_1$. En un maximum cette dérivée doit être nulle.\n",
    "\n",
    "$$\\frac{d}{dE_1}(\\Omega_1(E_1)\\Omega_2(E_2))=0$$\n",
    "\n",
    "Développons :\n",
    "\n",
    "$$\\Omega_2(E_2)\\frac{d}{dE_1}\\Omega_1(E_1)+\\Omega_1(E_1)\\frac{d}{dE_2}\\Omega_2(E_2)\\frac{dE_2}{dE_1}=0$$\n",
    "\n",
    "Comme $E_1+E_2=E$ est constante, on a $\\frac{dE_2}{dE_1}=-1$. Ainsi,\n",
    "\n",
    "$$\\Omega_2(E_2)\\frac{d}{dE_1}\\Omega_1(E_1)=\\Omega_1(E_1)\\frac{d}{dE_2}\\Omega_2(E_2)$$\n",
    "\n",
    "ce que l'on peut encore écrire\n",
    "\n",
    "$$\\frac 1 {\\Omega_1(E_1)}\\frac{d}{dE_1}\\Omega_1(E_1)=\\frac 1{\\Omega_2(E_2)}\\frac{d}{dE_2}\\Omega_2(E_2)$$\n",
    "\n",
    "ou encore\n",
    "\n",
    "$$\\frac{d}{dE_1}\\ln\\Omega_1(E_1)=\\frac{d}{dE_2}\\ln\\Omega_2(E_2)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Posons\n",
    "\n",
    "$$\\frac 1 {k_BT_i}=\\frac {d\\ln\\Omega_i(E_i)} {dE_i}$$\n",
    "\n",
    "où $k_b=1.3807\\times 10^{-23}JK^{-1}$ est la __constante de Boltzmann__. Le nombre $T_i$ est la __température__ du système $\\mathcal S_i$, de sorte que l'équilibre thermique s'écrit\n",
    "\n",
    "$$T_1=T_2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ceci nous suggère la\n",
    "\n",
    "__Définition__ : La température $T$ d'un système $\\mathcal S$ est donnée par\n",
    "\n",
    "$$\\frac 1 {k_BT}=\\frac {d\\ln\\Omega(E)} {dE}$$\n",
    "\n",
    "où $\\Omega$ est le nombre de micro-états du système, vu comme une fonction de son énergie $E$.\n",
    "\n",
    "__Remarque__ : Si l'apparition du nombre magique $k_B$ vous inquiète, il vous suffit de remarquer que la valeur de la constante en question dépend des unités choisies. Historiquement, la différence entre la température d'ébullition de l'eau et celle de sa fusion a été fixée par un certain M. Celsius à 100. Pourquoi pas ? Puis un certain M. Kelvin a décidé que l'unité de température serait une simple translation de l'échelle Celsius. Si M. Celsius avait décidé de mettre $6.022\\times 10^{23}$ à la place de 100, $k_B$ aurait évidemment une tout autre valeur.\n",
    "\n",
    "__Exercice__ : quelle valeur ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Probabilités"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considérons à nouveau deux systèmes que nous appellerons $\\mathcal R$ et $\\mathcal S$. Nous faisons les deux hypothèses suivantes :\n",
    "\n",
    "- Le système $\\mathcal R$ est \"énorme\", nous l'appellerons le __réservoir__. On peut retirer ou ajouter à ce système de l'énergie sans que sa température soit modifiée de façon appréciable.\n",
    "- Pour chaque énergie possible $\\varepsilon$ du système $\\mathcal S$ il y a un unique micro-état de $\\mathcal S$ associé. Avec des notations évidentes, $\\Omega_{\\mathcal S}(\\epsilon)=1$.\n",
    "\n",
    "Imaginez par exemple que $\\mathcal S$ est une casserole d'eau, plongée dans l'océan $\\mathcal R$. Comme envisagé plus haut, les deux systèmes peuvent échanger de l'énergie et nous supposons que l'énergie totale du système, $E$, est constante. À l'équilibre thermique les deux systèmes ont la même température qui est, avec nos hypothèses, celle du réservoir (de combien l'océan se réchauffe-t-il quand vous y versez une casserole d'eau chaude ?)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La probabilité $P(\\epsilon)$ que le système $\\mathcal S$ ait l'énergie $\\varepsilon$ est proportionnelle au nombre total de micro-états, qui est $\\Omega(E)=\\Omega_{\\mathcal R}(E-\\epsilon)\\times \\Omega_{\\mathcal S}(\\epsilon)$ :\n",
    "\n",
    "$$P(\\epsilon)=\\lambda\\Omega_{\\mathcal R}(E-\\epsilon)$$\n",
    "\n",
    "où $\\lambda>0$ est choisi pour que $P$ soit une probabilité (la somme sur toutes les valeurs possibles de $\\epsilon$ doit valoir 1). On a, en faisant un développement limité à l'ordre 1,\n",
    "\n",
    "$$\\ln \\Omega_{\\mathcal R}(E-\\epsilon)=\\ln\\Omega_{\\mathcal R}(E)-\\epsilon\\frac{d\\ln\\Omega_{\\mathcal R}(E)}{dE}+o(\\epsilon)$$\n",
    "\n",
    "ou encore\n",
    "\n",
    "$$\\ln\\Omega_{\\mathcal R}(E-\\epsilon)=\\ln\\Omega_{\\mathcal R}(E)-\\frac{\\epsilon}{k_BT}+o(\\epsilon)$$\n",
    "\n",
    "d'où\n",
    "\n",
    "$$\\Omega_{\\mathcal R}(E-\\epsilon)=\\Omega_{\\mathcal R}(E)e^{-\\epsilon/k_BT}(1+o(\\epsilon))$$\n",
    "\n",
    "Admettons que dans le cadre où nous nous sommes placés on peut négliger le terme $o(\\varepsilon)$. On obtient donc\n",
    "\n",
    "$$P(\\epsilon)=\\lambda e^{-\\epsilon/k_BT}$$\n",
    "\n",
    "Pour chaque valeur de $T$ (fixée par le réservoir), nous avons une distribution de probabilité sur l'ensemble des énergies du système $\\mathcal S$. Elle est appelée la __distribution de Boltzmann__."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Après ce long intermède théorique, voyons ce que cela donne pour notre jeu !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Retour au jeu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quelle est la valeur de la case $S[0]$ ? Eh bien c'est un entier naturel ... Mais encore ? Quelle est, disons, la probabilité que $S[0]=k$, où $k\\in\\mathbb N$ ? Regardons notre jeu comme la réunion de deux systèmes $\\mathcal R$ et $\\mathcal S$, où $\\mathcal S$ est le système formé par la case $S[0]$ et $\\mathcal R$, le réservoir, est la liste formée de toutes les autres cases. Nous pouvons maintenant appliquer ce que nous avons raconté dans la section précédente. Il faudrait en toute rigueur remplacer $N$ par $N-1$ dans les calculs qui vont suivre (le réservoir n'a que $N-1$ cases), mais je fais confiance au lecteur pour vérifier que cela ne change pas les résultats obtenus : nous allons de toute façon faire tendre $N$ vers l'infini."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Estimations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Souvenons-nous que dans notre jeu le nombre de micro-états est\n",
    "\n",
    "$$\\Omega=\\frac 1{(N-1)!}(\\frac E {E_0}+1)(\\frac E {E_0}+2)\\ldots(\\frac E {E_0}+N-1)$$ \n",
    "\n",
    "On a donc\n",
    "\n",
    "$$\\ln\\Omega=\\sum_{k=1}^{N-1}\\ln(\\frac E {E_0}+k)-\\ln (N-1)!$$\n",
    "\n",
    "De là,\n",
    "\n",
    "$$\\frac 1 {k_BT}=\\frac{d\\ln\\Omega}{dE}=\\frac 1 {E_0}\\sum_{k=1}^{N-1}\\frac 1{\\frac E {E_0}+k}=\\frac 1 {E_0}\\sum_{k=1}^{N-1}\\frac 1{Nq+k}$$\n",
    "\n",
    "Réécrivons tout cela sous la forme\n",
    "\n",
    "$$\\frac {E_0} {k_BT}=\\frac 1 {Nq}\\sum_{k=1}^{N-1}\\frac 1{1+\\frac k{Nq}}$$\n",
    "\n",
    "Nos yeux experts voient apparaître une somme de Riemann associée à l'intégrale\n",
    "\n",
    "$$\\int_0^{\\frac 1 q}\\frac{dx}{1+x}=\\ln(1+\\frac 1 q)$$\n",
    "\n",
    "Ainsi, en passant à la limite lorsque $N\\to\\infty$,\n",
    "\n",
    "$$P(\\epsilon)=\\frac \\lambda {E_0}e^{-\\epsilon\\ln(1+1/ q)/E_0}=\\frac \\lambda {E_0} \\frac{1}{(1+\\frac 1 q)^{\\epsilon/E_0}}$$\n",
    "\n",
    "Que vaut $\\lambda$ ? Eh bien les énergies possibles sont les multiples entiers de $E_0$. Nous avons donc $\\sum_{i=0}^\\infty P(iE_0)=1$. Or,\n",
    "\n",
    "$$\\sum_{i=0}^\\infty\\frac{1}{(1+\\frac 1 q)^i}=\\frac 1 {1-\\frac 1 {1+\\frac 1 q}}=q+1$$\n",
    "\n",
    "Ainsi, $\\lambda =\\frac {E_0}{q+1}$ :\n",
    "\n",
    "$$P(\\epsilon)=\\frac 1 {q+1} \\frac{1}{(1+\\frac 1 q)^{\\epsilon/E_0}}$$\n",
    "\n",
    "ou encore, pour s'en tenir aux valeurs possibles de l'énergie :\n",
    "\n",
    "$$P(kE_0)=\\frac 1 {q+1} \\frac{1}{(1+\\frac 1 q)^{k}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `proba` ci-dessous prend deux paramètres $k$ et $q$. Elle renvoie $P(k E_0)$, le paramètre $q$ étant le paramètre fixé pour notre système.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def proba(k, q):\n",
    "    return 1 / (1 + 1 / q) ** k * 1 / (q + 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voilà donc l'histogramme \"limite\" dont nous avons parlé plus haut, que nous obtiendrions après une infinité d'itérations de notre jeu."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ks = [k / 100 for k in range(2000)]\n",
    "ys = [proba(k, 3) for k in ks]\n",
    "plt.plot(ks, ys, 'k')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sur cet exemple, une case de la liste a une probabilité d'environ $\\frac 1 4$ de contenir la valeur 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Intermède rigolo - Quelle est la température de notre jeu ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quelle est la température de notre jeu ? Les calculs précédents montrent que\n",
    "\n",
    "$$T=\\frac{E_0}{k_B\\ln(1+\\frac 1 q)}$$\n",
    "\n",
    "Il va nous falloir prendre une décision quant à la valeur de $E_0$. Quelle est la masse d'une molécule d'azote ? La masse molaire du diazote est environ 28 grammes par mole."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_Avo = 6.022e23\n",
    "m_azote = 1 / N_Avo * 28 * 1e-3\n",
    "print(m_azote)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sachant que la vitesse d'une telle molécule dans l'air est de l'ordre de $500ms^{-1}$, son énergie cinétique est donc, en joules,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vmoy = 500\n",
    "E0 = 0.5 * m_azote * vmoy ** 2\n",
    "print(E0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prenons cette valeur pour $E_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def temp(q):\n",
    "    kb = 1.3807e-23\n",
    "    return E0 / kb / math.log(1 + 1 / q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tentons avec $q=1$. N'attendons aucun miracle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp(1) - 273.15"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Houla, 334 degrés Celsius, c'est chaud  :-) ! Mais ce n'est pas $10^{-50}$ ou $10^{312}$, on trouve un nombre qui ressemble vraiment à une température \"normale\". \n",
    "\n",
    "J'aurais adoré trouver 25 degrés Celsius, la température de l'air par un beau jour du joli mois de mai, mais avec notre modèle ultra-simpliste c'était peine perdue ... Si vous voulez vraiment épater le monde scientifique, changez les limitations de vitesse pour les molécules et prenez une vitesse de $350ms^{-1}$ :-)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vmoy = 350\n",
    "E0 = 0.5 * m_azote * vmoy ** 2\n",
    "temp(1) - 273.15"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__NB__ : Eh oui, $25^oC$. Mais c'est pas beau de tricher :-)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Le moment de vérité"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous y voilà. Le moment fatidique c'est celui où l'on compare théorie et pratique. Affichons sur un même graphique l'histogramme de la liste $S$ et la courbe de $P(\\epsilon)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stats3(S, q, N):\n",
    "    h = histogramme(S)\n",
    "    plt.plot(h, color='k')\n",
    "    xs = [proba(k, q) for k in range(len(h) + 1)]\n",
    "    plt.plot(xs, 'r')\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Allons-y.\n",
    "- En noir, l'histogramme de notre jeu.\n",
    "- En rouge, les probabilités prévues par la théorie."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 3\n",
    "N = 10000\n",
    "S = iterer(N, q, 500000)\n",
    "stats3(S, q, N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous sommes satisfaits au delà de tout ce que nous avions pu espérer :-). Et encore, j'ai pris un nombre d'itérations pas trop grand. Parce que si vous augmentez le nombre d'itérations, vous ne verrez plus du tout la courbe noire (qui a été tracée en premier par `matplotlib`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Bibliographie__ : Pour ce projet je me suis fortement inspiré du chaptire 4 du livre :\n",
    "\n",
    "Concepts in Thermal Physics, S. Blundell et K. Blundell, _Oxford University Press (2006)_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
