La course en direct : dichotomie contre Newton

Les Coulisses n°4 ont chronométré les deux méthodes et laissé ce tiroir ouvert. Voici la course, en scripts exécutables, sur l’équation qui ouvre et ferme le chapitre :

$x^5 - x - 1 = 0, \qquad \text{de racine } c = 1{,}16730397826141868\ldots$

Aucune formule ne donne $c$ : c’est le théorème d’Abel–Ruffini de l’ouverture du cours. Deux algorithmes se disputent le droit de l’écrire chiffre à chiffre.

Les deux concurrentes

La dichotomie part d’un encadrement $[a\,;b]$ où $f$ change de signe, coupe en deux, garde la moitié qui change encore de signe. Le milieu approche la racine à $\frac{b-a}{2^{n+1}}$ près après $n$ étapes : elle garantit une erreur, ce qui est rare et précieux.

La méthode de Newton–Raphson–Simpson part d’une estimation $x_0$ et itère

$\forall n \in \mathbb{N}, \qquad x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}.$

Elle ne garantit rien du tout, mais quand elle marche, elle double le nombre de décimales exactes à chaque tour. (Sur son vrai nom et sur la légende de la tangente, voir la page « Vrai ou légende ».)

Le script de la course

def f(x):
    return x**5 - x - 1

def f_prime(x):
    return 5 * x**4 - 1

def dichotomie(a, b, eps):
    """Renvoie le milieu et le nombre d'etapes, avec erreur garantie < eps."""
    n = 0
    while (b - a) / 2 > eps:
        m = (a + b) / 2
        if f(a) * f(m) <= 0:
            b = m
        else:
            a = m
        n = n + 1
    return (a + b) / 2, n

def newton(x0, eps):
    """Renvoie l'approximation et le nombre d'iterations."""
    x, n = x0, 0
    while abs(f(x)) > eps:
        x = x - f(x) / f_prime(x)
        n = n + 1
    return x, n

for eps in (1e-3, 1e-6, 1e-9, 1e-12):
    xd, nd = dichotomie(1.0, 2.0, eps)
    xn, nn = newton(1.2, eps)
    print("eps = %.0e : dichotomie %2d etapes -> %.12f | Newton %d iterations -> %.12f"
          % (eps, nd, xd, nn, xn))

Sortie :

eps = 1e-03 : dichotomie  9 etapes -> 1.166992187500 | Newton 2 iterations -> 1.167311020061
eps = 1e-06 : dichotomie 19 etapes -> 1.167304039001 | Newton 3 iterations -> 1.167303978357
eps = 1e-09 : dichotomie 29 etapes -> 1.167303978465 | Newton 3 iterations -> 1.167303978357
eps = 1e-12 : dichotomie 39 etapes -> 1.167303978261 | Newton 4 iterations -> 1.167303978261

Le verdict est sans appel : 39 étapes contre 4 itérations pour douze décimales.

Une précision de vocabulaire, qui vaut des points. Ici la condition d’arrêt porte sur l’erreur du milieu, $\frac{b-a}{2} \leqslant \varepsilon$, d’où les $19$ et $39$ étapes de la Coulisses. L’exercice 29 de l’Atelier arrête la boucle sur la longueur de l’intervalle, $b - a \leqslant \varepsilon$, et compte donc une étape de plus : $20$ et $40$. Les deux comptes sont exacts, ils ne mesurent simplement pas la même chose. Dites toujours laquelle des deux quantités vous majorez.

Le tableau de marche, tour par tour

Pour voir les deux régimes de convergence côte à côte.

racine = 1.1673039782614187        # valeur de reference, 17 chiffres

a, b = 1.0, 2.0
x = 1.2
print(" n | dichotomie : milieu | erreur  | Newton : x_n        | erreur")
for n in range(1, 7):
    m = (a + b) / 2
    if f(a) * f(m) <= 0:
        b = m
    else:
        a = m
    milieu = (a + b) / 2
    x = x - f(x) / f_prime(x)
    print("%2d | %19.15f | %.1e | %19.15f | %.1e"
          % (n, milieu, abs(milieu - racine), x, abs(x - racine)))

Sortie :

 n | dichotomie : milieu | erreur  | Newton : x_n        | erreur
 1 |   1.250000000000000 | 8.3e-02 |   1.169222886421862 | 1.9e-03
 2 |   1.125000000000000 | 4.2e-02 |   1.167311020061335 | 7.0e-06
 3 |   1.187500000000000 | 2.0e-02 |   1.167303978356634 | 9.5e-11
 4 |   1.156250000000000 | 1.1e-02 |   1.167303978261419 | 0.0e+00
 5 |   1.171875000000000 | 4.6e-03 |   1.167303978261419 | 0.0e+00
 6 |   1.164062500000000 | 3.2e-03 |   1.167303978261419 | 0.0e+00

Lisez les deux colonnes d’erreur : ce sont deux mondes.

• Dichotomie, convergence linéaire. L’erreur est divisée par 2 à chaque étape. Chaque décimale coûte $\log_2(10) \approx 3{,}32$ étapes, invariablement. C’est une progression géométrique de raison $\frac12$, exactement l’objet du chapitre A1.
• Newton, convergence quadratique. Les exposants de l’erreur se lisent $-3$, $-6$, $-11$ : ils doublent à chaque tour. Le nombre de décimales exactes double au lieu d’augmenter d’un tiers. À la quatrième itération, l’erreur est nulle à la précision de la machine : ce zéro n’est pas mathématique, il signifie seulement que les deux flottants les plus proches de $c$ sont épuisés.

L’origine de cette différence est dans les Coulisses : au voisinage du point fixe, la dichotomie multiplie l’erreur par $\frac12$, tandis que Newton l’élève au carré. C’est le même mécanisme que la pente nulle de l’itération de Héron, $f'(\sqrt{2}) = 0$, qui expliquait déjà sa vitesse deux mille ans plus tôt.

Trois pièges : quand la lièvre perd la course

La vitesse de Newton a un prix, et le script ci-dessous le fait payer. Chaque expérience est un contre-exemple authentique.

def essai(nom, g, g_prime, x0, tours=8):
    x = x0
    valeurs = []
    for _ in range(tours):
        d = g_prime(x)
        if d == 0:
            valeurs.append("derivee nulle : arret")
            break
        x = x - g(x) / d
        valeurs.append("%.6f" % x)
    print(nom, ":", " -> ".join(valeurs))

# Piege 1 : la tangente presque horizontale renvoie tres loin
essai("depart en 0.6 ", f, f_prime, 0.6)

# Piege 2 : un cycle qui ne converge jamais (x^3 - 2x + 2, depart 0)
essai("cycle 0 <-> 1 ", lambda x: x**3 - 2*x + 2,
                       lambda x: 3*x**2 - 2, 0.0)

# Piege 3 : convergence, mais vers l'autre racine
essai("depart en 0.4 ", lambda x: x**3 - x,
                       lambda x: 3*x**2 - 1, 0.4)

Sortie :

depart en 0.6  : -3.724545 -> -2.981696 -> -2.388871 -> -1.916727 -> -1.541404 -> -1.241686 -> -0.992738 -> -0.740821
cycle 0 <-> 1  : 1.000000 -> 0.000000 -> 1.000000 -> 0.000000 -> 1.000000 -> 0.000000 -> 1.000000 -> 0.000000
depart en 0.4  : -0.246154 -> 0.036457 -> -0.000097 -> 0.000000 -> 0.000000 -> 0.000000 -> 0.000000 -> 0.000000

• Piège 1, la tangente plate. En $x_0 = 0{,}6$, on a $f(0{,}6) = -1{,}52224$ et $f'(0{,}6) = 5 \times 0{,}6^4 - 1 = -0{,}352$ : la tangente est presque horizontale, elle coupe l’axe très loin, et le premier pas éjecte la suite en $-3{,}72$, à quatre unités de l’intervalle de départ. Elle rampe ensuite vers la racine, et après huit tours elle est encore négative, très loin derrière une dichotomie qui aurait déjà six chiffres binaires. Diviser par $f'(x_n)$ est le talon d’Achille : la méthode explose là où le dénominateur s’annule.
• Piège 2, le cycle. Sur $x^3 - 2x + 2$ en partant de $0$, la suite oscille indéfiniment entre $0$ et $1$. Elle ne diverge pas, elle ne converge pas : elle tourne en rond, pour toujours. Aucun test d’arrêt sur $|f(x)|$ ne s’y déclenche jamais, et la boucle while du premier script tournerait sans fin.
• Piège 3, la mauvaise racine. Sur $x^3 - x$, dont les racines sont $-1$, $0$ et $1$, un départ en $0{,}4$ franchit l’origine, repasse à $-0{,}246$, puis converge vers $0$. Rien n’assure que Newton attrape la racine que vous visiez : les bassins d’attraction des trois racines s’entrelacent, et le dessin de leur frontière est une fractale (Cayley l’a découvert dès 1879 en travaillant sur $z^3 - 1$ dans le plan complexe).

La dichotomie ne connaît aucun de ces trois accidents. Tant que $f$ est continue et change de signe aux bornes, le TVI garantit qu’une racine est dans l’intervalle, et l’intervalle rétrécit à coup sûr. Elle est lente et elle ne se trompe jamais.

Le vrai vainqueur : personne, ou plutôt les deux

C’est la morale annoncée par les Coulisses : robustesse contre vitesse est un faux duel. Les bibliothèques de calcul modernes ne choisissent pas, elles hybrident. La méthode de Brent, standard depuis les années 1970, maintient en permanence un encadrement à la dichotomie et tente à chaque tour un pas rapide (interpolation quadratique ou sécante) ; si le pas rapide sort de l’encadrement ou ne progresse pas assez, elle le rejette et fait une bissection.

Voici l’idée, dans sa version la plus simple : garder l’encadrement, essayer Newton, refuser s’il sort.

def newton_encadre(a, b, eps):
    """Newton quand son pas reste dans l'encadrement, bissection sinon."""
    x = (a + b) / 2
    n_newton, n_bissection = 0, 0
    for _ in range(200):
        d = f_prime(x)
        candidat = x - f(x) / d if d != 0 else None
        if candidat is not None and a <= candidat <= b:
            pas = candidat - x                 # pas rapide accepte
            x = candidat
            n_newton = n_newton + 1
        else:
            pas = (a + b) / 2 - x              # filet de securite
            x = (a + b) / 2
            n_bissection = n_bissection + 1
        if f(a) * f(x) <= 0:                   # l'encadrement reste valide
            b = x
        else:
            a = x
        if abs(pas) < eps:
            break
    return x, n_newton, n_bissection

for encadrement in [(1.0, 2.0), (0.0, 2.0), (-5.0, 5.0), (0.5, 2.0)]:
    x, nn, nb = newton_encadre(encadrement[0], encadrement[1], 1e-12)
    print("encadrement %-12s : x = %.15f  (%d pas de Newton, %d bissections)"
          % (str(encadrement), x, nn, nb))

Sortie :

encadrement (1.0, 2.0)   : x = 1.167303978261419  (7 pas de Newton, 0 bissections)
encadrement (0.0, 2.0)   : x = 1.167303978261419  (6 pas de Newton, 0 bissections)
encadrement (-5.0, 5.0)  : x = 1.167303978261419  (12 pas de Newton, 1 bissections)
encadrement (0.5, 2.0)   : x = 1.167303978261419  (5 pas de Newton, 0 bissections)

Une dizaine de pas au lieu de trente-neuf, et l’encadrement reste vrai du début à la fin. La troisième ligne est la plus parlante : sur l’intervalle très large $[-5\,;5]$, le premier pas de Newton part dans le décor, le filet se déclenche (une bissection), et la course reprend. Sans ce filet, le départ en $x_0 = 0$ aurait envoyé la suite exactement dans le piège n°1. On paie la sécurité en tenant un encadrement à jour, jamais en perdant la vitesse.

Deux remarques de lecture du script. Le test d’arrêt porte sur la taille du pas et non sur la largeur de l’encadrement : après un bon pas de Newton, l’encadrement contient toujours un bord resté loin, et attendre qu’il se referme reviendrait à payer les trente-neuf bissections. Et la mise à jour des bornes est faite quel que soit le type de pas : c’est elle qui garantit qu’à tout instant une racine reste piégée entre $a$ et $b$.

C’est ce qui tourne derrière scipy.optimize.brentq en Python, derrière fzero en MATLAB, et derrière la touche « résoudre » de votre calculatrice.

Sources

• Pièce 34 des Coulisses n°4, Deux façons d’attraper une racine, et le paragraphe 6 du cours F1.
• Richard Brent, Algorithms for Minimization without Derivatives, 1973, pour la méthode hybride.
• Arthur Cayley (1879) pour les bassins d’attraction de Newton sur $z^3 - 1$.
• Tous les scripts de cette page ont été exécutés ; les sorties sont recopiées telles quelles.