Équations non linéaires

Objectifs du chapitre

À la fin de ce chapitre, vous serez capable de:

• localiser une racine d'une équation non linéaire par le théorème des valeurs intermédiaires, et dire ce que ce théorème garantit — et ce qu'il ne garantit pas;
• programmer la bissection, majorer son erreur par $(b-a)/2^{k+1}$ et en déduire a priori le nombre d'itérations nécessaires pour une tolérance donnée;
• reformuler une équation en $x=g(x)$, tester la condition de contraction $|g^{\prime }|<1$ et prévoir si l'itération convergera ou divergera;
• écrire la méthode de Newton, démontrer sa convergence quadratique pour une racine simple, et reconnaître les trois situations où elle échoue;
• diagnostiquer une racine multiple à la perte de la convergence quadratique, et restaurer celle-ci par le Newton modifié;
• étendre le point fixe et Newton aux systèmes $\mathbf{F}(\mathbf{x})=0$: condition de contraction par le rayon spectral de la jacobienne, un système linéaire par pas de Newton, et jacobienne approchée quand elle coûte trop cher;
• choisir un critère d'arrêt en sachant que ni le résidu $|f(x_k)|$ ni l'incrément $|x_{k+1}-x_k|$ ne majorent l'erreur $|x_k-x^{\ast }|$.

Une équation, et aucune formule

Résoudre $a{x}^2+bx+c=0$ demande une formule; résoudre $x^3-x-1=0$ en demande une autre, due à Cardan, déjà pénible; résoudre une équation polynomiale de degré cinq n'en demande aucune, car le théorème d'Abel–Ruffini établit qu'il n'en existe pas par radicaux. Et dès que l'équation cesse d'être polynomiale — , , ou l'équation de Kepler qui donne la position d'une planète sur son orbite — la question d'une formule ne se pose même plus.

L'analyse numérique répond autrement: elle renonce à la formule et construit une suite $x_0,x_1,x_2,\dots$ qui converge vers la solution, en s'arrêtant quand la précision demandée est atteinte. Ce chapitre étudie quatre façons de construire cette suite, et surtout la question qui les sépare: à quelle vitesse?

Toutes les méthodes de ce chapitre sont illustrées sur la même équation, celle qui sert de fil rouge au cours:

Ses valeurs aux bornes valent $f(1)=-1$ et $f(2)=5$, et sa racine réelle est le nombre plastique

Nous l'avons calculée en arithmétique décimale à soixante chiffres, puis arrondie à la double précision: c'est la valeur de référence contre laquelle toutes les erreurs de ce chapitre sont mesurées. Deux constantes reviendront:

Existence et localisation d'une racine

Avant d'itérer, il faut savoir qu'il y a quelque chose vers quoi converger, et savoir à peu près où.

Ce théorème appartient à l'analyse réelle et sa démonstration — par dichotomie, précisément, ou par la borne supérieure de l'ensemble où $f$ est négative — figure au cours d'Analyse I, chapitre 5. Nous l'admettons ici, mais il faut remarquer que la démonstration usuelle de ce théorème est l'algorithme de la section suivante: la bissection n'est pas une astuce de calcul greffée sur le théorème, c'est sa preuve, exécutée.

L'unicité, elle, se déduit de la monotonie. Si $f$ est continue sur $[a,b]$, dérivable sur $]a,b[$, avec $f(a)f(b)<0$ et $f^{\prime }$ de signe constant, alors la racine est unique: deux racines distinctes $x_1<x_2$ donneraient, par le théorème de Rolle (Analyse I, chapitre 8), un point où s'annule, ce qui contredit le signe constant.

La bissection

L'idée tient en une phrase: si $f$ change de signe sur $[a,b]$, on coupe l'intervalle en deux, on regarde dans quelle moitié le changement de signe subsiste, et l'on recommence.

def bissection(f, a, b, tol=1e-10):
    """Racine de f sur [a, b], en supposant f(a) f(b) < 0."""
    fa = f(a)
    while (b - a) / 2 > tol:
        c = (a + b) / 2
        fc = f(c)
        if fa * fc <= 0:      # le changement de signe est dans [a, c]
            b = c
        else:                 # il est dans [c, b]

Chaque itération coûte une seule évaluation de $f$ — c'est le point de la variable fa, qui mémorise $f(a)$ au lieu de le recalculer. Le nombre total d'évaluations est donc égal au nombre d'itérations, plus une.

Démonstration. Montrons d'abord par récurrence que, pour tout $k$, l'intervalle $[a_k,b_k]$ vérifie les deux propriétés

Pour $k=0$ c'est l'hypothèse et une évidence. Supposons-les au rang $k$. La règle (2.4) choisit $[a_k,c_k]$ lorsque $f(a_k)f(c_k)<0$, et sinon; dans le premier cas la propriété de signe est immédiate, dans le second on a , donc est de même signe que (ou nul), et comme on en tire . Dans les deux cas la longueur est exactement divisée par deux, puisque est le milieu. La récurrence est établie.

Les suites $(a_k)$ et $(b_k)$ sont respectivement croissante et décroissante, et leur écart $(b-a)2^{-k}$ tend vers $0$: elles sont (Analyse I, chapitre 3) et convergent vers une limite commune . Par continuité de , et , donc

Or chaque terme de cette suite est $\le 0$, donc sa limite l'est aussi: $f(x^{\ast }{)}^2\le 0$, ce qui force $f(x^{\ast })=0$. La limite est bien une racine.

Il reste la majoration. Le milieu $c_k$ et la limite $x^{\ast }$ appartiennent tous deux à $[a_k,b_k]$, et $c_k$ en est le centre; leur distance ne dépasse donc pas la demi-longueur:

ce qui est (2.5). Enfin, $(b-a)2^{-(k+1)}\le \epsilon$ équivaut à $2^{k+1}\ge (b-a)/\epsilon$, c'est-à-dire à (2.6). $\square$

Deux remarques sur cette borne, qui en font la plus honnête du chapitre. D'abord elle est a priori: elle ne dépend ni de $f$, ni de la racine, ni du hasard des itérés. Sur $[1,2]$ et pour $\epsilon =10^{-6}$, elle donne

et ces vingt itérations suffiront quelle que soit la fonction. Aucune autre méthode de ce chapitre ne permet ce genre de promesse. Ensuite elle est une majoration et rien d'autre: l'erreur réelle oscille sous la borne sans la suivre, et il lui arrive d'augmenter d'une itération à la suivante.

La bissection a une vertu et un défaut, et les deux sont extrêmes. Sa vertu: elle converge toujours, avec une borne connue d'avance, sans hypothèse sur $f$ autre que la continuité — ni dérivabilité, ni proximité du point de départ. Son défaut: elle gagne exactement ${\log }_{10}2=0,301$ chiffre décimal par itération, soit un chiffre toutes les $3,32$ itérations, et cette vitesse ne s'améliore jamais, si régulière que soit la fonction. Elle ignore complètement la valeur de $f$: seul son signe est utilisé, ce qui revient à jeter la quasi-totalité de l'information que chaque évaluation a coûté. Toutes les méthodes qui suivent exploitent cette information, et vont plus vite — au prix de garanties nettement plus faibles.

Les itérations de point fixe

La bissection n'utilise que le signe de $f$. Pour aller plus vite, il faut se servir des valeurs, et la façon la plus simple de le faire consiste à réécrire l'équation sous une forme qui se lise elle-même comme une recette de calcul.

Il y a toujours plusieurs façons de construire une telle $g$, et c'est là que tout se joue. Pour l'équation témoin $x^3-x-1=0$, deux réécritures s'offrent immédiatement:

Elles sont mathématiquement équivalentes: les deux équations $x=g_1(x)$ et $x=g_2(x)$ ont exactement la même solution réelle, $x^{\ast }=1,324717957244746$. Numériquement, elles n'ont rien à voir l'une avec l'autre.

La figure 2.2 rend visible le mécanisme. L'itération se lit géométriquement ainsi: partant de $(x_k,x_k)$ sur la diagonale, on monte (ou descend) verticalement jusqu'à la courbe $y=g(x)$, ce qui donne l'ordonnée $g(x_k)=x_{k+1}$; puis on se déplace horizontalement jusqu'à la diagonale, ce qui reporte cette ordonnée en abscisse. Si , la courbe est plus plate que la diagonale et les marches rétrécissent; si , elle est plus raide et elles s'agrandissent. Le théorème suivant transforme cette image en démonstration.

Démonstration. Existence. Posons $h(x)=g(x)-x$, continue sur $[a,b]$. L'hypothèse de stabilité donne $g(a)\ge a$ et $g(b)\le b$, donc $h(a)\ge 0$ et . Le théorème 2.1 (ou son cas limite, si l'une des deux inégalités est une égalité) fournit un avec , c'est-à-dire .

Unicité. Supposons deux points fixes $x^{\ast }\ne y^{\ast }$ dans $[a,b]$. Le théorème des accroissements finis (Analyse I, chapitre 8) donne un $\xi$ strictement entre eux tel que

Comme $|x^{\ast }-y^{\ast }|>0$, on peut diviser: $1\le L$, ce qui contredit $L<1$. Le point fixe est donc unique.

Convergence. La stabilité assure que toute la suite reste dans $[a,b]$, par récurrence immédiate. Appliquons les accroissements finis entre $x_k$ et $x^{\ast }$: il existe ${\xi }_k$ entre eux avec

d'où $|x_{k+1}-x^{\ast }|\le L|x_k-x^{\ast }|$ et, par récurrence, la première inégalité de (2.9). Comme , et .

La seconde majoration, qui est celle qu'on peut réellement calculer puisqu'elle ne fait pas intervenir $x^{\ast }$, s'obtient en télescopant. Pour $j\ge 1$, le même argument donne $|x_{j+1}-x_j|\le L|x_j-x_{j-1}|$, donc . Pour , l'inégalité triangulaire donne alors

la dernière étape étant la somme d'une série géométrique de raison $L<1$. En faisant tendre $m$ vers l'infini à $k$ fixé, $x_m\to x^{\ast }$ et l'on obtient la seconde inégalité de (2.9). $\square$

En pratique, on ne vérifie presque jamais les hypothèses du théorème 2.3 sur un intervalle explicite: on se contente du critère local, qui en découle par continuité.

Démonstration. Supposons $|g^{\prime }(x^{\ast })|=\lambda <1$ et posons $L=(1+\lambda )/2$, de sorte que $\lambda <L<1$. Par continuité de , il existe tel que sur . Cet intervalle est stable: pour , les accroissements finis donnent , donc . Les deux hypothèses du théorème 2.3 sont réunies sur , d'où la convergence.

Si $|g^{\prime }(x^{\ast })|>1$, choisissons $\mu >1$ et $\delta >0$ tels que $|g^{\prime }(x)|\ge \mu$ sur . Tant que , la relation (2.10) donne : l'erreur est par un facteur au moins . Si , l'erreur croît géométriquement et la suite quitte nécessairement au bout d'un nombre fini d'étapes — elle ne peut donc pas converger vers , puisqu'une suite convergente finirait par rester dans .

C'est très exactement ce que le tableau de l'exemple 2.3 montre: $|g_1^{\prime }(x^{\ast })|=0,19<1$ et la suite converge; $|g_2^{\prime }(x^{\ast })|=5,26>1$ et elle diverge, y compris en partant à de la cible. C'est l'illustration la plus économique de tout le chapitre, et elle mérite d'être gardée en tête: le comportement d'une méthode itérative n'est pas une propriété du problème, c'est une propriété de la reformulation choisie.

La lecture pratique de (2.11) est une affaire de chiffres. Si $|e_k|\approx 10^{-d_k}$, alors l'ordre $p$ donne $d_{k+1}\approx p{d}_k-{\log }_{10}C$: Ordre 1 avec un taux : on ajoute chiffres par tour, un nombre constant. Ordre 2: on le nombre de chiffres. C'est la différence entre additionner et multiplier, et sur quatre itérations elle est spectaculaire.

Une itération de point fixe avec $g^{\prime }(x^{\ast })\ne 0$ est d'ordre exactement 1, par (2.10): le rapport $|e_{k+1}|/|e_k|$ tend vers , ce qui est précisément la colonne de droite du premier tableau de l'exemple 2.3. Toute la question devient alors: La réponse est oui, et elle porte un nom.

La méthode de Newton

Revenons à $f(x)=0$ et changeons de point de vue. Au lieu de chercher une réécriture astucieuse, remplaçons $f$ par ce que nous savons le mieux d'elle près de $x_k$: sa tangente.

Si $f$ est dérivable, la formule de Taylor au premier ordre (Analyse I, chapitre 9) donne

En négligeant le reste et en cherchant le zéro de ce qui reste — un polynôme du premier degré, dont on sait résoudre l'équation — on obtient le point où la tangente coupe l'axe:

def newton(f, fp, x0, tol=1e-12, kmax=50):
    """Renvoie la liste des iteres; s'arrete quand le residu passe sous tol."""
    x = x0
    suite = [x]
    for k in range(kmax):
        d = fp(x)
        if d == 0.0:
            return suite          # tangente horizontale: on ne peut pas continuer
        x = x - f(x) / d
        suite.append(x)

Deux évaluations par itération — $f(x_k)$ et $f^{\prime }(x_k)$ — contre une seule pour la bissection et le point fixe. Nous verrons qu'elles sont largement rentabilisées, mais c'est une dépense réelle, et le fait que $f^{\prime }$ doive être disponible en est une autre: dans bien des applications, est le résultat d'une simulation et sa dérivée n'existe que sous forme approchée. C'est précisément ce qui motivera la méthode de la sécante.

Démonstration. Comme $f^{\prime }(x^{\ast })\ne 0$ et $f^{\prime }$ est continue, il existe $\rho >0$ et $m>0$ tels que sur ; l'itération y est donc bien définie. Posons aussi , fini par continuité de sur le compact .

Soit $x_k\in J$. La formule de Taylor avec reste de Lagrange (Analyse I, chapitre 9), appliquée à $f$ entre $x_k$ et $x^{\ast }$, fournit un ${\xi }_k$ strictement compris entre et tel que

C'est ici que tout se joue: l'équation (2.14) est exacte, le reste de Lagrange n'est pas une approximation. Divisons-la par $f^{\prime }(x_k)\ne 0$ et réorganisons:

soit, en reconnaissant $x_{k+1}=x_k-f(x_k)/f^{\prime }(x_k)$,

La relation (2.15) est l'énoncé du théorème sous sa forme brute: l'erreur nouvelle est proportionnelle au carré de l'ancienne. Il reste à en tirer la convergence.

Posons $K=M/(2m)$, de sorte que $|e_{k+1}|\le K|e_k{|}^2$ pour $x_k\in J$. Choisissons

Si $x_k\in I_{\delta }$, alors $|e_k|\le \delta$ et

donc $x_{k+1}\in I_{\delta }$ également. Par récurrence, toute la suite reste dans $I_{\delta }$ et $|e_k|\le 2^{-k}|e_0|\to 0$: la convergence est acquise. (Remarquez que cette première majoration n'est que linéaire — elle sert uniquement à établir la suite converge; la vitesse réelle vient ensuite.)

Enfin, puisque $x_k\to x^{\ast }$ et que ${\xi }_k$ est coincé entre $x_k$ et $x^{\ast }$, on a aussi . En divisant (2.15) par et en passant à la limite, la continuité de et de donne

ce qui est (2.13) au module près. $\square$

Là où Newton échoue

Le théorème 2.5 est local: il affirme l'existence d'un $\delta$, sans le donner, et ne dit rien du comportement hors de $I_{\delta }$. Voici les trois manières dont les choses tournent mal, chacune sur un exemple réel et calculé.

Échec 1: la tangente horizontale. Sur l'équation témoin, $f^{\prime }(x)=3x^2-1$ s'annule en $\pm 1/\sqrt{3}=\pm 0,5773503$. Partons de , où est petite et ne l'est pas. La tangente y est presque horizontale et la projette très loin:

La suite continue ensuite $-3,364865$, $-2,280955$, $-1,556277$, $-1,043505$, remontant lentement depuis là où elle a été expédiée. Elle finira par revenir — $f$ n'a qu'une racine réelle — mais elle aura dépensé une dizaine d'itérations à réparer un seul pas.

Le cas limite est plus instructif encore. Partons exactement de $x_0=1/\sqrt{3}$ tel que la machine le représente. En arithmétique exacte, $f^{\prime }(x_0)=0$ et l'itération n'est pas définie. En double précision, est évalué à — pas zéro, du bruit d'arrondi — et la division donne

Aucune exception, aucun message: le programme continue tranquillement à seize ordres de grandeur de la racine. C'est le chapitre 1 qui parle: le test if fp(x) == 0.0: ne protège de rien, car la dérivée ne vaut pratiquement jamais zéro exactement. Un vrai garde-fou compare $|f^{\prime }(x_k)|$ à un seuil, ou surveille la taille du pas.

Échec 2: le cycle. Prenons $f(x)=x^3-2x+2$, dont l'unique racine réelle vaut $-1,769292354238631$, et partons de $x_0=0$:

La suite vaut $0,1,0,1,0,1,\dots$ indéfiniment. Elle ne diverge pas, elle ne converge pas, et surtout chaque itération est parfaitement licite: la dérivée ne s'annule jamais, le pas ne devient pas grand, aucun garde-fou naïf ne se déclenche. Seul un plafond d'itérations arrête le programme. Notez que le cycle est stable: en partant de $x_0=0,01$ on tourne un moment autour du même cycle avant de s'en échapper, ce qui rend le phénomène robuste et non anecdotique.

Échec 3: le mauvais point de départ, tout simplement. Prenons $f(x)=\arctan x$, dont l'unique racine est $0$. L'itération s'écrit $x_{k+1}=x_k-(1+x_k^2)\arctan x_k$. Depuis :

convergence quadratique impeccable. Depuis $x_0=1,5$:

divergence, en changeant de signe à chaque pas. La frontière entre les deux comportements est le point $\bar{x}=1,3917452$ solution de $(1+x^2)\arctan x=2x$: en deçà l'itération converge, au-delà elle s'éloigne. Un dixième d'écart sur $x_0$ sépare ici la réussite de la catastrophe, sur une fonction pourtant analytique, monotone et à racine simple.

Racines multiples

Le théorème 2.5 exige $f^{\prime }(x^{\ast })\ne 0$. Quand cette hypothèse tombe, la conclusion tombe avec elle — et pas d'un peu.

Démonstration. Écrivons $f(x)=(x-x^{\ast }{)}^{m}h(x)$ avec $h$ de classe ${\mathcal{C}}^1$ et $h(x^{\ast })\ne 0$. Alors

de sorte que, pour $x\ne x^{\ast }$ voisin,

La fonction d'itération de Newton, $g(x)=x-f(x)/f^{\prime }(x)$, vérifie donc

En divisant par $x-x^{\ast }$ et en faisant tendre $x$ vers $x^{\ast }$, le crochet tend vers $1-h(x^{\ast })/(mh(x^{\ast }))=1-1/m$, ce qui est (2.17) — et qui n'est autre que . Pour ce nombre est dans : strictement inférieur à 1, donc le point fixe reste attractif, mais non nul, donc la convergence n'est plus que linéaire. Pour on retrouve et le théorème 2.5.

Pour le Newton modifié, la fonction d'itération devient $G(x)=x-mf(x)/f^{\prime }(x)$, et le même calcul donne

Donc $G^{\prime }(x^{\ast })=0$, et le terme linéaire de l'erreur disparaît exactement comme dans le cas simple: la convergence redevient au moins quadratique. $\square$

Sur l'exemple 2.5, avec $m=2$, l'itération (2.18) depuis $x_0=2$ donne

soit des erreurs $1$, $1,11\cdot 10^{-1}$, $1,95\cdot 10^{-3}$, $6,33\cdot 10^{-7}$, : les rapports valent , , — la convergence quadratique est de retour, avec sa constante qui se stabilise. La dernière ligne est déjà sous l'effet de l'arrondi et son rapport n'a plus de sens: nous ne l'imprimons pas.

La méthode de la sécante

Newton demande $f^{\prime }$. Or il arrive très souvent qu'on ne l'ait pas: $f$ peut être le résultat d'une simulation, d'une table, d'un solveur imbriqué. L'idée de la sécante est de remplacer la tangente par la droite qui passe par les deux derniers points calculés.

Chaque itération ne demande qu'une seule nouvelle évaluation de $f$ — celle en $x_k$; $f(x_{k-1})$ a été calculée au tour précédent et se recycle. C'est la moitié du coût de Newton, et cela suffit à changer le verdict de la comparaison.

Nous admettons ici la convergence — sa démonstration reprend, avec deux points au lieu d'un, l'argument de contraction du théorème 2.5 et n'apporte rien de neuf — mais l'ordre, lui, se démontre en quelques lignes et fait l'objet de l'exercice 2.6, où vous établirez d'abord la relation d'erreur

analogue de (2.15), puis en déduirez $\phi$ comme racine positive de $p^2=p+1$. C'est l'un des rares endroits d'un cours d'analyse numérique où le nombre d'or apparaît sans être invité, et la raison en est transparente: (2.22) fait intervenir les deux erreurs précédentes, donc l'exposant suit une récurrence de Fibonacci.

Une variante mérite d'être nommée, parce qu'elle est souvent confondue avec la sécante: la méthode de la fausse position (regula falsi) utilise la même formule (2.20), mais choisit les deux points de façon à conserver un encadrement de la racine — comme la bissection. Elle est donc toujours convergente, ce qui est un vrai gain; en contrepartie, l'une des deux bornes reste souvent figée et la convergence retombe à l'ordre 1. C'est le compromis inverse de celui de la sécante, et il explique pourquoi les bibliothèques modernes préfèrent des hybrides (Brent, Illinois) qui gardent l'encadrement et l'ordre superlinéaire.

Les critères d'arrêt

Nous savons construire des suites convergentes. Reste la question la plus pratique du chapitre, et celle où l'on se trompe le plus: quand s'arrêter? Trois candidats se présentent, et aucun n'est ce qu'on croit.

Candidat 1: le résidu, $|f(x_k)|\le \epsilon$. Il est tentant car il ne demande rien d'autre qu'une évaluation déjà faite. Mais le résidu et l'erreur sont reliés par la pente:

Le facteur $1/|f^{\prime }(x^{\ast })|$ est le conditionnement du problème «trouver la racine», au sens exact du chapitre 1: il dit de combien une perturbation de $f$ déplace la racine. Si $f^{\prime }(x^{\ast })$ est grande, le résidu est pessimiste et tout va bien. S'il est petit — racine multiple ou presque multiple — le résidu ment dans le mauvais sens. Sur $f(x)=(x-1{)}^7$, un résidu de correspond à une erreur de : on est à dix pour cent de la racine avec un résidu de sept chiffres.

Candidat 2: l'incrément, $|x_{k+1}-x_k|\le \epsilon$. C'est le critère que tout le monde écrit, et c'est celui contre lequel ce chapitre doit vous mettre en garde.

Chiffrons-le sur l'équation témoin, avec une itération parfaitement légitime: la relaxation

dont la fonction d'itération $g(x)=x-\omega f(x)$ a pour dérivée $g^{\prime }(x^{\ast })=1-\omega f^{\prime }(x^{\ast })=1-0,004265=0,995735$. Le facteur (2.24) vaut donc . Voici ce que le programme observe:

Le rapport converge vers $234,5$, la valeur prédite par (2.24). Concrètement: un programme qui s'arrête sur $|{\Delta }_k|<10^{-8}$ le fait ici à l'itération $k=2590$, avec $x_k=1,324720297448863$ et un pas de — et une , deux cent trente-quatre fois plus grande que ce que le critère laissait croire. Le programme annonce huit chiffres et en a cinq.

Le cas extrême est instructif: la suite $x_k={\sum }_{i=1}^{k}1/i$ a des pas $1/k$ qui tendent vers zéro, et pourtant elle diverge. Un pas qui tend vers zéro ne prouve strictement rien sur la convergence. L'incrément n'est un estimateur honnête de l'erreur que lorsqu'on sait, par ailleurs, que $g^{\prime }(x^{\ast })$ est loin de 1 — ce qui est le cas pour Newton et la sécante, et ce qui ne l'est pas pour une itération lente.

Candidat 3: l'encadrement. Seule la bissection (et ses hybrides) fournit une majoration garantie de l'erreur, la demi-longueur de l'intervalle courant. C'est le seul critère d'arrêt de ce chapitre qui soit une véritable borne et non une estimation, et c'est la raison pour laquelle on paie sa lenteur.

En pratique, on combine, et l'on écrit toujours un critère mixte — relatif pour les grandes valeurs, absolu près de zéro, exactement comme la comparaison à tolérance du chapitre 1:

def arret(x_new, x_old, f_new, tol_x=1e-10, tol_f=1e-12):
    """Vrai si l'on peut s'arreter. Jamais de test d'egalite entre flottants."""
    pas_petit = abs(x_new - x_old) <= tol_x * max(1.0, abs(x_new))
    residu_petit = abs(f_new) <= tol_f
    return pas_petit and residu_petit

et l'on y ajoute toujours le plafond kmax, avec un code de retour qui distingue «convergé» de «abandonné après kmax itérations». Un solveur qui ne peut pas signaler son échec est un solveur qui ment.

Systèmes d'équations non linéaires

Tout ce qui précède porte sur une équation à une inconnue. La pratique en demande presque toujours davantage: l'équilibre d'une structure articulée, le point de fonctionnement d'un circuit à diodes, la composition d'un mélange réactif à l'équilibre, et — nous le verrons au chapitre 10 — chaque pas d'un schéma implicite pour un système d'équations différentielles. On cherche alors ${\mathbf{x}}^{\ast }\in {\mathbb{R}}^n$ tel que

soit $n$ équations scalaires couplées. La bissection ne survit pas au passage: il n'existe pas de «changement de signe» d'une fonction vectorielle, ni d'intervalle à couper en deux. Les deux autres idées du chapitre, en revanche, se généralisent presque mot pour mot — le point fixe et Newton — à condition de remplacer la dérivée par une matrice.

Le système témoin de cette section est l'intersection d'un cercle et d'une hyperbole:

Il a quatre solutions, symétriques deux à deux, et celle du premier quadrant sous la diagonale se calcule exactement: en élevant $xy=1$ au carré et en substituant, $x^4-4x^2+1=0$, d'où $x^{\ast 2}=2+\sqrt{3}$ et , c'est-à-dire

Avoir une solution exacte est un luxe que la pratique n'offre pas, et que nous prenons ici pour la même raison qu'au début du chapitre: mesurer de vraies erreurs. Le jacobien en ce point vaut $\det J_{\mathbf{F}}({\mathbf{x}}^{\ast })=2(x^{\ast 2}-y^{\ast 2})=4\sqrt{3}\ne 0$: la racine est simple.

Le point fixe dans ${\mathbb{R}}^n$

On réécrit $\mathbf{F}(\mathbf{x})=0$ sous la forme $\mathbf{x}=\mathbf{G}(\mathbf{x})$ et l'on itère ${\mathbf{x}}^{(k+1)}=\mathbf{G}({\mathbf{x}}^{(k)})$, en notant désormais l'indice d'itération en exposant, comme le fera le chapitre 4. La valeur absolue devient une norme, et la démonstration du théorème 2.3 se transporte sans changement, à un détail près qu'il faut voir.

Démonstration. Le télescopage de la démonstration du théorème 2.3 s'écrit à l'identique avec une norme: $\Vert {\mathbf{x}}^{(j+1)}-{\mathbf{x}}^{(j)}\Vert \le L^j\Vert {\mathbf{x}}^{(1)}-{\mathbf{x}}^{(0)}\Vert$, puis, pour ,

La suite est donc de Cauchy. Ici l'argument change de nature: sur un intervalle, l'existence d'un point fixe venait du théorème des valeurs intermédiaires, qui n'a pas d'analogue dans ${\mathbb{R}}^n$. On utilise à la place la complétude de ${\mathbb{R}}^n$ (Analyse II, chapitre 2): toute suite de Cauchy converge, vers un ${\mathbf{x}}^{\ast }$ qui appartient à $D$ puisque $D$ est fermé. La contraction rend $\mathbf{G}$ continue, donc en passant à la limite dans ${\mathbf{x}}^{(k+1)}=\mathbf{G}({\mathbf{x}}^{(k)})$ on obtient . L'unicité et la première inégalité de (2.28) viennent de , la seconde du passage à la limite ci-dessus.

Pour le dernier point, voici le détail annoncé. Il n'y a pas de théorème des accroissements finis avec égalité pour une fonction vectorielle — il n'existe en général aucun point $\mathbf{\xi }$ tel que $\mathbf{G}(\mathbf{x})-\mathbf{G}(\mathbf{y})=J_{\mathbf{G}}(\mathbf{\xi })(\mathbf{x}-\mathbf{y})$. On le remplace par la forme intégrale, valable parce que le segment $[\mathbf{y},\mathbf{x}]$ est contenu dans convexe:

d'où $\Vert \mathbf{G}(\mathbf{x})-\mathbf{G}(\mathbf{y})\Vert \le {\sup }_D\Vert J_{\mathbf{G}}\Vert \Vert \mathbf{x}-\mathbf{y}\Vert \le L\Vert \mathbf{x}-\mathbf{y}\Vert$: c'est l' des accroissements finis, et elle suffit.

Comme en dimension un, on ne vérifie presque jamais ces hypothèses sur un domaine explicite; on se sert du critère local, où la dérivée $g^{\prime }(x^{\ast })$ devient la jacobienne et sa valeur absolue le rayon spectral $\rho$, plus grand module des valeurs propres (annexe A).

Démonstration. Posons $\rho =\rho (J_{\mathbf{G}}({\mathbf{x}}^{\ast }))<1$ et $L=(1+\rho )/2$. Le théorème de Householder — énoncé au chapitre 4 comme partie du théorème 4.5 et admis à cet endroit — fournit une norme induite telle que $\Vert J_{\mathbf{G}}({\mathbf{x}}^{\ast })\Vert <L$. Par continuité de , on a encore sur une boule fermée de centre et de rayon assez petit, pour cette norme. est fermée et convexe, et elle est stable: . Le théorème 2.8 s'applique. Le taux asymptotique , plutôt que , vient de la formule de Gelfand de l'annexe A, appliquée à la récurrence linéarisée ; nous l'admettons.

Le passage de la norme au rayon spectral n'est pas une coquetterie. Une jacobienne peut avoir toutes ses normes usuelles supérieures à 1 et un rayon spectral inférieur à 1: l'itération converge alors, alors qu'aucune des normes ne permettait de le prévoir. C'est exactement la situation que le chapitre 4 retrouvera pour les méthodes itératives linéaires, dont l'itération est un point fixe affine ${\mathbf{x}}^{(k+1)}=B{\mathbf{x}}^{(k)}+\mathbf{f}$ de jacobienne constante $B$.

La méthode de Newton pour les systèmes

L'idée de la tangente se transporte telle quelle: près de ${\mathbf{x}}^{(k)}$, on remplace $\mathbf{F}$ par son approximation affine $\mathbf{F}({\mathbf{x}}^{(k)})+J_{\mathbf{F}}({\mathbf{x}}^{(k)})(\mathbf{x}-{\mathbf{x}}^{(k)})$ et l'on annule celle-ci. Il n'y a plus de division par , mais un à résoudre.

On écrit parfois ${\mathbf{x}}^{(k+1)}={\mathbf{x}}^{(k)}-J_{\mathbf{F}}({\mathbf{x}}^{(k)}{)}^{-1}\mathbf{F}({\mathbf{x}}^{(k)})$, et c'est très bien sur le papier. Dans un programme, : on résout le système (2.29) par l'élimination de Gauss avec pivot partiel du chapitre 3, ce qui coûte opérations, contre pour calculer l'inverse puis pour le multiplier par . Voici la méthode, avec une petite élimination écrite pour l'occasion:

def resout(A, b):
    """Elimination de Gauss avec pivot partiel (chapitre 3), puis remontee."""
    n = len(b)
    M = [A[i][:] + [b[i]] for i in range(n)]
    for k in range(n):
        p = max(range(k, n), key=lambda i: abs(M[i][k]))
        M[k], M[p] = M[p], M[k]
        for i in range(k +

Le coût d'une itération se lit sur le code: une évaluation de $\mathbf{F}$ ($n$ fonctions scalaires), une évaluation de $J_{\mathbf{F}}$ ($n^2$ dérivées partielles), et une résolution en $\frac{2}{3}{n}^3+O(n^2)$ opérations. Pour grand, c'est la factorisation qui domine, et l'on verra plus bas comment l'économiser.

Démonstration. Nous admettons un lemme de perturbation, conséquence du lemme de Neumann que le chapitre 4 invoque pour (4.16): si $\Vert \mathbf{x}-{\mathbf{x}}^{\ast }\Vert$ est assez petit, $J_{\mathbf{F}}(\mathbf{x})$ est inversible et $\Vert J_{\mathbf{F}}(\mathbf{x}{)}^{-1}\Vert \le 2\Vert J_{\mathbf{F}}({\mathbf{x}}^{\ast }{)}^{-1}\Vert$. C'est l'analogue de la constante de la démonstration du théorème 2.5: loin d'une jacobienne singulière, l'inverse reste borné.

Le reste suit la démonstration scalaire, avec la formule de Taylor intégrale au lieu du reste de Lagrange, puisque celui-ci n'existe pas pour une fonction vectorielle. Posons ${\mathbf{e}}^{(k)}={\mathbf{x}}^{(k)}-{\mathbf{x}}^{\ast }$ et $J_k=J_{\mathbf{F}}({\mathbf{x}}^{(k)})$. Comme ,

En soustrayant ${\mathbf{x}}^{\ast }$ de (2.29) et en multipliant par $J_k$,

Les deux points où la jacobienne est évaluée sont à distance $(1-t)\Vert {\mathbf{e}}^{(k)}\Vert$ l'un de l'autre, et ${\int }_0^1(1-t)\mathrm{d}t=\frac{1}{2}$; le caractère lipschitzien donne donc

ce qui est (2.30) dès que ${\mathbf{x}}^{(k)}$ est dans la boule du lemme. Posons $K=\gamma \Vert J_{\mathbf{F}}({\mathbf{x}}^{\ast }{)}^{-1}\Vert$ et choisissons $\delta$ inférieur au rayon de cette boule et à . Si , alors : la suite reste dans la boule et converge, et (2.30) donne l'ordre 2 — exactement l'argument du théorème 2.5.

La constante de (2.30) fait apparaître $\Vert J_{\mathbf{F}}({\mathbf{x}}^{\ast }{)}^{-1}\Vert$, et ce n'est pas un hasard: c'est l'analogue de $1/|f^{\prime }(x^{\ast })|$, le conditionnement du problème «trouver la racine» de la relation (2.23). Une jacobienne presque singulière au point cherché — l'équivalent vectoriel d'une racine presque double — ralentit Newton et dégrade la précision atteignable, pour les mêmes raisons qu'en dimension un.

La figure 2.4 montre ce que fait (2.29) géométriquement. Chaque équation $F_i(\mathbf{x})=0$ définit une courbe du plan; son approximation affine en ${\mathbf{x}}^{(k)}$ définit une droite, parallèle à la tangente de la courbe de niveau de $F_i$ qui passe par ${\mathbf{x}}^{(k)}$. Le pas de Newton va à l'intersection des deux droites. En dimension , ce sont hyperplans, et leur intersection est la solution du système linéaire. On comprend aussi d'un coup d'œil les échecs: si les deux droites sont presque parallèles — jacobienne presque singulière, au sens exact du chapitre 4 — leur intersection part très loin, exactement comme la tangente presque horizontale de l'échec 1.

Critères d'arrêt. Tout ce que la section précédente a dit vaut encore, avec deux nuances. D'abord, pour Newton près d'une racine simple, ${\mathbf{e}}^{(k)}=-{\mathbf{\delta }}^{(k)}+O(\Vert {\mathbf{e}}^{(k)}{\Vert }^2)$: le pas est un excellent estimateur de l'erreur de l'itéré précédent, et c'est le critère du code ci-dessus, relatif pour les grandes composantes et absolu près de zéro. Ensuite, le résidu dépend désormais de l' de chaque équation: multiplier la seconde équation de (2.27) par ne change pas la solution mais change complètement la norme du résidu, donc le moment où un critère sur se déclenche. Avant de comparer des résidus entre eux, on met les équations à l'échelle — on les rend sans dimension, ou du même ordre de grandeur au point de départ. Et l'on garde le plafond , le drapeau d'échec, et l'amortissement du Newton amorti, qui en dimension n'est plus une option: aucun encadrement ne vient au secours d'une itération partie trop loin.

Quand la jacobienne coûte trop cher

Les $n^2$ dérivées partielles sont souvent le vrai obstacle: $\mathbf{F}$ peut sortir d'un code de simulation que personne ne sait dériver à la main. Trois remèdes, du plus simple au plus fin.

La jacobienne par différences finies. On approche la colonne $j$ par une différence avant, au sens du chapitre 8:

Le choix $h\approx \sqrt{{\epsilon }_{\text{mach}}}$ est celui qui équilibre l'erreur de troncature de la différence avant et l'erreur d'arrondi, comme le chapitre 8 le démontre. Coût: $n$ évaluations supplémentaires de $\mathbf{F}$ par itération. La jacobienne n'est plus exacte qu'à $10^{-8}$ près environ, donc la convergence n'est plus rigoureusement quadratique; en pratique on ne voit pas la différence: sur l'exemple 2.8, avec un pas de , il faut toujours cinq itérations.

La jacobienne gelée. On calcule et l'on factorise $J_{\mathbf{F}}({\mathbf{x}}^{(0)})$ une fois, puis on la réutilise à chaque pas: c'est la méthode de la corde, déjà rencontrée à la question 2.8 en dimension un. Chaque itération ne coûte plus qu'une évaluation de $\mathbf{F}$ et une descente-remontée en $2n^2$ opérations, mais la convergence redevient linéaire. Sur l'exemple, depuis $(2;0)$, il faut 22 itérations pour un résidu de $10^{-12}$, contre 5. Le compromis usuel est intermédiaire: on recalcule la jacobienne toutes les quelques itérations, ou seulement quand la convergence ralentit.

Les méthodes de quasi-Newton. La méthode de la sécante remplaçait $f^{\prime }(x_k)$ par une pente calculée sur les deux derniers itérés. Sa généralisation, due à Broyden, met à jour une approximation $B_k$ de la jacobienne de sorte qu'elle reproduise le dernier accroissement observé — la condition de la sécante $B_{k+1}{\mathbf{s}}_k={\mathbf{y}}_k$, avec et —, en la modifiant le moins possible:

La correction est de rang un, ce qui permet de mettre à jour la factorisation en $O(n^2)$ au lieu de la refaire, et chaque itération ne demande qu'une évaluation de $\mathbf{F}$. La convergence est superlinéaire (théorème de Dennis et Moré, que nous admettons: sa démonstration est longue et n'apporte pas d'idée nouvelle au regard du théorème 2.7). Sur l'exemple, partie de $B_0=J_{\mathbf{F}}({\mathbf{x}}^{(0)})$, la méthode de Broyden atteint un résidu de en 7 itérations: deux de plus que Newton, sans aucune dérivée après le départ.

Synthèse

• Localiser d'abord. Le théorème des valeurs intermédiaires donne l'existence d'une racine dès que $f$ change de signe, jamais l'unicité — qui vient de la monotonie — et il ne voit pas les racines de multiplicité paire, qui ne traversent pas l'axe.
• La bissection est la seule méthode à garantie inconditionnelle: $|c_k-x^{\ast }|\le (b-a)/2^{k+1}$, connu de calculer, soit 20 itérations pour sur un intervalle de longueur 1, et environ 3,32 itérations par chiffre décimal. Son erreur effective peut augmenter d'un pas au suivant; c'est la longueur de l'encadrement qui décroît, pas la distance à la racine.

Exercices

Vous pouvez afficher le corrigé directement sous chaque énoncé après avoir cherché la solution.

Références

• Quarteroni, A., Sacco, R. et Saleri, F., Méthodes numériques — algorithmes, analyse et applications, Springer, Milan, chap. 6 (résolution des équations non linéaires, méthodes de point fixe et de Newton).
• Rappaz, J. et Picasso, M., Introduction à l'analyse numérique, Presses polytechniques et universitaires romandes, Lausanne, chap. 2.
• Burden, R. L. et Faires, J. D., Numerical Analysis, Cengage, Boston, chap. 2 (Solutions of Equations in One Variable), avec le traitement des racines multiples et de la méthode de Steffensen.
• Quarteroni, A., Saleri, F. et Gervasio, P., Calcul scientifique, Springer, Milan, chap. 2.
• Brent, R. P., Algorithms for Minimization without Derivatives, Prentice-Hall, Englewood Cliffs, chap. 4 (l'hybride bissection–sécante–interpolation quadratique utilisé par les bibliothèques actuelles).
• Schatzman, M., Analyse numérique — une approche mathématique, Dunod, Paris, chap. 4.
• Dennis, J. E. et Schnabel, R. B., Numerical Methods for Unconstrained Optimization and Nonlinear Equations, SIAM, Philadelphie, chap. 5 à 8 (Newton pour les systèmes, jacobienne par différences finies, méthode de Broyden).