Une curiosité numérique
Complément numérique

Cet article complète l’article paru dans le n°548 d’Au fil des maths .

François Boucher

© APMEP Juin 2023
⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅♦⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅

\[\sum_{k=1}^{+\infty}\frac{(-1)^{k-1}}{2k-1}=1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+ \cdots=\frac{\pi}{4}\cdotp\]

On pose \(\displaystyle s_n=4\sum_{k=1}^n\frac{(-1)^{k-1}}{2k-1}\) et on a l’encadrement (si \(n\) est pair) : \(s_{n}< \pi < s_{n}+\dfrac{4}{2n+1}\cdotp\)

Un programme de calcul simple — dans lequel on accumule la somme en commençant par les termes les plus petits — fournit toutefois les valeurs approchées suivantes :

def Leibniz(lim) :

  s=0.
  n=lim-1
  sgn=-4
  while n>0:
    s=s+sgn/n
    n=n-2
    sgn=-sgn

  return(s)

On obtient le tableau suivant où les valeurs approchées de \(\pi\) sont a priori connues.

On note bien que, si \(n=10^p\), la \(p\)-ième décimale n’est pas exacte à une unité par défaut près et que la différence \(\overline{\pi} – \overline{s_n}\) vérifie \(0<\overline{\pi} – \overline{s_n} \approx 10^{-p}<2\times10^{-p}\). On a indiqué dans la colonne des valeurs approchées de \(s_n\) la (\(3p+1\))-ième décimale calculée.

Démontrons, ainsi que permet de le conjecturer la petite expérience numérique précédente que \(\pi-s_n \sim \dfrac{1}{n}\) (si \(n\) est pair). On remplace \(n\) par \(2n\) pour traduire l’hypothèse \(n\) pair.

Cette dernière inégalité montre que, dans le cas \(2n=10^p\), l’erreur sur la \(p\)-ième décimale est bien de \(1\) par défaut, et que les \(2p\) décimales suivantes sont exactes. Ainsi, avec \(p=5\), on peut calculer une valeur approchée de \(\pi\) avec quinze décimales exactes à l’aide d’un système de calcul utilisant des flottants IEEE754 doubles.

Reprenons le calcul précédent avec \(n=50\,000\) en calculant en multiple précision (par exemple, avec le module Decimal de Python ou le logiciel Maple). On travaille avec \(\dfrac{\pi}{2}\) et on choisit \(n\) de sorte que \(2n\) soit une puissance de dix.

from decimal import getcontext, Decimal

getcontext().prec = 100 # nombre de décimales du significande

def conv(z):
""" convertit un type numeric en type Decimal"""
  return(Decimal(str(z)))

  def Leibniz(lim):
  """ calcule (*$\sum_1^{k=2lim} 2\frac{(-1)^{k-1}}{2*k-1}$*)"""
  
  # initialisations
  n=2*lim-1
  s=conv(0.)
  sgn=conv(-2)
  
  # boucle de sommation inversée
  while(n>0):
    s+=sgn/conv(n)
    n-=2
    sgn=-sgn
  return(s)

On obtient alors : \(\overline{s}=1,570\,7{\color{red}8}6\,326\,794\,89{\color{red}7}\,619\,231\,321\,{\color{red}1}91\,639\,75{\color{red}2}\,{\color{red}{05}}2\,098\) (en rouge, les décimales erronées).

Avec \(n=5\,000\,000\), la durée de calcul reste acceptable.

\[\begin{aligned}
\overline{s}=1,&570\,796\,{\color{red}2}26\,794\,896\,619\,23{\color{red}2}\,321\,691\,639\,75{\color{red}1}\,{\color{red}{39}}2\,098\,584\,699\,6{\color{red}{93}}\,\\
&{\color{red}6}52\,910\,487\,47{\color{red}0}\,{\color{red}{911}}\,153\,908\,203\,{\color{red}{648}}\,{\color{red}{31}}4\,499\,31{\color{red}3}\,{\color{red}{747}}\,{\color{red}{136}}\,{\color{red}1}71\,058
\end{aligned}\]

Un développement asymptotique

Démontrons le développement annoncé dans le n°548  (\(n\) pair) \[\frac{\pi}{2}-2\sum_{k=0}^{n-1}\dfrac{(-1)^{k}}{2k+1}= \sum_{k=0}^{m-1}\dfrac{\mathrm{E}_{2k}}{(2n)^{2k+1}}+\mathop{\mathrm{O}}_{n \infty}\left(\frac{1}{n^{2m+1}}\right)\] valide pour tout \(m \geqslant 1\). Les nombres d’Euler \(E_{2n}\) ont été définis dans la version papier  paru dans Au fil des maths.

On a besoin de quelques matériaux qui s’obtiennent plus aisément en faisant appel aux séries entières, bien plus facilement avec les fonctions analytiques (variable complexe) mais c’est plus cher, et, plus laborieusement mais élémentairement, avec des développements limités. La théorie fournit existence et unicité, et les outils de calcul, de quoi faire un bon sujet de concours (par exemple ENSI M 1983 épreuve 1). Le raisonnement par récurrence est d’usage constant.

1. les polynômes d’Euler : pour \(x \in \mathbb{R}\), l’application \(t \longmapsto \dfrac{2\mathrm{e}^{xt}}{\mathrm{e}^t+1}\) possède un développement en série entière (dit exponentiel) de la forme \(\displaystyle \sum_{n=0}^{+\infty}\dfrac{\mathrm{E}_n(x)}{n{!}}t^n\) avec un rayon de convergence \(R \geqslant 1\) (on peut montrer que \(R=\pi\)), et dans lequel \(\mathrm{E}_n(x)\) est un polynôme unitaire en \(x\) de degré \(n\).

   \[\mathrm{E}_0(x)=1,\ \mathrm{E}_1(x)=x-\frac{1}{2},\ \mathrm{E}_2(x)=x^2-x,\ \mathrm{E}_3(x)=x^3-\frac{3}{2}x^2+\frac{1}{4},\ \mathrm{E}_4(x)=x^4-2x^3+x\]

2. Les \(\mathrm{E}_n(x)\) vérifient de nombreuses propriétés, dont les plus utiles pour la suite sont (en omettant les quantifications) :
   • \(\displaystyle \mathrm{E}_n(x)=x^n-\dfrac{1}{2}\sum_{k=0}^{n-1}\binom{n}{k}\mathrm{E}_k(x)\) ;
   • \(\mathrm{E}_n(0)+\mathrm{E}_n(1)=0\) pour \(n \geqslant 1\), \(\mathrm{E}_0(x)=1\), \(\mathrm{E}’_{n+1}(x)=(n+1)\mathrm{E}_n(x)\) qui caractérisent la famille \((\mathrm{E}_n(x))_{n \in \mathbb{N}}\) ;
   • lien avec les nombres d’Euler : \(E_n=2^n\mathrm{E}_n\left(\dfrac{1}{2}\right)\) ;
   • \(\mathrm{E}_n(x)+\mathrm{E}_n(x+1) = 2x^n\) elle-aussi caractéristique ;
   • \(\mathrm{E}_n(1-x)=(-1)^n\mathrm{E}_n(x)\) qui donne \(\mathrm{E}_{2n+1}\left(\dfrac{1}{2}\right)=0\) et \(\mathrm{E}_{2n}(0)=\mathrm{E}_{2n}(1)=0\) ;
   • En utilisant ce qui précède, on démontre (par récurrence !) qu’il y a quatre types de représentations graphiques entre \(0\) et \(1\) pour les polynômes \(\mathrm{E}_n\) selon \(n\) mod 4 :

     

     En particulier, on a donc \(\displaystyle \max_{[0;1]}\left|\mathrm{E}_{2n}\right|=\left|\mathrm{E}_{2n}\left(\dfrac{1}{2}\right)\right|=\dfrac{\mathrm{E}_{2n}}{2^{2n}}\cdotp\)

3. On définit une fonction périodique \(\widetilde{\mathrm{E}}_n\), de période 2 en posant
   \(\widetilde{\mathrm{E}}_n(x)=\mathrm{E}_n(x)\) pour \(0 \leqslant x < 1\) et \(\widetilde{\mathrm{E}}_n(1+x)=-\widetilde{\mathrm{E}}_n(x)\) pour tout \(x\).
   On démontre que, pour \(n\geqslant 1\), \(\widetilde{\mathrm{E}}_n\) est de classe \(\mathcal{C}^{n-1}\) sur \(\mathbb{R}\) (les raccords se font bien en \(0\) et \(1\)).

Démontrons ensuite la version de la formule de Boole retenue : soit \(x \in \mathbb{R}\) et \(m\) entier non nul ; si \(f\) est de classe \(\mathcal{C}^{\infty}\) sur \([x;x+1]\) alors pour \(h \in [0;1]\) et tout entier \(m\geqslant 1\), on a \[f(x+h)=\sum_{k=0}^{m-1} \dfrac{1}{k{!}}\mathrm{E}_k(h)\frac{1}{2}\left\{f^{(k)}(x+1)+f^{(k)}(x)\right\} +\frac{1}{2}\int_0^1\dfrac{\widetilde{\mathrm{E}}_{m-1}(h-t)}{(m-1){!}}f^{(m)}(x+t)\,\mathrm{d}t.\]

On en déduit une formule sommatoire : on prend \(h=\dfrac{1}{2}\) et \(x=p \in \mathbb{N}^{\ast}\) et \(m=2M+1\), il vient \[\begin{aligned}
2(-1)^pf\left(p+\frac{1}{2}\right) &= \sum_{k=0}^{2M}\frac{1}{k{!}}\mathrm{E}_k\left(\frac{1}{2}\right)\left\{(-1)^pf^{(k)}(p)-(-1)^{p+1}f^{(k)}(p+1)\right\} +(-1)^p\int_0^1\frac{\widetilde{\mathrm{E}}_{2M}\left(\dfrac{1}{2}-t\right)}{(2M){!}}f^{(2M+1)}(p+t)\,\mathrm{d}t \\
&= \sum_{h=0}^{M}\frac{1}{2^{2h}(2h){!}}\mathrm{E}_{2h}\left\{(-1)^pf^{(2h)}(p)-(-1)^{p+1}f^{(2h)}(p+1)\right\} +\int_p^{p+1}\frac{\widetilde{\mathrm{E}}_{2M}\left(\dfrac{1}{2}-u\right)}{(2M){!}}f^{(2M+1)}(u)\,\mathrm{d}u
\end{aligned}\] puis, en supposant valides toutes les convergences utiles \[\begin{aligned}
2\sum_{p=n}^{+\infty}(-1)^pf\left(p+\frac{1}{2}\right) = \sum_{h=0}^{M}\frac{1}{2^{2h}(2h){!}}\mathrm{E}_{2h} & \left\{(-1)^nf^{(2h)}(n)-\lim_{p \infty}(-1)^{p+1}f^{(2h)}(p+1)\right\} \\
& +\int_n^{+\infty}\frac{\widetilde{\mathrm{E}}_{2M}\left(\dfrac{1}{2}-u\right)}{(2M){!}}f^{(2M+1)}(u)\,\mathrm{d}u
\end{aligned}\] en prenant \(f(x)=\dfrac{1}{x}\), on valide les convergences et on obtient \[2\sum_{p=n}^{+\infty}(-1)^p\frac{1}{2p+1}=(-1)^n\sum_{h=0}^{M}\frac{\mathrm{E}_{2h}}{(2n)^{2h+1}}- \frac{1}{2}\int_n^{+\infty}\widetilde{\mathrm{E}}_{2M}\left(\frac{1}{2}-u\right)\frac{2M+1}{u^{2M+2}}\,\mathrm{d}u\cdotp\]

Nous avons vu que \(|\widetilde{\mathrm{E}_{2M}}|\) est majorée par \(\dfrac{|\mathrm{E}_{2M}|}{2^{2M}}\) ; on en déduit que \[\left|\frac{1}{2}\int_n^{+\infty}\widetilde{\mathrm{E}}_{2M}\left(\frac{1}{2}- u\right)\frac{2M+1}{u^{2M+2}}\,\mathrm{d}u\right| \leqslant \frac{|\mathrm{E}_{2M}|}{2^{2M+1}}\int_n^{+\infty}\frac{2M+1}{u^{2M+2}}\,\mathrm{d}u =\frac{|\mathrm{E}_{2M}|}{(2n)^{2M+1}}\]

En supposant \(n\) pair, on a donc \[2\sum_{p=n}^{+\infty}(-1)^p\frac{1}{2p+1}=\sum_{h=0}^{M}\frac{\mathrm{E}_{2h}}{(2n)^{2h+1}}+r_{2M}\] avec \(0 \leqslant r_{2M} \leqslant \dfrac{|\mathrm{E}_{2M}|}{(2n)^{2M+1}}\), ce qui, à \(M\) fixé, valide le développement asymptotique.

Calcul de 1 000 décimales de \(\pi\) avec la formule de Leibniz

Le problème est donc de trouver le \(n\) et le \(M\) convenables ; la condition suffisante est \[\dfrac{2|\mathrm{E}_{2M}|}{(2n)^{2M+1}} < \frac{1}{2}10^{-1\,000}.\]

On cherche \(n\) et \(M\) avec \(n \gg M\) car le calcul des \(\mathrm{E}_{2k}\), bien que non coûteux en temps, l’est en mémoire ; on peut démontrer que \(|\mathrm{E}_{2k}| \sim \dfrac{4^{k+1}(2k){!}}{\pi^{2k+1}}\) qu’on va utiliser pour estimer l’ordre de grandeur de \(|\mathrm{E}_{2k}|\). En utilisant des logarithmes décimaux, on obtient le tableau

Une exploration numérique donne comme valeurs possibles \(M=118\) et \(2n=10^6\).

On propose un programme en Python avec le module Decimal. Le programme est naïf au sens où on ne cherche pas à optimiser les calculs ; ainsi on pourrait ne pas calculer la liste de tous les \(\mathrm{E}_{2h}\) avant de calculer le terme correctif. L’intention est de privilégier la lisibilité globale.

def pi(n,M,d) :
  """ calcule (*$\pi$*) avec d décimales """
  """ n et M satisfaisant la CS """

  def Euler(M):
    """ construit la liste (*$[\mathrm{E}_{2h}, 0\leqslant h \leqslant M]$*)"""
    N=2*M+1
    L=[]

    for h in range(0,N):
      L.append([0 for k in range(0,h+1)])
      L[0][0]=1
      for h in range(1,N):
      for k in range(1,h+1):
        L[h][k]=L[h][k-1]+L[h-1][h-k]
    \mathrm{E}=[(1-2*(h%2))*L[2*h][2*h] for h in range(M+1)]

    return(\mathrm{E})

  def conv(z):
    """ convertit un type numeric en type Decimal"""
    return(Decimal(str(z)))

  def Leibniz(lim):
    """ calcule (*$\sum^1_{k=2lim} 4\frac{(-1)^{k-1}}{(2k-1)}$*)"""
    s=conv(0.)
    n=2*lim-1
    sgn=conv(-4)

    while(n>0):
      s+=sgn/conv(n)
      n-=2
      sgn=-sgn
    return(s)

  def cor(n,M):
    """ calcule (*$\sum_{h=0}^{M}2\frac{\mathrm{E}_{2h}}{(2n)^{2h+1}}$*)"""
    T=Euler(M)
    h=conv(2*n)
    C=conv(T[0])/h
    H=h*h

    for k in range(1,M+1):
      h=h*H
      C+=conv(T[k])/h
    C=2*C

    return(C)

  getcontext().prec = d # nombre de décimales des significandes
  return(Leibniz(n)+cor(n,M))

pi(500000,118,1010) donne

Decimal(’ 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067
  982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381
  964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141
  273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511
  609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011
  949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669
  405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689
  258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816
  096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177
  669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909
  2164201989380952572’)

décimales toutes exactes (a posteriori).

D’autre part, l’équivalent donné de \(\mathrm{E}_{2M}\) prouve que, à \(n\) fixé, la série \(\displaystyle \sum_M \frac{\mathrm{E}_{2M}}{(2n)^{2M+1}}\) est très grossièrement divergente ; ce qui n’empêche pas ses sommes partielles de fournir des termes correctifs efficaces.

⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅♦⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅

François Boucher, à la retraite depuis quelques années, continue de s’intéresser aux mathématiques et à leur enseignement. Il fait également partie de l’équipe d’Au fil des maths.