Résolution numérique d'équations sur les réels¶
Zéro d'une fonction¶
Soit $E$ une partie non vide de $\mathbb{R}$
Définition 1 :
On appelle zéro d'une application $f : E \rightarrow \mathbb{R}$ tout $x_r\in E$ vérifiant $f(x_r) = 0$
Définition 2 :
On appelle point fixe d'une application $g : E \rightarrow\mathbb{R}$ tout $x_r\in E$ vérifiant $g(x_r)=x_r$
Remarque :
Il est équivalent de s'intéresser à la résolution d'équations $f(x)=0$ ou à la recherche de points fixes. L'équation $f(x)=0$ est équivalente à $f(x) + x = x$. On pose $g(x) = f(x) + x$ donc on obtient $g(x) = x$
Représentation graphique¶
Les zéros d'une fonction¶
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
def f(x):
return x*np.sin(np.pi/x)
x = np.linspace(0.05, 0.5, 300)
y = f(x)
plt.grid()
plt.plot(x, y)
plt.show()
Les points fixes
def f(x):
return x**2 - 2*x -3 # les zéros sont -1 et 3
def g(x):
return f(x) + x
x = np.linspace(-2, 5, 300)
yf = f(x)
yg = g(x)
plt.grid()
plt.plot(-1, 0, 'o')
plt.plot(3, 0, 'o')
plt.plot(x, yf, 'b-', label='$f(x) = x^2-2x-3$')
plt.plot(x, yg, 'r--', label="$g(x) =f(x)+x$")
plt.plot(x, x, 'g', label='$h(x) = x$')
plt.legend(loc='upper center', shadow=True, fontsize='x-large')
plt.show()
Méthode de la dichotomie¶
Théorème : Soient $a, b \in \mathbb{R}, a<b$. Soit $f:[a, b]\rightarrow\mathbb{R}$ une application continue telle que $f(a)$ et $f(b)$ sont non nuls et de signes opposés. Alors il existe $x_r\in]a, b[$ tel que $f(x_r)=0$. Si de plus $f$ est strictement monotone, alors $x_r$ est unique.
Principe¶
La méthode d'approximation de $x_r$ par la dichotomie est la suivante :
On pose $g_0=a$ et $d_0=b$, soit $n\in \mathbb{N}$ et $f(g_n)f(d_n) \leq 0$
On pose $m=\dfrac{g_n+d_n}{2}$
Si $f(g_n)f(m) \leq 0$, on pose $g_{n+1}=g$ et $d_{n+1}=m$
Sinon on pose $g_{n+1}=m$ et $d_{n+1} = d$
Ces deux suites convergent vers une solution $x_r$ de l'équation $f(x)=0$
Critère d'arrêt¶
On se donne $\epsilon > 0$ et on s'arrète dès que $|d_n - g_n| < \epsilon$, attention un critère d'arrêt comme $f(m)=0$ ou $(d_n - g_n) = 0$ serait une très mauvaise idée
Implantation de l'algorithme de dichotomie¶
from math import sqrt
def dichotomie(f, a, b, epsilon):
g, d = a, b
fg, fd = f(g), f(d)
cpt = 0
while abs(d - g) > epsilon:
m = (g + d) / 2.
cpt +=1
fm = f(m)
if fm * fg < 0:
d, fd = m, fm
else:
g, fg = m, fm
return g, d, cpt # On renvoie les deux bornes de notre encadrement g < xr < d
Déterminer sur votre feuille les 4 premiers itérés dans l’intervalle [1, 3] avec la fonction $f(x) = x^2 - 2$
Avant la première itération $g=1\quad fg=-1\quad d=3\quad fd=7$
$\begin{array}{|l|c|c|c|c|c|c|c|}\hline
Itération & |d-g| & m & fm & fm\times fg & g & fg & d & fd\\\hline
1 & 2 & 2 & 2 & négatif & 1 & -1 & 2 & 2\\\hline
2 & 1 &1.5 & 0.25 & négatif & 1 & -1 & 1.5 & 0.25\\\hline
3 & 0.5 &1.25 & -0.4375 & positif & 1.25 & -0.4375 & 1.5 & 0.25\\\hline
4 & 0.25 &1.375 & -0.109375 & positif & 1.375 & -0.109375 & 1.5 & 0.25\\\hline
5 & 0.125 &1.4375 & 0.06640625 & négatif & 1.375 & -0.109375 & 1.4375 & 0.06640625\\\hline
\end{array}
$
sqrt(2), dichotomie(lambda x: x**2 - 2, 1, 3, 1e-6)
Invariant de boucle¶
- Terminaison
L'intervalle $[a, b]$ est divisé par 2 à chaque itération, au bout de la $k$ième itération (avec $k\in \mathbb{N}$), on a $d-g = \dfrac{b-a}{2^k}$ donc pour $k$ suffisamment grand $d-g < \epsilon$ ce qui termine la boucle
- Sémantique
Pour être sûr de trouver une valeur numérique approchée de $x_r$ il faut qu' à la fin de chaque itération l'assertion $P(n) : f(g_n)f(d_n) < 0$ soit vraie
Preuve : Raisonnement par récurrence
Initialisation : Nous avons vérifié dans la question précédente que $P(0)$ est vraie
Hérédité : Supposons $P(n)$ vraie pour un certain rang $n$ $(n\in \mathbb{N})$ et montrons que $P(n+1)$ est vraie
- Soient $\quad m=\dfrac{g_n+d_n}{2}\quad$ et $\quad f(m)f(g_n) < 0$
donc d'après l'algorithme de dichotomie : $d_{n+1} = m$ et $g_{n+1} = g_n$ donc $f(m)f(g_n) = f(g_{n+1})f(d_{n+1}) < 0$ donc P(n+1) : vraie
- Soient $\quad m=\dfrac{g_n+d_n}{2}\quad$ et $\quad f(m)f(g) \geq 0\quad$ donc $f(m)$ et $f(g_n)$ de même signe, on sait par hypothèse de récurrence que $f(g_n)$ et $f(d_n)$ sont de signe opposé donc $f(m)$ et $f(d_n)$ sont aussi de signe opposé.
donc d'après l'algorithme de dichotomie : $d_{n+1} = d_n$ et $g_{n+1} = m$, donc $f(g_{n+1})f(d_{n+1}) < 0$ donc P(n+1) : vraie
Ce qui termine la récurrence et avec la condition d'arrêt fait la preuve de notre algorithme.
Complexité¶
Soit $n\in \mathbb{N}$ le nombre d'itérations. À partir de l'étude de la terminaison on peut écrire: $$b_n - a_n = \frac{b - a}{2^n} \quad \text{donc quand}\quad (b_n-a_n) \leq \epsilon \quad\text{on a} \quad \frac{b - a}{2^n} \leq \epsilon \quad\text{donc}\quad 2^n \geq \frac{b - a}{\epsilon} \quad\text{donc}\quad n \geq log_2\left(\frac{b-a}{\epsilon}\right)$$
Méthode de Newton¶
Critères d'arrêt¶
Nous avons le choix entre deux types de critères d'arrêt pour interrompre le processus itératif d'approximation de la racine de la fonction $f$ (appelée $x_r$). On note $\epsilon$ la tolérance fixée par le calcul approché de $x_r$ et $e_n = x_r - x_n$ l'erreur absolue. On rapelle que la suite récurrente utilisée est de la forme. $$x_{n+1} = x_n-\frac{f(x_n)}{f'(x_n)}$$
- Contrôle du résidu les itérations s'achèvent dès que $|f(x_n)| < \epsilon$
- $f'(x_n) \simeq 1$ donc $|x_r - x_n| \simeq \epsilon$
- $f'(x_n) \ll 1$ danger car $|x_r - x_n|$ peut-être grand par rapport à $\epsilon$
- $f'(x_n) \gg 1$ alors $|x_r - x_n| \ll \epsilon$ et le test est dit trop restrictif
- Contrôle de l'incrément les itérations s'achèvent dès que $|x_{n+1} - x_n| < \epsilon$
On note $e_{n+1} = x_r - x_{n+1}$ et soit $\phi$ une fonction définie sur $[x_r; x_n]$ continue et dérivable sur $[x_r; x_n]$ tel que $\phi(x) = x - \dfrac{f(x)}{f'(x)}$
Soit $x_r \in \mathbb{R}$ la racine de $f$ donc $f(x_r) = 0$ $(f'(x_r)\neq0)$ donc $\phi(x_r) = x_r$ donc $x_r$ point fixe de la fonction $\phi$. Approcher les zéros de $f$ se ramène donc au problème de la détermination des points fixes de la fonction $\phi$
donc $e_{n+1} = \phi(x_r) - \phi(x_n)$ et d'après l'égalité des accroissements finis il existe $c \in [x_r; x_n]$ (avec $x_r < x_n$) tel que $e_{n+1}=x_r - x_{n+1}=\phi'(c)(x_r-x_n)$
Or $x_{n+1}-x_n = (x_r-x_n) + (x_{n+1}-x_r) = (x_r-x_n) - (x_r - x_{n+1})$
et $(x_r-x_n) + (x_r - x_{n+1}) = e_n - e_{n+1} = e_n - \phi '(c)(x_r-x_n) = e_n - \phi '(c)\times e_n = e_n(1 - \phi '(c))$
donc $$x_{n+1} - x_n = e_n(1-\phi '(c))$$
Soit $n \in \mathbb{N}$ un nombre d'itérations suffisamment grand pour que $\phi '(c) \simeq \phi '(x_r)$ et on sait que $\phi '(x_r) = 0$
$$\text{donc}\quad x_{n+1} - x_n \simeq e_n < \epsilon$$
L'erreur que l'on commet lorsque l'on adopte ce critère est donc inférieure à la tolérance fixée.
Implantation de l'algorithme de Newton¶
def newton(f, g, a, epsilon):
x, cpt = a, 0 # Premier terme de la suite récurrente x0
y = x - f(x)/g(x) # On calcule x(n+1)
while abs(y-x)>epsilon: # Critère d'arrêt |x(n+1) - x(n) < epsilon|
x = y # On conserve x(n) avant de calculer x(n+1) pour tester la condition d'arrêt
y = y - f(y) / g(y)
cpt += 1
print y
return y, cpt
Avec Newton la convergence est quadratique, le nombre de décimales exactes double (au moins) à chaque étape. Attention pour les valeurs supérieurs à 0.2 Tracer la courbe 1/x - 0.7 Que se passe-t-il si l'erreur est grande ?
newton(lambda x: 1./x -7, lambda x: -1./x**2, 0.001, 1e-10)
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
def graph_newton(f, g, a, epsilon, I=[-10, 10]):
x = float(a)
y = x - f(x)/g(x)
X, Y = [x, x, y, y], [0, f(x), 0, f(y)]
while abs(y-x) > epsilon:
x = y
y = y - f(y) / g(y)
X += [y, y]
Y += [0, f(y)]
disc = np.linspace(I[0], I[1], 100)
plt.grid()
plt.plot(disc, f(disc))
plt.plot(X, Y, 'r--*')
plt.savefig("newton.png")
plt.show()
graph_newton(lambda x: x**3-x**2+8*x-8, lambda x: 3*x**2-2*x+8, 6, 1e-7, [-2, 6])
graph_newton(lambda x: x**2-2, lambda x: 2*x, 10, 1e-10, [-4, 10])
graph_newton(lambda x: 1-x*np.log(x), lambda x: -(np.log(x)+1), 5, 1e-10, [0, 5])
graph_newton(lambda x: np.cos(x), lambda x: - np.sin(x) , 0.5, 1e-10, [-1, 3])
Pour la fonction $f : x \mapsto xe^{-x}$ la racine est $0$. La méthode de Newton converge pour $x_0 \in [0,1]$ et diverge pour $x_0 > 1$
graph_newton(lambda x: x*np.exp(-x), lambda x: np.exp(-x)*(1-x) , 0.5, 1e-10, [-1, 1.5])
graph_newton(lambda x: x*np.exp(-x), lambda x: np.exp(-x)*(1-x) , 1.5, 1e-10, [-1, 5])
Une fonction qui donne un cycle $f : x \mapsto (x+2)(x^2+0.46)$
graph_newton(lambda x: (x+2)*(x**2+0.46), lambda x: 3*x**2+4*x+0.46, -0.6, 1, [-1, 1])
Peut-il y avoir un cycle avec la fonction $f: x\mapsto \sin x$ ?
La fonction $f:x \mapsto \sin x$ est impaire donc il faut qu'il ait un $x_n \in I$, tel que $x_{n+1} = -x_n$ la méthode de Newton permet alors de construire un parallèlogramme avec les tangentes.
$x_{n+1} = x_n - \dfrac{f(x_n)}{f'(x_n)}$ donc $-x_n = x_n - \dfrac{\sin x_n}{\cos x_n}$ donc $2x_n - \tan x_n = 0$ donc on doit résoudre
très difficile de résoudre cette équation à la main, donc on utilise l'algorithme de newton (et voila la boucle est bouclée, si vous n'étiez pas encore persuadé de son utilité il n'a pas fallu attendre longtemps)
newton(lambda x: np.tan(x)-2*x, lambda x: 1/np.cos(x)**2 - 2, 1.16, 1e-15)
Une solution à $10^{-15}$ près est $x_r = 1.1655611852072112$
graph_newton(lambda x: np.sin(x), lambda x: np.cos(x), 1.1655611852072112, 1e-7, [-2, 2])
On remarquera que la valeur numérique $x_r$ n'est qu'une valeur approchée puisqu'on finit par sortir du cycle et converger vers 0.
Le 2-cycle précédent est-il unique ?
Soit la fonction $g : x \mapsto x - \tan(x)$ avec $x\in [0, \frac{\pi}{2}[$
L'ensemble des valeurs de $x$ menant à un 2-cycle avec la méthode de Newton lors de la recherche des zéros de la fonction sinus est donné par les solutions de l'équation :
On pose $h_1 : x \mapsto \tan(x)$ et $h_2 : x \mapsto -\tan(x - \tan(x))$. Pour résoudre l'équation précédente on cherche s'il y en a les points d'intersections des courbes des fonctions $h_1$ et $h_2$ pour $x\in [0, \frac{\pi}{2}[$ Chaque point d'intersection étant le point de départ d'un 2-cycle.
def f_h1(x):
y = np.tan(x)
y[np.diff(y) >= 0.5] = np.nan # si y[n+1] - y[n] >= 0.5 alors y[n] = not a number détection d'asymptote
return y
def f_h2(x):
y = -np.tan(x - np.tan(x))
y[np.diff(y) >= 20] = np.nan
return y
x = np.linspace(1, np.pi/2, 1e5)
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
ax.axis([1, 1.55, -1, 20])
ax.plot(x, f_h2(x))
ax.plot(x, f_h1(x), 'r')
plt.grid()
plt.show()
$h_1$ en rouge intersecte $h_2$ un grand nombre de fois. Il y a donc de nombreuses solutions sur l'intervalle $[0, \frac{\pi}{2}[$ à l'équation : $\tan(x) + \tan(x-\tan(x)) = 0$. En conclusion il y a de nombreux 2-cycles avec la méhode de Newton lors de la recherche de zéros de la fonction $sinus$. Il en est de même avec les k-cycles ($k\in \mathbb{N}$).
christophe.casseau@ensam.eu