Les scripts du chapitre

Le laboratoire du cours donne deux programmes courts ; ils sont ici complétés, exécutés et commentés. Le devoir en propose un troisième, à trous ; il est ici livré entier. Chaque bloc se copie tel quel dans votre éditeur, et chaque sortie affichée est une sortie réelle, recopiée telle quelle.

1. Archimède, la version du cours

Le principe est celui de la section 6.1 : le demi-périmètre du polygone régulier à $n$ côtés inscrit dans le cercle de rayon $1$ vaut $n\sin\frac{\pi}{n}$, et ce nombre tend vers $\pi$ parce que $\frac{\sin h}{h} \longrightarrow 1$. On part de l’hexagone, où $\sin\frac{\pi}{6} = \frac12$ se connaît exactement, et l’on double le nombre de côtés grâce à la formule de l’angle moitié

$$\forall \theta \in \left[0\,;\frac{\pi}{2}\right], \quad \sin\frac{\theta}{2} = \sqrt{\frac{1 - \sqrt{1 - \sin^2\theta}}{2}},$$

valable sur cet intervalle parce que $\cos\theta$ y est positif. La valeur de $\pi$ n’est jamais utilisée : le programme la fabrique.

from math import sqrt

s, n = 0.5, 6            # hexagone : sin(pi/6) = 1/2
for etape in range(30):
    print(f"{n:>12} cotes : {n * s:.15f}")
    s = sqrt((1 - sqrt(1 - s * s)) / 2)   # angle moitie, version naive
    n = 2 * n

Le cours s’arrête à quinze étapes ; poussons à trente, et regardons ce qui arrive.

           6 cotes : 3.000000000000000
          12 cotes : 3.105828541230250
          24 cotes : 3.132628613281237
          48 cotes : 3.139350203046872
          96 cotes : 3.141031950890530
         192 cotes : 3.141452472285344
        (...)
       12288 cotes : 3.141592618640789
       24576 cotes : 3.141592645321216
       49152 cotes : 3.141592645321216
       98304 cotes : 3.141592645321216
      196608 cotes : 3.141592645321216
      393216 cotes : 3.141593669849427
      786432 cotes : 3.141592303811738
     1572864 cotes : 3.141608696224804
     3145728 cotes : 3.141586839655041
     6291456 cotes : 3.141674265021758
    12582912 cotes : 3.141674265021758
    25165824 cotes : 3.143072740170040
    50331648 cotes : 3.159806164941135
   100663296 cotes : 3.181980515339464
   201326592 cotes : 3.354101966249685
   402653184 cotes : 4.242640687119286
   805306368 cotes : 6.000000000000000
  1610612736 cotes : 0.000000000000000
  3221225472 cotes : 0.000000000000000

Le piège des flottants, raconté

Trois phases, parfaitement lisibles.

La montée. Jusqu’à la treizième ligne, tout se passe comme les mathématiques le promettent : chaque doublement divise l’erreur par $4$, et la meilleure valeur atteinte est $3{,}141\,592\,645$, à $8{,}3 \cdot 10^{-9}$ de $\pi$. C’est la cinquième ligne qui mérite un salut au passage : $96$ côtés, $3{,}141\,03$, le polygone d’Archimède lui-même, au III^e siècle avant notre ère.

Le gel. Aux lignes $13$ à $16$, le programme affiche quatre fois exactement le même nombre. Il ne progresse plus, et il ne le dit pas.

L’effondrement. À partir de la dix-septième ligne, les valeurs se dégradent, d’abord discrètement, puis massivement : $3{,}143$, $3{,}16$, $3{,}18$, $3{,}35$, $4{,}24$, puis $6$, puis $0$. Le programme finit par annoncer que le demi-périmètre d’un polygone à trois milliards de côtés est nul.

Le coupable est dans la formule, et il est unique : la soustraction $1 - \sqrt{1 - s^2}$. Quand $s$ devient petit, $s^2$ devient minuscule, $\sqrt{1-s^2}$ devient un nombre très proche de $1$, et l’ordinateur retranche alors deux quantités presque égales. Les chiffres significatifs communs s’annulent, et il ne reste que les derniers, ceux que l’arrondi a déjà abîmés : c’est ce qu’on appelle une élimination catastrophique. À l’étape $24$, il ne survit plus que la première décimale ; à l’étape $28$, $s$ vaut environ $7 \cdot 10^{-9}$, si bien que $1 - s^2$ s’arrondit à $1$, que la différence vaut exactement $0$, et que tout s’éteint au pas suivant.

Ce n’est pas la méthode d’Archimède qui échoue : c’est son écriture. La preuve tient dans le programme suivant.

2. Archimède, la forme conjuguée

Reprenons la formule de l’angle moitié et multiplions haut et bas par la quantité conjuguée $1 + \cos\theta$, en écrivant $\cos\theta = \sqrt{1 - \sin^2\theta}$ :

$$\frac{1 - \cos\theta}{2} = \frac{(1-\cos\theta)(1+\cos\theta)}{2(1+\cos\theta)} = \frac{\sin^2\theta}{2(1+\cos\theta)}, \qquad \text{donc} \qquad \sin\frac{\theta}{2} = \frac{\sin\theta}{\sqrt{2\,(1 + \cos\theta)}}.$$

Mathématiquement, c’est la même formule. Numériquement, ce n’est plus du tout le même programme : la soustraction dangereuse a disparu, remplacée par une addition, et l’addition de deux nombres positifs ne perd jamais de chiffres significatifs.

from math import sqrt

s, n = 0.5, 6
for etape in range(30):
    print(f"{n:>12} cotes : {n * s:.15f}")
    c = sqrt(1 - s * s)                   # cos(theta), positif sur [0 ; pi/2]
    s = s / sqrt(2 * (1 + c))             # forme conjuguee : plus de soustraction
    n = 2 * n

Les dernières lignes de la sortie :

     3145728 cotes : 3.141592653589270
     6291456 cotes : 3.141592653589662
    12582912 cotes : 3.141592653589759
    25165824 cotes : 3.141592653589784
    50331648 cotes : 3.141592653589790
   100663296 cotes : 3.141592653589792
   201326592 cotes : 3.141592653589793
   402653184 cotes : 3.141592653589793
   805306368 cotes : 3.141592653589793
  1610612736 cotes : 3.141592653589793
  3221225472 cotes : 3.141592653589793

Le programme converge sagement jusqu’à la précision maximale d’un flottant, $3{,}141\,592\,653\,589\,793$, et s’y tient : les cinq dernières lignes sont identiques parce qu’il n’y a plus rien à gagner, non parce que quelque chose a cassé. Deux gels d’apparence semblable, deux causes opposées ; savoir les distinguer est une compétence à part entière.

Une note de contrôle, qui relie cette page à celle des sorties de l’Atelier : la ligne $393\,216$ côtés, la dix-septième, est exactement le polygone avec lequel Viète a calculé neuf décimales de $\pi$ en 1579. La version naïve s’y trompe dès la sixième décimale, elle n’en donne que cinq d’exactes ; la version conjuguée y donne dix décimales exactes.

3. Voir le théorème de Fourier, avec un curseur

Le script du cours trace trois sommes partielles du signal carré. En voici la version complète, avec un curseur qui fait varier $N$ de $1$ à $60$ et une ligne repère à la hauteur du sursaut. Elle demande matplotlib et ouvre une fenêtre interactive.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider


def somme_partielle(t, N):
    """Somme des N premieres sinusoides du signal carre."""
    return sum(4 / np.pi * np.sin((2 * k + 1) * t) / (2 * k + 1)
               for k in range(N))


t = np.linspace(-np.pi, np.pi, 2000)

figure, axes = plt.subplots(figsize=(9, 5))
plt.subplots_adjust(bottom=0.22)
axes.plot(t, np.sign(np.sin(t)), color="0.7", lw=1, label="signal carre")
courbe, = axes.plot(t, somme_partielle(t, 1), lw=1.8, label="somme partielle")
axes.axhline(1.17898, ls="--", lw=0.8, color="crimson")
axes.set_ylim(-1.45, 1.45)
axes.legend(loc="lower right")

curseur = Slider(plt.axes([0.15, 0.08, 0.7, 0.03]), "N", 1, 60,
                 valinit=1, valstep=1)


def redessiner(valeur):
    courbe.set_ydata(somme_partielle(t, int(valeur)))
    figure.canvas.draw_idle()


curseur.on_changed(redessiner)
plt.show()

Tirez le curseur lentement, et regardez le coin du créneau. Le plateau se raidit, les ondulations se resserrent, mais le petit pic accroché au saut ne baisse pas. Il se rapproche du bord, il s’affine, sa hauteur ne cède rien. Pour le mesurer sans le regarder, il suffit de remplacer la fenêtre par une boucle :

for N in (1, 3, 5, 15, 50, 200):
    u = np.linspace(1e-6, np.pi, 200001)
    print(N, round(float(somme_partielle(u, N).max()), 5))

Sortie :

1 1.27324
3 1.18836
5 1.18233
15 1.17935
50 1.17901
200 1.17898

C’est le phénomène de Gibbs de la pièce 58 des Coulisses, obtenu chez vous. La suite des maxima ne tend pas vers $1$, mais vers

$$\frac{2}{\pi}\int_0^{\pi}\frac{\sin t}{t}\,\mathrm{d}t = 1{,}178\,980\ldots,$$

soit un dépassement de $8{,}95\,\%$ du saut total, définitif. Et regardez l’intégrande : c’est le sinus cardinal, la fonction dont la limite en $0$ a débloqué toutes les dérivées du chapitre, et dont l’étude complète fait l’exercice 24 de l’Atelier. Tout se tient.

4. Le laboratoire du sauveteur

L’annexe B du devoir La lumière qui calcule propose un script à deux trous : la fonction $t$ de la question 4a, dont la formule est imprimée dans l’énoncé, et la mise à jour de la borne dans la dichotomie. Le voici complet et prêt à exécuter, avec les paramètres du devoir.

import math

a, b, d = 50, 30, 80
v1, v2 = 6, 1.5


def t(x):
    return math.sqrt(a*a + x*x)/v1 + math.sqrt((d - x)**2 + b*b)/v2


def tp(x):            # la derivee de Q5b
    return x/(v1*math.sqrt(a*a + x*x)) \
        - (d - x)/(v2*math.sqrt((d - x)**2 + b*b))


for x in [0, 20, 40, 50, 60, 70, 75, 80]:
    print(x, round(t(x), 2))

lo, hi = 0, d         # dichotomie : t' < 0 a gauche, > 0 a droite
for _ in range(50):
    m = (lo + hi)/2
    if tp(m) < 0:
        lo = m
    else:
        hi = m
print("x* =", round(m, 3), " t(x*) =", round(t(m), 4))

Une frontière que cette page ne franchit pas. La boucle d’affichage du haut redonne les huit valeurs de la table de la question 4b, et cette table est une question notée, sans corrigé : sa sortie n’est donc pas reproduite ici. Faites tourner le script, et confrontez ses lignes $x = 0$, $50$, $70$ et $80$ aux quatre valeurs déjà imprimées dans l’énoncé : si elles tombent juste, votre fonction t est correcte, et les quatre autres lignes le sont aussi. C’est le principe même d’une valeur de contrôle.

La dichotomie, en revanche, affiche exactement les deux nombres que l’annexe B du devoir imprime déjà, et les voici :

x* = 73.657  t(x*) = 35.2796

Deux mots sur la méthode, parce qu’elle est plus subtile qu’elle n’en a l’air. La dichotomie ne cherche pas le minimum de $t$ en comparant des valeurs de $t$ : elle cherche le zéro de $t'$, en comparant des signes. C’est possible parce que la partie II du devoir démontre que $t'$ est strictement croissante et change de signe entre les bornes, ce qui garantit un zéro et un seul. Sans ce théorème, l’encadrement lo, hi ne voudrait rien dire : rien n’interdirait à plusieurs zéros de se cacher dans l’intervalle, et la boucle en attraperait un au hasard. Le script ne remplace pas la démonstration, il en dépend.

Cinquante itérations divisent l’intervalle de départ, de largeur $80$, par $2^{50}$ : la largeur finale est de l’ordre de $7 \cdot 10^{-14}$ mètre, bien au-delà de ce que la précision des flottants sait garantir. Vingt itérations suffiraient largement pour les trois décimales affichées.

Le tracé, en supplément

Pour voir la vallée plutôt que la lire, quelques lignes suffisent. Le script trace la fonction durée sur tout l’intervalle.

import numpy as np
import matplotlib.pyplot as plt

a, b, d = 50, 30, 80
v1, v2 = 6, 1.5


def t(x):
    return np.sqrt(a*a + x*x)/v1 + np.sqrt((d - x)**2 + b*b)/v2


x = np.linspace(0, d, 1000)
plt.figure(figsize=(8, 4.5))
plt.plot(x, t(x), lw=2)
plt.xlabel("x (m) : point d'entree dans l'eau")
plt.ylabel("t(x) (s)")
plt.grid(alpha=0.3)
plt.show()

Vous retrouverez la figure du devoir, et surtout son asymétrie : la vallée descend en pente raide depuis la gauche et remonte très mollement vers la droite. La raison est physique, et elle vaut d’être méditée avant de se jeter à l’eau : une erreur de dix mètres vers la gauche se paie en mètres de nage, au tarif lent ; la même erreur vers la droite se paie en mètres de course. Le sauveteur qui hésite a intérêt à entrer dans l’eau trop loin plutôt que trop tôt.

En changeant les cinq paramètres du haut, le même script devient le banc d’essai de la question libre du devoir : inventez vos vitesses, vos distances, et regardez la brisure se déplacer.

Sources

• Cours F4, section 6, Le laboratoire : programme d’Archimède et note de marge « Le piège des flottants » ; programme des sommes de Fourier. Coulisses n°7, pièce 58, pour les valeurs du sursaut de Gibbs. Atelier F4, exercice 24, pour le sinus cardinal.
• Devoir F4 n°1, La lumière qui calcule, annexe B Le laboratoire du sauveteur, et partie II pour le théorème dont la dichotomie dépend.
• Tous les scripts de cette page ont été exécutés avant publication, sous Python 3.13 avec numpy et matplotlib ; les blocs de sortie sont recopiés de la console, sans retouche.