Preuve assitée par ordinateur

Au mois de Mai dernier H. A. Helfgott  a publié une preuve de la version faible de la conjecture de Goldbach qui a fait grand bruit. La preuve repose sur un résultat antérieur (le Théorème de Vinogradov) affirmant que  la conjecture est vraie pour tout $n\geq n_0$  avec $n_0=e^{3100}\approx 2 \times 10^{1346}$ mais l'auteur améliore la borne $n_0$ jusqu'à l'amener à $10^{30}$ rendant ainsi possible une vérification exhaustive des cas manquant ... en d'autre termes cette preuve repose  en partie sur l'utilisation d'un ordinateur !

La question de l'utilisation d'un outil "mécanique" pour effectuer une preuve (en entier ou en partie) est un sujet assez polémique dans le monde des mathématiques, comme ce fut le cas pour  la preuve du théorème des quatre couleurs proposée  en 1976 par K. Appel et W. Haken . Pour vous inviter à réfléchir sur le sujet je vous propose d'étudier la  preuve  assistée par ordinateur d'un résultat étonnant sur la série $\sum_{k\in{\mathbb Z}}e^{-k^2/2}$ :

$$\sum_{k=-10}^{10}e^{-k^2/2}\approx 2.5066283...\approx \sqrt{2\pi}~~~ \textrm{à}~~ 10^{-8}~~\textrm{près}$$
pourtant $\sum_{k\in{\mathbb Z}}e^{-k^2/2}\neq \sqrt{2\pi}$  plus précisément on peut  montrer que :
$$10^{-8}<\sum_{k=-\infty}^{+\infty}e^{-k^2/2}-\sqrt{2\pi}<2\times10^{-8}$$


La formule d'Euler-MacLaurin (que j'ai donné dans un billet précédent) permet de montrer facilement que la série ne vaut pas $\sqrt{2\pi}$. En appliquant cette formule pour $a=-\infty$ et $b=\infty$ on obtient :
$$ \sum_{k=-\infty}^{\infty}e^{-k^2/2}=\int_{-\infty}^\infty e^{-t^2/2}dt -\int_{-\infty}^\infty \{t\}te^{-t^2/2}dt>\sqrt{2\pi}$$
en remarquant que  :

• la première intégrale est l'intégrale de Gauss qui vaut  $\sqrt{2\pi}$
• la deuxième intégrale, bien que très petite en valeur absolue, est non nulle

Pour montrer le deuxième point l'idée est d'évaluer l'intégrale par une méthode numérique  tout en contrôlant l'erreur de telle sorte qu'on puisse affirmer que le terme de reste  est bien de l'ordre de $10^{-8}$.

Les méthodes d'intégration numérique les plus simples consistent à découper l'intervalle d'intégration en m sous-intervalles sur lesquels on approche la fonction par un polynôme de degré 0 (méthode des rectangles) 1 (méthode des trapèzes) ou 2 ( méthode de Simpson ) . Cette dernière méthode est particulièrement efficace puisque l'erreur commise par la  valeur approchée I de l'intégrale est majorée par la formule:

$$\left| I - \int_{a}^bf(t)dt \right| \leq \frac{M\times (b-a)^5}{2880m^4}$$
où $M=\sup_{[a,b]}\vert f''''(t)\vert $ le maximum de la valeur absolue de la dérivée d'ordre 4 de la fonction intégrée sur $[a,b]$. Pour l'appliquer ici il faut donc

• choisir un intervalle $[a,b]$ assez large  pour que le reste de l'intégrale soit négligeable
• choisir le niveau de discrétisation de $[a,b]$ pour que l'erreur de calcul de l'intégrale sur cet intervalle soit faible

Le choix de l'intervalle est assez simple en utilisant que pour $n\geq 1$   :
$$ \int_{n}^\infty e^{-t^2/2}dt\leq \int_{n}^\infty te^{-t^2/2}dt=\left[-e^{-t^2/2}\right]_n^\infty=e^{-n^2/2}$$
 On peut donc se ramener à $[a,b]=[-7,7]$ en sachant que le reste est  négligeable
$$ \left\vert \int_{[-7,7]^c}t\{t\} e^{-t^2/2}dt\right\vert\leq 2\int_{7}^\infty te^{-t^2/2}dt=2 e^{-7^2/2}<2\times 10^{-10}$$
Il faut ensuite calculer les dérivées de $f(t)= t\{t\}e^ {- {{t^2}\over{2}} }$ :

• $f'(t)=e^ {- {{t^2}\over{2}} }\,\left((1-t^2 )\{t\}+t\right)$
• $f''(t)=e^ {- {{t^2}\over{2}} }\,\left((t^3-3\,t)\{t\}-2\,t^2+2\right)$
• $f'''(t)=-e^ {- {{t^2}\over{2}} }\,\left((t^4-6\,t^2\,+3\,)\{t\}-3\,t^3+9\,t\right)$
• $f''''(t)=e^ {- {{t^2}\over{2}} }\,\left((t^5\,-10\,t^3\,+15\,t\,)\{t\}-4\,t^4+24\,t^2-12\right)$

Pour faire ces obtenir ces différentes formules  j'ai utilisé le logiciel de calcul formel WxMaxima , il faut être astucieux pour dériver les expressions contenant $u(t)=\{t\}$, dont la dérivée (en dehors des points entiers) vaut 1 :

f(t):=t*u(t)*exp(-t^2/2);
define(f1(t),subst(1,diff(u(t),t),diff(f(t),t)));
define(f2(t),subst(1,diff(u(t),t),diff(f1(t),t)));
define(f3(t),subst(1,diff(u(t),t),diff(f2(t),t)));
define(f4(t),subst(1,diff(u(t),t),diff(f3(t),t)));
factor(f1(t));factor(f2(t));factor(f3(t));factor(f4(t));

si maintenant on trace le graphe de $f''''$  (par exemple avec Scilab!)  on se rend compte qu'on a   $M\leq 14$   :

donc pour avoir une erreur inférieur à $10^{-9}$ il suffit de prendre   :
$$m\geq (14\times 14^5\times10^9/2880)^{1/4}\approx  1271.5807  ....$$
On peut alors calculer la valeur approchée de l'intégrale avec un logiciel comme Scilab  et le script suivant :

function y=f(x) 
     t=pmodulo(x,1)
     y=t*x*exp(-x^2/2) 
endfunction 
x=[-7:0.01:7];
y=feval(x,f);  
I=intsplin(x,y) 

Qui nous donne la valeur $I=-1.339\times 10^{-8}$, en tenant compte de l'erreur  commise  en dehors de l'intervalle $[-7,7]$  et de l'erreur de la méthode de Simpson on peut donc affirmer que

$$\int_{-\infty}^\infty \{t\}te^{-t^2/2}dt\approx 1.34\times 10^{-8}~~\textrm{à}~~10^{-9} ~~\textrm{près}$$
ce qui suffit à justifier le résultat.

La morale de cette histoire

Cette  démonstration repose donc sur un calcul numérique que je n'aurai pas pu faire à la main,  cependant l'emploi d'un programme pour effectuer cette vérification nécessite aussi l'utilisation de mathématiques pour en justifier la validité ....  pourquoi s'en passer à notre époque où le calcul numérique sur ordinateur est à la porté  de tous? D'un autre coté  l'utilisation systématique d'ordinateurs pour les calculs  peut s'avérer dangereuse  : ici  déduire la limite de la série à partir de la valeur approchée obtenue d'une somme partielle conduit à une mauvaise conjecture ! Pour autant c'est un peu de cette manière qu' Euler a pu conjecturer la solution du problème de Bâle $\sum_{k=1}^{+\infty}{1\over k^2}={\pi^2\over 6}$ ... mais il lui fallu encore près de 10 ans pour en donner une preuve rigoureuse!