FEUILLE DE TP 6

Méthode de Newton¶

La méthode de Newton nous permet de calculer numériquement les zéros d'une fonction $f:[a,b]\rightarrow \mathbb{R}$, c'est-à-dire les solutions de l'équation $f(x) = 0$.

Supposons $f$ de classe $C^1$ et $x_0 \in (a,b)$ un point arbitraire.

L'équation de la tangente au graphe de la fonction $f$ au point d'absisse $x_0$ est donnée par

\begin{equation} y = f(x_0) + f'(x_0)(x-x_0)\,. \end{equation}

Pour obtenir une approximation de la solution, on pose $y = 0$ et si $f'(x_0) \neq 0$ on obtient

\begin{equation} x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}\,. \end{equation}

On peut alors répéter le même raisonnement en considérant $x_1$ au lieu de $x_0$ pour calculer l'approximation de $x_2$ et ainsi de suite.

La méthode de Newton consiste donc à construire une suite d'approximations $x_k$ de la solution $x^*$ de l'équation $f(x) = 0$, définie par

\begin{equation} x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}\,, \qquad k\in \mathbb{N}\, . \qquad (1) \end{equation}
In [1]:
%matplotlib inline
%config InlineBackend.figure_format = 'retina'

import numpy as np
import matplotlib.pyplot as plt

color_bleu = 'xkcd:dusky blue'
color_vert = 'xkcd:grassy green'
color_orange = 'xkcd:orangish'
color_rouge = 'xkcd:blood orange'
color_bleuclair = 'xkcd:sky blue'
color_vertclair = 'xkcd:grass green'
color_marron = 'xkcd:burnt sienna'
color_violet = 'xkcd:blue violet'

liste_color = [
    color_bleu,
    color_vert,
    color_orange,
    color_rouge,
    color_bleuclair,
    color_vertclair,
    color_marron,
    color_violet
]

Calcul de valeurs approchées de nombres irrationnels¶

Nous allons utiliser la méthode de Newton pour calculer des approximations de quelques nombres irrationnels. On considère

\begin{equation} f_0(x) = x^3-2~ \text{ avec } ~ x\in [1,4] \, , \quad f_1(x) = \sin^2(x)~ \text{ avec } ~ x\in [2,4] \end{equation}

En appliquant la méthode de Newton à $f_0$ et $f_1$, on espère obtenir respectivement une approximation de $2^{1/3}$ et $\pi$.

Question 1

Définissez les fonctions $f_0$ et $f_1$, ainsi que leurs dérivés $f'_0$ et $f'_1$, et tracez-les sur leur intervalle de définition.

In [2]:
liste_f = [
    {
        'name': 'f_0',
        'fonction': lambda x: x**3-2,
        'derive': lambda x: 3*x**2,
        'a': 1,
        'b': 4,
        'zero': 2**(1/3)
    },
    {
        'name': 'f_1',
        'fonction': lambda x: np.sin(x)**2,
        'derive': lambda x: 2*np.cos(x)*np.sin(x),
        'a': 2,
        'b': 4,
        'zero': np.pi
    }
]
In [3]:
taille = 6
fig = plt.figure(figsize=(taille*len(liste_f), taille))

for k, dico in enumerate(liste_f):
    ax = fig.add_subplot(1, len(liste_f), k+1)
    ax.grid(False)
    a, b = dico['a'], dico['b']
    f, df = dico['fonction'], dico['derive']
    x = np.linspace(a, b, 128)
    ax.plot(x, f(x), linewidth=2, color=color_bleu, label="fonction")
    ax.plot(x, df(x), linewidth=2, color=color_orange, label="dérivée")
    ax.scatter([dico['zero']], [0], color=color_rouge, s=50)
    ax.axhline(0, color=color_vert, linestyle='dashed')
    ax.legend()
    ax.set_title(f"Fonction ${dico['name']}$", fontsize=16)

Question 2

Avant d'implémenter la méthode de Newton (1) on a besoin de choisir un critère d’arrêt: à chaque itération $k$ de la méthode ce critère nous permet de décider s'il faut caluler l'approximation suivante ou arrêter l'algorithme.

On choisira d’implémenter ici un critère d’arrêt donné par:

\begin{equation} |f(x_k)|< \epsilon \quad \text{ OU } \quad k> N_{max} \end{equation}

où $\epsilon>0$ est un réel fixé et $N_{max} \in \mathbb{N}$ est le nombre maximal d'itérations.

Proposez une fonction newton prenant en arguments

  • une fonction $f$
  • sa dérivée $f'$
  • une condition initiale $x_0$
  • une tolérance $\epsilon>0$ pour le critère d’arrêt
  • un nombre d'itération maximal $N_{max}$ pour le critère d’arrêt

et qui donne en sortie

  • le nombre d'itérations nécessaires pour vérifier le critère d'arrêt,
  • la solution approchée $x^*$ calculée avec la méthode de Newton,
  • xL un ndarray contenant tous les termes de la suite $(x_n)$ construite par la méthode de Newton.
In [4]:
def newton(f, df, x0, epsilon, Nmax):
    """
    Algorithme de Newton pour résoudre f(x) = 0 
    
    Parameters
    ----------
    f: function
        fonction f
    df: function
        dérivé de f
    x0: float 
        approximation initiale  
    epsilon: float
        tolerance dans le critere d'arret
    Nmax: int
        nombre d'iterations maximal

    Returns
    -------
    niter: float
        nombre d'iteration
    xstar: float
        solution approchée
    xL: ndarray
        la liste des itérés de la méthode
    """
    x = x0
    fx, dfx = f(x), df(x)
    niter = 0
    xL = [x0]
    while abs(fx) >= epsilon and niter < Nmax:
        x -= fx / dfx
        fx, dfx = f(x), df(x)
        niter += 1
        xL.append(x)
    return niter, x, np.asanyarray(xL)

Question 3

  1. Reprenez les graphiques de la question 1.
  2. Calculez la solution approchée $x^\star$ donnée par la méthode de Newton pour chacun des cas et représentez la à l'aide d'un scatter sur la figure. Pour l'initialisation, vous prendrez $x_0=4$.
  3. Sur une autre figure, tracez pour chaque fonction le nuages de points $(|e_n|, |e_{n+1}|)$, où $e_n=x_n-\ell$ (Pour $f_0$, $\ell=2^{1/3}$ et pour $f_1$, $\ell=\pi$). Vous pouvez utiliser une échelle logarithmique ! Les courbes sont-elles cohérentes avec la théorie ? Quelle est l'ordre de la méthode ?
In [5]:
x0 = 4
epsilon, Nmax = 1.e-14, 100
taille = 6
fig = plt.figure(figsize=(taille*len(liste_f), 2*taille))

for k, dico in enumerate(liste_f):
    ax = fig.add_subplot(2, len(liste_f), k+1)
    ax.grid(False)
    a, b = dico['a'], dico['b']
    f, df = dico['fonction'], dico['derive']
    x = np.linspace(a, b, 128)
    ax.plot(x, f(x), linewidth=2, color=color_bleu, label="fonction")
    ax.plot(x, df(x), linewidth=2, color=color_orange, label="dérivée")
    ax.axhline(0, color=color_vert, linestyle='dashed')
    ax.legend()
    ax.set_title(f"Fonction ${dico['name']}$", fontsize=16)
    niter, xstar, xL = newton(f, df, x0, epsilon, Nmax)
    ax.scatter([xstar], [0], color=color_rouge, s=50)
    
    eL = abs(xL-dico['zero'])
    ax_e = fig.add_subplot(2, len(liste_f), len(liste_f)+k+1)
    ax_e.set_xscale('log', base=2)
    ax_e.set_yscale('log', base=2)
    ax_e.scatter(
        eL[:-1], eL[1:],
        color=color_orange, label='Newton'
    )
    eL = eL[eL>0]
    p = np.polyfit(np.log(eL[:-1]), np.log(eL[1:]), 1)[0]
    ax_e.plot(
        eL[:-1], eL[1]*(eL[:-1]/eL[0])**p,
        color='black', alpha=0.5, label=f'pente ${p:5.3f}$'
    )
    ax_e.legend()
    ax_e.set_title(f"${niter}$ itérations", fontsize=16, color='black')
    ax_e.set_xlabel(r"$|e_{n}|$", color='black')
    ax_e.set_ylabel(r"$|e_{n+1}|$", color='black')

Un exemple avec plusieurs solutions¶

Si on applique la méthode de Newton à une fonction $f$ telle que $f(x) = 0$ a plusieurs solutions, on peut s'attendre à deux résultats possibles: soit la méthode ne converge pas, soit elle converge vers un seul des zéros de $f$. En particulier la méthode peut converger vers des zéros différents si on choisit de points initiaux $x_0$ différents.

L'ensemble des $x_0$ tels que la méthode converge vers la solution $x^*$ s'appelle le bassin d'attraction de $x^*$.

Question 4

Soit $f$ la fonction telle que $f(x) =x^3 - 4 x+ 1$.

  • Combien l’équation $f(x) = 0$ admet-elle de solutions ? (vérification graphique)
  • Recherchez ces solutions avec la fonction newton en initialisant successivement la méthode avec différentes valeurs initiales $x_0$.
  • Tracez la fonction $f(x)$ et vérifiez le résultat.
  • Ajoutez sur la figure la ligne brisée $(x_0, 0)$, $(x_0, f(x_0))$, $(x_1, 0)$, $(x_1, f(x_1))$, $\ldots$.
In [6]:
def f(x):
    return x**3 - 4*x + 1

def df(x): 
    return 3*x**2 - 4

fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(1, 1, 1)
ax.grid(False)
x = np.linspace(-3, 3, 128)
ax.plot(x, f(x), linewidth=2, color=color_bleu, label="fonction")
ax.plot(x, df(x), linewidth=2, color=color_orange, label="dérivée")
ax.axhline(0, color=color_vert, linestyle='dotted')
ax.set_title(r'Fonction $x^3-4x+1$', fontsize=16);
In [7]:
fig = plt.figure(figsize=(6, 6))
ax = fig.add_subplot(1, 1, 1)
ax.grid(False)
x = np.linspace(-3, 3, 128)
ax.plot(x, f(x), linewidth=2, color=color_bleu, label="fonction")
ax.plot(x, df(x), linewidth=2, color=color_orange, label="dérivée")
ax.axhline(0, color=color_vert, linestyle='dotted')
ax.set_title(r'Fonction $x^3-4x+1$')

print("-"*28)
print("|  x0  |       x*      | n |")
print("-"*28)
k = 0
for x0 in [-2.5, 0, 2.5]:
    k += 1
    n, x, xL = newton(f, df, x0, 1.e-14, 20)
    print(f"| {x0:4.1f} | {x:13.10f} | {n:1d} |")
    ax.scatter([x], [0], label=f"$x_0={x0}$", color=liste_color[2*k])
    abscisse = xL.repeat(2)
    ordonnee = f(abscisse)
    ordonnee[::2] = 0
    ax.plot(abscisse, ordonnee, linewidth=1, color=liste_color[2*k])
print("-"*28)
ax.legend();

----------------------------
|  x0  |       x*      | n |
----------------------------
| -2.5 | -2.1149075415 | 5 |
|  0.0 |  0.2541016884 | 4 |
|  2.5 |  1.8608058531 | 6 |
----------------------------

Méthode de Newton modifiée : la fausse position¶

Lorsque le calcul de la dérivée est impossible ou trop long, il est possible de remplacer le terme $f'(x_n)$ dans la méthode de Newton par le taux d'accroissement $$ \frac{f(x_n)-f(x_{n-1})}{x_n-x_{n-1}}. $$

La méthode revient alors, en partant de deux données initiales $x_0$ et $x_1$, à construire la suite $(x_n)$ définie par $$ x_{n+1} = x_n - \frac{f(x_n)(x_n-x_{n-1})}{f(x_n)-f(x_{n-1})} = \frac{x_{n-1}f(x_n)-x_nf(x_{n-1})}{f(x_n)-f(x_{n-1})}. $$

Cette méthode s'appelle la fausse position. Elle possède les mêmes propriétés que la méthode de Newton (suite qui n'est pas toujours constructible, convergence locale) mais est un peu moins rapide.

Reprenez les questions sur la méthode de Newton avec cette nouvelle méthode. Quel ordre de convergence observez-vous dans le cas d'un zéro simple et dans celui d'un zéro double ?

In [8]:
def fausse_position(f, x0, x1, epsilon, Nmax):
    """
    Algorithme de la fausse position pour résoudre f(x) = 0 
    
    Parameters
    ----------
    f: function
        fonction f
    df: function
        dérivé de f
    x0: float 
        approximation initiale  
    epsilon: float
        tolerance dans le critere d'arret
    Nmax: int
        nombre d'iterations maximal

    Returns
    -------
    niter: float
        nombre d'iteration
    xstar: float
        solution approchée
    """
    x, xold = x1, x0
    fx, fxold = f(x), f(xold)
    niter = 0
    xL = [x0, x1]
    while abs(fx) >= epsilon and niter < Nmax:
        xold, x = x, (x*fxold-xold*fx)/(fxold-fx)
        fxold, fx = fx, f(x)
        niter += 1
        xL.append(x)
    return niter, x, np.asanyarray(xL)
In [11]:
x0 = 4
x1 = 3.75
epsilon, Nmax = 1.e-14, 100
taille = 6
fig = plt.figure(figsize=(taille*len(liste_f), 2*taille))

for k, dico in enumerate(liste_f):
    ax = fig.add_subplot(2, len(liste_f), k+1)
    ax.grid(False)
    a, b = dico['a'], dico['b']
    f, df = dico['fonction'], dico['derive']
    x = np.linspace(a, b, 128)
    ax.plot(x, f(x), linewidth=2, color=color_bleu, label="fonction")
    ax.plot(x, df(x), linewidth=2, color=color_violet, label="dérivée")
    ax.axhline(0, color=color_vert, linestyle='dashed')
    ax.legend()
    ax.set_title(f"Fonction ${dico['name']}$", fontsize=16)

    ax_e = fig.add_subplot(2, len(liste_f), len(liste_f)+k+1)
    ax_e.set_xscale('log', base=2)
    ax_e.set_yscale('log', base=2)

    # newton
    niter, xstar, xL = newton(f, df, x0, epsilon, Nmax)
    ax.scatter([xstar], [0], color=color_orange, s=50)
    eL = abs(xL-dico['zero'])
    ax_e.scatter(
        eL[:-1], eL[1:],
        color=color_orange, label='Newton'
    )
    eL = eL[eL>0]
    p = np.polyfit(np.log(eL[:-1]), np.log(eL[1:]), 1)[0]
    ax_e.plot(
        eL[:-1], eL[1]*(eL[:-1]/eL[0])**p,
        color=color_orange, alpha=0.5, label=f'pente ${p:5.3f}$'
    )
    # fausse position
    niter, xstar, xL = fausse_position(f, x0, x1, epsilon, Nmax)
    ax.scatter([xstar], [0], color=color_vert, s=50)
    eL = abs(xL-dico['zero'])
    ax_e.scatter(
        eL[:-1], eL[1:],
        color=color_vert, label='fausse position'
    )
    eL = eL[eL>0]
    p = np.polyfit(np.log(eL[:-1]), np.log(eL[1:]), 1)[0]
    ax_e.plot(
        eL[:-1], eL[1]*(eL[:-1]/eL[0])**p,
        color=color_vert, alpha=0.5, label=f'pente ${p:5.3f}$'
    )
    
    ax_e.legend()
    ax_e.set_title("Convergence", fontsize=16, color='black')
    ax_e.set_xlabel(r"$|e_{n}|$", color='black')
    ax_e.set_ylabel(r"$|e_{n+1}|$", color='black')
In [ ]: