La formule de Stirling

Quand on me demande ce qui m'a donné envie de faire des mathématiques il me vient toujours  deux formules à l'esprit : la série $\sum_{k\geq 1}{1\over k^2}={\pi^2\over 6}$ (question aussi appelée "problème de Bâle") et la formule de Stirling:

$$n!=1\times 2\times \dots\times n\sim_\infty \left({n\over e}\right)^n\sqrt{2\pi n}
\Leftrightarrow
\lim_{n\to\infty} {n!\over \left({n\over e}\right)^n\sqrt{2\pi n}}=1$$

le fait qu'un objet  défini uniquement à partir de nombres entiers puisse conduire à une formule contenant $\pi$ (et même combinant $e,\pi$ et $\sqrt{~}$ dans le cas de $n!$)  m'a tout de suite fasciné, et en rentrant à l'université  je voulais absolument comprendre "la" démonstration  de telles curiosités mathématiques. Vingt-cinq ans après je suis toujours impressionné de la richesse de ces deux problèmes, l'étude de leurs démonstrations conduit toujours à l'introduction de méthodes  mathématiques fondamentales que je décris souvent sur ce blog. Laissez moi vous en donner un petit aperçu dans le cas de formule de Stirling.

Comparaisons séries intégrales

la première méthode pour analyser le comportement asymptotique de $n!$  consiste à se ramener à l'étude d'une série en passant au logarithme :

$$\ln(n!)=\underbrace{\ln(1)}_{=0}+\ln(2)+\dots+\ln(n)=\sum_{k=2}^n\ln(k)$$

en regardant cette série comme une somme d'aires de rectangles de base $[k,k+1]$ et de hauteur $\ln(k)$ (ou $\ln(k+1)$)  on obtient immédiatement que :
  $$ \ln(n!)= \sum_{k=2}^n\ln(k)\geq\int_1^n\ln(t) dt= n(\ln(n)-1)+1=\ln\left(\left({n\over e}\right)^ne\right)\geq \sum_{k=1}^{n-1}\ln(k)=\ln(n!)-\ln(n)$$
on en déduit en repassant aux exponentielles :

$$ \left({n\over e}\right)^ne\leq  n!\leq \left({n\over e}\right)^nen$$

Cet encadrement donne correctement  le terme $(n/e)^n$  mais on peut aussi obtenir le facteur $\sqrt{n}$ en remplaçant les rectangles par des trapèzes judicieusement choisis (voir pdf) :

• en considérant les trapèzes de largeur  $[k,k+1]$ et de petite (resp grande) base  $\ln(k)$ (resp. $\ln(k+1)$) on obtient que $n!\leq \left({n\over e}\right)^ne\sqrt{n}$ via :
  $$\int_1^n \ln(t) dt=n\ln(n)-n+1\geq \sum_{k=1}^{n-1} {\ln(k)+\ln(k+1)\over 2}= \ln(n!)-{\ln(n)\over 2}$$
• en considérant les trapèzes ayant pour cotés   $[k,k+1]$ et son image par la tangente au graphe de $x\mapsto \ln(x)$ en $x=k$ (équation $T(x)=\ln(k)+{1\over k}(x-k)$) dont l'aire est donnée par :
  $$ \int_k ^{k+1}T(x) dx=\ln(k)+{1\over 2k}$$
  on obtient alors que $n!\geq \left({n\over e}\right)^n\sqrt{2en}$ via l'inégalité :
  $$\int_1^n \ln(t) dt=n\ln(n)-n+1\leq \sum_{k=1}^{n-1} \ln(k)+{1\over 2k}= \ln(n!)-\ln(n)+{\ln(n)+1-\ln(2)\over 2}$$

Ces encadrements  donnent déjà une approximation correcte du facteur $C=\sqrt{2\pi}\approx  2.5066283  $
$$2.331644 \approx\sqrt{2 e} \leq C\leq e\approx 2.7182818 $$
pourtant ils ne nous assurent pas  que $n!\sim_\infty C(n/e)^n \sqrt{n}$ (en effet $n!$ pourrait osciller entre les deux bornes), Pour cela il faut une méthode de comparaison série/intégrale  plus efficace. Si on cherche un peu on finit forcément par aboutir à la formule d'Euler-MacLaurin d'ordre 2 qui donne directement le résultat :

$$\begin{align*}
\sum_2^n \ln(k)&=  \int_{1}^n \ln(t) dt+ {\ln(n)-\ln(1)\over 2}-\int_{1}^n -{\beta(t)\over t^2}  dt\\
&= \underbrace{n(\ln(n)-1)+ {\ln(n)\over 2}}_{=\ln\left(\left({n\over e}\right)^n\sqrt{n}\right)}
+\underbrace{1+\int_{1}^\infty {\beta(t)\over t^2}  dt}_{=\ln(C)}
-\underbrace{\int_{n}^\infty {\beta(t)\over t^2}  dt}_{=O\left({1\over n}\right)}\\
\Rightarrow n!&=\left({n\over e}\right)^n\sqrt{n}\times C
\times e^{O\left({1\over n}\right)}
\end{align*}
$$
pour les justifications : la suite d'intégrales  $\int_{1}^n {\beta(t)\over t^2}  dt$ converge en $n\to\infty$ et le "reste"  est bien en $O\left({1\over n}\right)$ car $\beta(t)= {\{t\}^2-\{t\}\over 2}$  est bornée  $(\{t\}$  est la partie décimale de $t$). La formule de Stirling  est alors démontrée puisque $e^{O\left({1\over n}\right)}$ tend vers 1, mis à part qu'il reste à voir que  la valeur de  la constante  est :
$$C=1+\int_{1}^\infty {\beta(t)\over t^2}  dt=\sqrt{2\pi}$$
cela se fait le plus souvent à l'aide de calculs d'intégrales permettant de relier les équivalents de $(n!)^2$  et $(2n)!$ car
$$ {(2n)!\over (n!)^2}\sim_\infty \left({2n\over e}\right)^{2n}\sqrt{2n} C\times
\left({e\over n}\right)^{2n}{1\over nC^2}={4^n\over C}\sqrt{2\over n}$$
Le moyen le plus célèbre d'y arriver est le calcul des intégrales de Wallis . Une série de démonstration par récurrence permet de trouver :

• d'un coté $W_n=\int_0^{\pi/2}\sin(t)^n dt\sim_\infty \sqrt{\pi\over 2n}$
• de l'autre $W_{2n}={(2n)!\over 4^n(n!)^2}{\pi\over 2}\sim_\infty \sqrt{\pi^2\over 2nC^2}$
• au final on a $\sqrt{\pi\over 4n}=\sqrt{\pi^2\over 2nC^2}\Rightarrow C=\sqrt{2\pi}$

Fonction Gamma et méthode de Laplace

Comme on l'a vu ci-dessus écrire une  formule "discrète" sous forme d'intégrale permet  d'accéder à des méthodes d'analyse  très efficaces pour l'étude asymptotique. C'est pour ainsi dire la base de toute la théorie des nombres moderne. Pour $n!$  la "bonne" formule intégrale est la suivante :

$$n!=\int_0^\infty t^n e^{-t} dt=\Gamma(n+1)$$

Cette fonction $\Gamma$ est fondamentale dans de nombreuses parties des mathématiques. La formule se démontre facilement en répétant des intégrations par parties pour diminuer l'exposant de $t^n$ jusqu'à arriver à $\int_0^\infty e^{-t}dt=1$ :

$$\begin{align*}
I_n&=\int_0^\infty t^n e^{-t} dt\\
&=\underbrace{\left[-t^ne^{-t}\right]_0^\infty}_{=0}
-\int_0^\infty -nt^{n-1}e^{-t} dt\\
&=nI_{n-1}=\dots=n\times(n-1)\times\dots\times 1\times I_0=n!
\end{align*}$$
La recherche  d'un équivalent pour $n!$  amène donc à se poser des questions sur l'étude du comportement d'une intégrale  par rapport à un paramètre.
$$n!=\int_0^\infty t^n e^{-t} dt=\int_0^\infty e^{n\ln(t)-t} dt=\int_0^\infty e^{n\phi(t)} dt$$
si on trace le graphe de $e^{n\ln(t)-t}$  on obtient une courbe en cloche  qui tend vers 0  aux deux extrémités de l'intervalle d'intégration,  avec un maximum en $t=n$ autour duquel on doit avoir que
$$n\ln(t)-t=n\left(\ln(t) -{t\over n}\right)=n\phi(t)=n\phi(n)+0+b(t-n)^2+...$$
 donc quitte à faire un  changement de variable $t\to u$ tel que $n\ln(t)-t=n\phi(n)-n u^2$ on peut se ramener à une intégrale Gaussienne :
$$n!=\int_{\mathbb R} e^{n\phi(n)} e^{-n u^2}g(u) du$$
le terme $ e^{n\phi(n)}$ peut être sorti de l'intégrale, car constant en $u$, il reste ensuite à faire un nouveau changement de variable  $s=\sqrt{n}u$  pour obtenir que
$$n!=e^{n\phi(n)} \int_{\mathbb R}  e^{-s^2}g\left({s\over \sqrt{n}}\right){ds\over \sqrt{n}}$$
maintenant si la fonction $g$  est bornée  sur ${\mathbb R}$  on obtiendra par convergence dominée que
$$ n!\sim_\infty e^{n\phi(n)} \int_{\mathbb R}  e^{-s^2}g(0){ds\over \sqrt{n}}= {e^{n\phi(n)}\over \sqrt{n}}\sqrt{\pi}g(0)$$
Cette stratégie  est très générale  et porte le nom de méthode de Laplace   :

Théorème :  soit $\phi$ une fonction $C^2$ ayant un unique maximum sur $[a,b]$  en $t_0$ alors
$$\int_a^b e^{n\phi(t)}\, dt\sim_\infty \sqrt{\frac{2\pi}{n|\phi''(t_0)|}}e^{n\phi(t_0)} $$

En prenant $\phi(t)=\ln(t) -{t\over n}$ on a
$$ \phi'(t)={1\over t} -{1\over n}=0\Rightarrow t_0=n\Rightarrow \phi''(t_0)=-{1\over n^2}$$
donc
$$n!\sim_\infty e^{n(\ln(n)-1)}\sqrt{2\pi\over {1\over n^2}n}= \left({n\over e}\right)^n\sqrt{2\pi n} $$

sauf que dans cet exemple précis $\phi(t)$ et $t_0$ dépendent de $n$  ce qui n’apparaît pas dans les hypothèses de la méthode de Laplace! Pourtant ce n'est pas grave, la preuve du théorème peut prendre en compte le cas où $\phi$ dépend de $n$ , mais on peut aussi se passer de la méthode de Laplace et construire à la main les différents changements de variables :

• on commence par poser $t=ne^z$  (donc $dt=ne^zdz$)
  $$n!=\int_0^\infty e^{n\ln(t)-t} dt=\int_0^\infty e^{n\ln(n)+n(z-e^z)} ne^z dz={n^{n+1}\over e^n}\int_0^\infty e^{-n(e^z-z-1)} e^z dz$$
• la fonction $\Phi(z)={\rm sign}(z)\sqrt{e^z-z-1}$ est bien $C^1$ puisque en $z=0$ les fonctions $\Phi$ et $\Phi'$ se raccordent continûment comme le montre le DL:
  $$\Phi(z)={\rm sign}(z)\sqrt{z^2/2+O(z^3)}= {z\over \sqrt{2}}+O(z^2)\Rightarrow \Phi'(0)={1\over \sqrt{2}}$$
  on peut donc faire le changement de variable $u=\Phi(z)$ (donc $dz={du \over (\Phi^{-1})'(u)}$) :
  $$n!={n^{n+1}\over e^n} \int_{-\infty}^\infty  e^{-nu^2} {e^{\Phi^{-1}(u)}\over (\Phi^{-1})'(u)} du$$
• enfin le changement de variable $s=\sqrt{n}u$ donne
  $$n!={n^{n+1}\over e^n\sqrt{n}} \int_{-\infty}^\infty  e^{-s^2} {e^{\Phi^{-1}(u/\sqrt{n})}\over (\Phi^{-1})'(u/\sqrt{n})} ds
  \sim_\infty {n^{n}\over e^n} \sqrt{\pi n} {e^{\Phi^{-1}(0)}\over (\Phi^{-1})'(0)} $$
• comme $\Phi(0)=0$  et $(\Phi^{-1})'(\Phi(z))={1\over \Phi'(z)}$ on en déduit que :
  $$ {e^{\Phi^{-1}(0)}\over (\Phi^{-1})'(0)} = {e^0\over {1\over \sqrt{2}}}=\sqrt{2}$$