TP2 Méthode de Newton¶

Dans ce TP on va s'intéresser à trouver le minimum d'une fonction de deux manières différentes. Ce genre de problème trouve des applications en statistique, analyse de données et machine learning.

La fonction sur laquelle on va appliquer ces méthodes est

$$ f(x) = \textrm{erf}(x-2)+\textrm{erf}(-x-2),\qquad x\in \mathbb{R},$$

où la fonction $\textrm{erf}$ est définie par

$$ \textrm{erf}(x) = \frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2} dt.$$

Ainsi la dérivée de $f$ est donnée par

$$ f'(x) = \frac{2}{\sqrt{\pi}}\Big( e^{-(x-2)^2} - e^{-(x+2)^2} \Big) $$

Tout d'abord, importons les librairies dont nous aurons besoin et la fonction erf de la librairie scipy.special.

Avec la méthode de Newton on va chercher là où la dérivée $f'$ s'annule, c'est-à-dire on va résoudre l'équation $g(x) = 0$, où $g(x) = f'(x)$.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf
from numpy import pi

On peut donc appeller la plupart des fonctions mathématiques de la façon suivante

In [2]:

np.exp(0), np.cos(2*pi)

Out[2]:

(np.float64(1.0), np.float64(1.0))

Pour construire une grille de point sur laquelle on va évaluer une fonction, on peut utiliser np.linspace(début, fin, nbre_de_point) ou np.arange(début, fin, pas_entre_les_points). Par exemple pour l'intervalle $[0,1]$.

In [3]:

np.linspace(0, 1, 5)

Out[3]:

array([0.  , 0.25, 0.5 , 0.75, 1.  ])

In [4]:

np.arange(0, 1, 0.1)

Out[4]:

array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9])

Enfin, pour définir une fonction mathématique on peut utiliser la syntaxe suivante. Par exemple, pour une fonction qui donne les valeurs de la fonction $h(x) = x^2 + 1$, on peut utiliser

In [5]:

def h(x):
    return x**2 + 1

In [6]:

h(0)

Out[6]:

Pour faire une boucle for, on procède de la façon suivante:

In [7]:

for i in range(3):
    print(i)

0
1
2

In [8]:

for i in range(2,6):
    print(i)

1) Tracer $f$ et $f'$ sur l'intervalle $[-5,5]$.¶

In [9]:

def f(x):
    return erf(x-2)+erf(-x-2)

def df(x):
    return (2/np.sqrt(pi))*(np.exp(-(x-2)**2) - np.exp(-(x+2)**2))

In [10]:

x = np.linspace(-5, 5,100)

In [11]:

plt.plot(x,f(x), label = "fct")
plt.plot(x,df(x), label="dérivée")
plt.legend()

Out[11]:

<matplotlib.legend.Legend at 0x7ad49c269f10>

No description has been provided for this image

2) La méthode de Newton.¶

Un énoncé du théorème sur la méthode de Newton est le suivant:

Théorème:¶

Soient $g$ une fonction de classe $C^2$ sur un intervalle $I=[a-r, a+r]$, pour laquelle $g'(x)\neq 0$ sur $I$, et où $a\in I$ tel que $g(a)=0$. On pose

$$ M=\max_{x\in I} \Big|\frac{g''(x)}{g'(x)}\Big|\qquad\text{et}\qquad h = \min\Big(r,\frac{1}{M}\Big).$$

Alors, pour tout $x_0 \in [a-h, a+h]$, la suite de points définie par

$$ x_{n+1} = x_n -\frac{g(x_n)}{g'(x_n)},\qquad n=0,1,\dots $$

converge vers $a$.

La méthode de Newton consiste donc à construire la suite de points $(x_n)_n$ qui nous donnera une approximation du point $a$ où $g$ s'annule ($g(a) = 0$, et donc $f'(a)=0$ pour notre problème). On remarquera qu'on a la dérivé $g'$ qui intervient (et donc $f''$ pour notre problème). Il faut donc être capable de l'évaluer même si on ne connait pas $f'$ explicitement. Pour cela on peut utiliser l'approximation

$$ g'(x) \simeq \frac{g(x+\epsilon) - g(x-\epsilon)}{2\epsilon}.$$

On pourra prendre $\epsilon = 10^{-3}$ par défaut.

a) Ecrire une fonction qui prend en argument une fonction, un point $x$, et $\epsilon$, puis renvoie la valeur approchée de la dérivée en ce point. Comparer la formule exacte de $f'$ à celle calculé avec votre fonction.¶

In [12]:

def derive(f, x, epsilon=1e-3):
    return (f(x + epsilon) - f(x - epsilon)) / (2 * epsilon)

In [13]:

derive(f, x);

In [14]:

plt.plot(x,df(x), label="vraie")
plt.plot(x,derive(f, x), label="approx")

Out[14]:

[<matplotlib.lines.Line2D at 0x7ad4132a5400>]

On ne voit pas la différence.

b) Ecrire une fonction qui prend en argument une fonction et le point de départ $x_0$, et qui renvoie la suite de point $x_0,\dots, x_N$.¶

In [15]:

def newton(g, x0, N, epsilon = 1e-3):
    '''
    g = fct pour resoudre g(a) = 0
    x0 = point de départ
    N = nombre d'itération
    epsilon = precision pour approcher la dérivée de g
    '''
    
    X = np.zeros(N+1) # pour enregistrer toutes les valeurs de la suites
    x_tmp = x0
    X[0] = x0
    
    for j in range(1,N+1): # nous donne tous les indices entre 1 et N
        x_tmp = x_tmp -  g(x_tmp) / derive(g, x_tmp, epsilon)
        X[j] = x_tmp
        
    return X

c) Illustrer graphiquement la méthode de Newton pour $f'$. On pourra prendre $N=5$ et $x_0=1.8$¶

In [16]:

N = 4
x0 = 1.8
X = newton(df, x0, N)

In [17]:

plt.plot(x, df(x), label="f'")
plt.plot(X, df(X), 'o', label="f'")

#illustration des tangentes sur 3 points.

plt.hlines(0, x.min(), x.max(), color="black") #droite y=0

plt.plot([X[0], X[1]], [df(X[0]), 0])
plt.plot([X[1], X[2]], [df(X[1]), 0])
plt.plot([X[2], X[3]], [df(X[2]), 0])

plt.xlim(-1,2);

In [18]:

X[-1]

Out[18]:

np.float64(-0.001199050661786799)

d) Jouer avec $x_0$ pour illustrer la convergence ou non de la méthode de Newton.¶

In [19]:

N = 4
x0 = 1.9
X = newton(df, x0, N)

In [20]:

plt.plot(x, df(x), label="f'")
plt.plot(X, df(X), 'o', label="f'")

plt.hlines(0, x.min(), x.max(), color="black") #droite y=0

plt.plot([X[0], X[1]], [df(X[0]), 0])
plt.plot([X[1], X[2]], [df(X[1]), 0])
plt.plot([X[2], X[3]], [df(X[2]), 0])

Out[20]:

[<matplotlib.lines.Line2D at 0x7ad4133c5af0>]

Ici le point de départ n'est pas approprié, la suite par vers $-\infty$.