Résolution numérique d’équations non linéaires

Cours ENPC — Pratique du calcul scientifique

Introduction

Objectif

But de ce chapitre

Étudier les méthodes numériques pour trouver les solutions de \[ \mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{0}}, \] où \(\mathbf{\boldsymbol{f}}\colon \mathbf R^n \to \mathbf R^n\) (nombre d’inconnues = nombre d’équations).

Applications

  • Calculer \(\sqrt{2}\) en résolvant \(x^2 = 2\) ;

  • Résoudre des problèmes aux limites à deux points, par exemple

    \[ u''(t) = f\bigl(t, u(t)\bigr), \qquad u(0) = u(1) = 0. \]

  • Résoudre des équations aux dérivées partielles non linéaires, par exemple \[ \nabla \cdot \bigl(\alpha(u)\nabla u \bigr) = f. \]

Remarque

Un système linéaire de la forme \(\mathsf A \mathbf{\boldsymbol{x}} = \mathbf{\boldsymbol{b}}\) peut se mettre sous cette forme avec \[ \mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}) = \mathsf A \mathbf{\boldsymbol{x}} - \mathbf{\boldsymbol{b}}. \] Avec cette notation, la méthode de Richardson pour les systèmes linéaires s’écrit \[ \mathbf{\boldsymbol{x}}_{k+1} = \mathbf{\boldsymbol{x}}_k - \omega \mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}_k). \] Questions :

  • Une telle méthode converge-t-elle pour des équations non linéaires générales ?

  • Quelle est la vitesse de convergence ?

  • Quelle est la valeur optimale de \(\omega\) ?

Approche la plus simple : la méthode de dichotomie

Cadre et idée

Supposons que \(f\colon [a, b] \to \mathbf R\) est continue et \(f(a) \leqslant 0 \leqslant f(b).\)

\({\color{green} \rightarrow}\) Alors il existe une racine de \(f\) dans \([a, b]\). Soit \(x = \frac{a + b}{2}\).

  • Si \(f(a)f(x) \leqslant 0\), alors il existe une racine de \(f\) dans \([a, x]\) ;

  • Si \(f(a)f(x) > 0\), alors il existe une racine de \(f\) dans \([x, b]\).

function bisection(f, a, b, ε = 1e-10)
    @assert f(a) * f(b) ≤ 0
    while abs(b - a) ≥ ε
        x = (a + b) / 2
        a, b = f(a) * f(x) ≤ 0 ? [a, x] : [x, b]
    end
    return (a + b) / 2
end

println(cbrt(2))
bisection(x -> x^3 - 2, 0, 2)
1.2599210498948732
1.259921049902914

Convergence

On a \[ |b_i - a_i| = \frac{|b - a|}{2^i} \] Les deux suites \((a_i)\) et \((b_i)\) convergent vers une racine \(x_*\) de \(f\).

\(~\)

Commentaires

  • Fonctionne uniquement dans le cas unidimensionnel ;

  • Robuste mais lent…

Approche la plus simple : la méthode de dichotomie

Quantification de la vitesse de convergence

Définition

On dit qu’une suite \((\mathbf{\boldsymbol{x}}_k)_{k\in \mathbf N}\) converge vers \(\mathbf{\boldsymbol{x}}_*\)

  • linéairement si \[ \lim_{k\to\infty} \frac{{\lVert {\mathbf{\boldsymbol{x}}_{k+1} - \mathbf{\boldsymbol{x}}_*} \rVert}}{{\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert}} \in (0, 1). \]

  • superlinéairement si \[ \lim_{k\to\infty} \frac{{\lVert {\mathbf{\boldsymbol{x}}_{k+1} - \mathbf{\boldsymbol{x}}_*} \rVert}}{{\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert}} = 0. \]

  • avec ordre \(q\) si \(\mathbf{\boldsymbol{x}}_k \to \mathbf{\boldsymbol{x}}_*\) et \[ \lim_{k\to\infty} \frac{{\lVert {\mathbf{\boldsymbol{x}}_{k+1} - \mathbf{\boldsymbol{x}}_*} \rVert}}{{\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert}^q} \in (0, \infty). \] Si \(q = 2\), on dit que la suite converge quadratiquement.

Exemples

  • La plupart des méthodes itératives pour les équations linéaires convergent linéairement.

  • La méthode de dichotomie converge linéairement.

Méthodes à point fixe

Méthodes à point fixe

Idée

Réécrire

\[\mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{0}} \qquad \Leftrightarrow \qquad \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{x}}\]

et utiliser le schéma itératif

\[ \mathbf{\boldsymbol{x}}_{k+1} = \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_k) \]

Exemples (plus de détails dans la section suivante)

  • Méthode de la corde : pour \(\alpha \neq 0\) \[ \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{x}} - \frac{\mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}})}{\alpha} \] \({\color{green} \rightarrow}\) Dans le cas linéaire, c’est l’itération de Richardson avec \(\omega = \alpha^{-1}\).

  • Méthode de Newton-Raphson : \[ \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{x}} - \mathsf J_f(\mathbf{\boldsymbol{x}})^{-1} \mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}), \] où \(\mathsf J_f(\mathbf{\boldsymbol{x}})\) est la matrice jacobienne de \(\mathbf{\boldsymbol{f}}\) en \(\mathbf{\boldsymbol{x}}\).

Proposition : convergence sous hypothèse de Lipschitz

Si \(\mathbf{\boldsymbol{F}}\) est lipschitzienne avec \(L < 1\), \[ \forall (\mathbf{\boldsymbol{x}}, \mathbf{\boldsymbol{y}}) \in \mathbf R^n \times \mathbf R^n, \qquad {\lVert {\mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) - \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{y}})} \rVert} \leqslant L {\lVert {\mathbf{\boldsymbol{x}} - \mathbf{\boldsymbol{y}}} \rVert} \] alors

  • Il existe un unique point fixe \(\mathbf{\boldsymbol{x}}_*\) de \(\mathbf{\boldsymbol{F}}\).

  • La suite \((\mathbf{\boldsymbol{x}}_k)_{k\in \mathbf N}\) converge exponentiellement vers \(\mathbf{\boldsymbol{x}}_*\) : \[ \forall k \in \mathbf N, \qquad {\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert} \leqslant L^k {\lVert {\mathbf{\boldsymbol{x}}_0 - \mathbf{\boldsymbol{x}}_*} \rVert}. \]

Esquisse de preuve

  • Montrer que \((\mathbf{\boldsymbol{x}}_k)\) est de Cauchy, donc convergente vers un certain \(\mathbf{\boldsymbol{x}}_* \in \mathbf R^n\).

  • Montrer que \(\mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_*) = \mathbf{\boldsymbol{x}}_*\).

  • Déduire la convergence exponentielle de \(\mathbf{\boldsymbol{x}}_{k+1} - \mathbf{\boldsymbol{x}}_* = \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_k) - \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_*)\).

Affaiblissement de l’hypothèse de Lipschitz globale

Proposition : convergence locale sous condition de Lipschitz locale

Supposons que \(\mathbf{\boldsymbol{x}}_*\) est un point fixe de \(\mathbf{\boldsymbol{F}}\) et \[ \forall \mathbf{\boldsymbol{x}} \in B_{\delta}(\mathbf{\boldsymbol{x}}_*), \qquad {\lVert {\mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) - \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_*)} \rVert} \leqslant L {\lVert {\mathbf{\boldsymbol{x}} - \mathbf{\boldsymbol{x}}_*} \rVert} \tag{1}\] pour \(L < 1\). Si \(\mathbf{\boldsymbol{x}}_0 \in B_{\delta}(\mathbf{\boldsymbol{x}}_*)\), alors \(\mathbf{\boldsymbol{x}}_k \in B_{\delta}(\mathbf{\boldsymbol{x}}_*)\) pour tout \(k \in \mathbf N\) et \[ {\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert} \leqslant L^k {\lVert {\mathbf{\boldsymbol{x}}_0 - \mathbf{\boldsymbol{x}}_*} \rVert}. \]

\(~\)

Proposition : critères de convergence locale basés sur la jacobienne

Supposons que \(\mathbf{\boldsymbol{F}}\) est différentiable en un point fixe \(\mathbf{\boldsymbol{x}}_*\).

  • Il existe \(\delta > 0\) tel que Équation 1 est satisfait ssi \({\lVert {\mathsf J_F(\mathbf{\boldsymbol{x}}_*)} \rVert} < 1\).

  • Il existe \(\delta > 0\) tel que Équation 1 est satisfait dans une certaine norme induite ssi \(\rho\bigl(\mathsf J_F(\mathbf{\boldsymbol{x}}_*)\bigr) < 1\).

Méthodes à point fixe : convergence superlinéaire

Proposition : convergence superlinéaire

Supposons que \(\mathbf{\boldsymbol{x}}_*\) est un point fixe de \(\mathbf{\boldsymbol{F}}\) et que \(\mathsf J_F(\mathbf{\boldsymbol{x}}_*) = 0\). Si \(\mathbf{\boldsymbol{x}}_k \to \mathbf{\boldsymbol{x}}_*\) lorsque \(k \to \infty\), alors \[ \lim_{k \to \infty} \frac{{\lVert {\mathbf{\boldsymbol{x}}_{k+1} - \mathbf{\boldsymbol{x}}_*} \rVert}}{{\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert}} = 0. \]

Remarque

Sous cette hypothèse, il existe \(\delta > 0\) tel que \((\mathbf{\boldsymbol{x}}_k)_{k\geqslant 0}\) est une suite convergeant vers \(\mathbf{\boldsymbol{x}}_*\) pour tout \(\mathbf{\boldsymbol{x}}_0 \in B_{\delta}(\mathbf{\boldsymbol{x}}_*)\).

Preuve

On a \[ \frac{{\lVert {\mathbf{\boldsymbol{x}}_{k+1} - \mathbf{\boldsymbol{x}}_*} \rVert}}{{\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert}} = \frac{{\lVert {\mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_{k}) - \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_*)} \rVert}}{{\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert}} = \frac{{\lVert {\mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_{k}) - \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}_*) - \mathsf J_F(\mathbf{\boldsymbol{x}}_*) (\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*)} \rVert}}{{\lVert {\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_*} \rVert}}. \] Puisque \(\mathbf{\boldsymbol{x}}_k - \mathbf{\boldsymbol{x}}_* \to \mathbf{\boldsymbol{0}}\) lorsque \(k \to \infty\), le membre de droite converge vers 0 par définition de la matrice jacobienne.

Exemples de méthodes à point fixe

La méthode de la corde

Description de la méthode

La méthode de la corde est l’itération \[ \mathbf{\boldsymbol{x}}_{k+1} = \mathbf{\boldsymbol{x}}_k - \mathsf A^{-1} \mathbf{\boldsymbol{f}}(x_k), \] avec \(\mathsf A \in \mathbf R^{n \times n}\) un paramètre matriciel inversible.

\(~\)

Remarques

  • C’est une méthode à point fixe avec \(\mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{x}} - \mathsf A^{-1} \mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}})\).

  • Clairement \(\mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{x}}\) si et seulement si \(\mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{0}}\).

  • Lorsque \(n = 1\), l’itération s’écrit \[ x_{k+1} = x_k - \frac{f(x_k)}{\alpha}, \] \({\color{green} \rightarrow}\) Dans le cas linéaire, c’est l’itération de Richardson.

Étude de convergence lorsque \(n=1\)

Dans ce cas \[ F(x) = x - \frac{f(x)}{\alpha}, \qquad F'(x) = 1 - \frac{f'(x)}{\alpha} \] On déduit :

  • L’itération converge localement près du point fixe si \[ \left| 1 - \frac{f'(x_*)}{\alpha} \right| < 1. \] \({\color{green} \rightarrow}\alpha\) doit être du même signe que \(f'(x_*)\) et suffisamment grand : \[ \left|\alpha\right|>\frac{\left|f'(x_*)\right|}{2} \]

  • La convergence est superlinéaire si \[ \alpha = f'(x_*) \]

Méthodes à point fixe : interprétation géométrique avec \(F\)

\[ F(x) = x - \frac{1}{\alpha}\tanh(x-1), \qquad x_* = 1, \]

\[ F'(x_*) = 1 - \frac{1}{\alpha\,\cosh^2(x-1)}=1-\frac{1}{\alpha} \]

Méthodes à point fixe : interprétation géométrique avec \(F\)

\[ F(x) = x - \frac{1}{\alpha}\tanh(x-1), \qquad x_* = 1, \]

\[ F'(x_*) = 1 - \frac{1}{\alpha\,\cosh^2(x-1)}=1-\frac{1}{\alpha} \]

Méthodes à point fixe : interprétation géométrique avec \(F\)

\[ F(x) = x - \frac{1}{\alpha}\tanh(x-1), \qquad x_* = 1, \]

\[ F'(x_*) = 1 - \frac{1}{\alpha\,\cosh^2(x-1)}=1-\frac{1}{\alpha} \]

Méthodes à point fixe : interprétation géométrique avec \(F\)

\[ F(x) = x - \frac{1}{\alpha}\tanh(x-1), \qquad x_* = 1, \]

\[ F'(x_*) = 1 - \frac{1}{\alpha\,\cosh^2(x-1)}=1-\frac{1}{\alpha} \]

La méthode de la corde : interprétation géométrique avec \(f\)

Interprétation géométrique

À chaque itération :

  • Approximer \(f(x)\) par la droite \[ \widehat f_k(x) = f(x_k) + \alpha (x - x_k) \]

  • Soit \(x_{k+1}\) la racine de \(\widehat f_k\).

La méthode de la corde : interprétation géométrique

La méthode de la corde : expérience numérique

function chord(f, α, x; ε = 1e-12, maxiter = 10^7)
    for i in 1:maxiter
        x = x - f(x) / α
        if abs(f(x)) < ε
            return x, i
        end
    end
    error("Failed to converge!")
end

Problème jouet

Appliquons la méthode de la corde à \[ f(x) = x^2 - 2. \] Dans ce cas \[ F(x) = x - \frac{x^2 - 2}{\alpha}, \qquad F'(x) = 1 - \frac{2x}{\alpha}. \] Ici \(x_* = \sqrt{2}\) et \(f'(x_*) = 2\sqrt{2}\).

chord(x -> x^2 - 2, 2√2, 1)
(1.4142135623729684, 4)

\(~\)

chord(x -> x^2 - 2, √2 - .01, 1)
"Failed to converge!"

\(~\)

chord(x -> x^2 - 2, √2 + .01, 1)
(1.414213562373444, 1893)

\(~\)

chord(x -> x^2 - 2, 10, 1)
(1.4142135623728198, 85)

\(~\)

chord(x -> x^2 - 2, 100, 1)
(1.4142135623727494, 975)

La méthode de Newton-Raphson

Motivation

  • La méthode de la corde est optimale avec \(\mathsf A = \mathsf J_f(\mathbf{\boldsymbol{x}}_*)\) ;

  • Cependant, \(\mathbf{\boldsymbol{x}}_*\) est inconnu a priori…

  • Méthode de Newton-Raphson : utiliser \[ \mathsf A = \mathsf A(k) = \mathsf J_f(\mathbf{\boldsymbol{x}}_k). \] ce qui conduit à l’itération \[ \mathbf{\boldsymbol{x}}_{k+1} = \mathbf{\boldsymbol{x}}_k - \mathsf J_f(\mathbf{\boldsymbol{x}}_k)^{-1} \mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}_k) \]

Version unidimensionnelle

En dimension un, l’itération de Newton-Raphson s’écrit \[ x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)} \]

Remarques

  • La méthode nécessite de résoudre un système linéaire à chaque étape.

  • L’itération est bien définie seulement si \(\mathbf{\boldsymbol{J}}_f(\mathbf{\boldsymbol{x}}_k)\) est non singulière.

  • La méthode est une itération à point fixe avec \[ \mathbf{\boldsymbol{F}}(\mathbf{\boldsymbol{x}}) = \mathbf{\boldsymbol{x}} - \mathsf J_f(\mathbf{\boldsymbol{x}})^{-1} \mathbf{\boldsymbol{f}}(\mathbf{\boldsymbol{x}}). \]

  • Dimension 1 : Si \(f'(x_*) \neq 0\), alors \[ F'(x_*) = 0. \] \({\color{green} \rightarrow}\) convergence superlinéaire locale.

  • Dimension \(n\) : Si \(\mathsf J_f(\mathbf{\boldsymbol{x}}_*)\) est non singulière, alors \[ \mathsf J_F(\mathbf{\boldsymbol{x}}_*) = \mathsf 0. \] \({\color{green} \rightarrow}\) convergence superlinéaire locale.

La méthode de Newton-Raphson : analyse de convergence

Théorème : Convergence quadratique

Supposons que \(f \in C^2(\mathbf R)\), que \(f'(x_*) \neq 0\), et que l’itération de Newton-Raphson converge vers \(x_*\). Alors \[ \lim_{k \to \infty} \frac{{\lvert {x_{k+1} - x_*} \rvert}}{{\lvert {x_{k} - x_*} \rvert}^2} = \left| \frac{f''(x_*)}{2 f'(x_*)} \right|. \]

Preuve

On remarque que \[ x_{k+1} - x_* = x_{k} - \frac{f(x_k)}{f'(x_k)} - x_* = \frac{1}{f'(x_k)} \underbrace{\Bigl(f'(x_k)(x_{k} - x_*) - f(x_k)\Bigr)}_{= \frac{1}{2} f''(\xi_k) (x_{k} - x_*)^2}. \] pour un certain \(\xi_k\) entre \(x_k\) et \(x_*\). On conclut que \[ x_{k+1} - x_* = \frac{f''(\xi_k) (x_k - x_*)^2}{2f'(x_k)}. \]

La méthode de Newton-Raphson : expérience numérique

function newton_raphson(f, df, x, maxiter; ε = 1e-12)
    for i in 1:maxiter
        x -= f(x) / df(x)
        if abs(f(x)) < ε
            return x, i
        end
    end
    return "Failed to converge!"
end;

Problème jouet

Appliquons la méthode de Newton-Raphson à \[ f(x) = (x^2 - 2)(x - 3)^4. \]

  • Racines simples \(\sqrt{2}\) et \(-\sqrt{2}\).

  • Racine quadruple \(3\).

    \({\color{green} \rightarrow}\) la convergence superlinéaire ne s’applique pas.

println(√2)
1.4142135623730951

\(~\)

# Importation du package pour la différentiation automatique
using Zygote

g(x) = (x^2 - 2) * (x - 3)^4
dg = g'
maxiter = 100
newton_raphson(g, dg, 1, maxiter)
(1.4142135623730951, 6)

\(~\)

newton_raphson(g, dg, 100, maxiter)
(3.000478179164197, 50)

\(~\)

newton_raphson(g, dg, -100, maxiter)
(-1.414213562373095, 25)

La méthode de Newton-Raphson : illustration géométrique

Interprétation géométrique

À chaque itération :

  • Approximer \(f(x)\) par sa tangente en \(x_k\) \[ \widehat f_k(x) = f(x_k) + f'(x_k) (x - x_k) \]

  • Soit \(x_{k+1}\) la racine de \(\widehat f_k\).

La méthode de Newton-Raphson : illustration géométrique