{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Approximation des racines d'une équation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Marc Lorenzi - 28 avril 2018" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "from math import *" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams['figure.figsize'] = (8, 8)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Résolution par dichotomie" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $f:[a,b]\\to\\mathbb R$ une fonction continue telle que que $f(a)$ et $f(b)$ sont de signes opposés. Le théorème des valeurs intermédiaires affirme qu'il existe $\\xi\\in[a,b]$ tel que $f(\\xi)=0$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Posons $a_0=a$ et $b_0=b$. Soit $c=\\frac{a+b} 2$. Si $f(a)$ et $f(c)$ sont de signes contraires alors il existe une racine de $f$ sur $[a,c]$. Posons alors $a_1=a_0$ et $b_1=c$. Sinon, il existe une racine de $f$ sur $[c,b]$ et on pose $a_1=c$ et $b_1=b_0$. La fonction $f$ possède donc une racine $\\xi_1$ sur l'intervalle $[a_1,b_1]$ et $f(a_1)$ et $f(b_1)$ sont de signes contraires. \n", "\n", "\n", "On peut évidemment recommencer à loisir cette opération ... On définit ainsi une suite $[a_n,b_n]$ de segments emboîtés tels que $f$ possède une racine $\\xi_n$ sur l'intervalle $[a_n,b_n]$ et $f(a_n)$ et $f(b_n)$ sont de signes contraires.\n", "\n", "De plus, $b_{n+1}-a_{n+1}=\\frac 1 2 (b_n-a_n)$ donc $b_n-a_n=\\frac 1 {2^n} (b_0-a_0)$. La longueur des segments tend vers 0 lorsque $n$ tend vers l'infini. Par le théorème des segments emboîtés, les suites $(a_n)$ et $(b_n)$ sont adjacentes. Notons $\\xi$ leur limite commune. Pour tout $n$, $a_n\\le\\xi_n\\le b_n$ donc, par le théorème d'encadrement des limites, $\\xi_n$ tend vers $\\xi$. Or, $f(\\xi_n)=0$ et $f$ est continue en $\\xi$. Donc $f(\\xi)=0$ : $\\xi$ est une racine de $f$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction `racine` prend en paramètres une fonction $f$, deux réels $a$ et $b$ et un réel $\\varepsilon > 0$. Elle suppose que $f$ est continue sur $[a,b]$ et que $f(a)$ et $f(b)$ sont de signes contraires. Tant que $|b-a|>\\varepsilon$, elle coupe l'intervalle $[a, b]$ en son milieu $c$ et elle recommence sur le demi-intervalle adéquat.\n", "\n", "La fonction renvoie une approximation d'une racine de $f$ ainsi que le nombre d'itérations effectuées pour la trouver." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def racine(f, a, b, epsilon):\n", " c = (a + b) / 2\n", " cpt = 0\n", " while abs(b - a) > epsilon:\n", " cpt = cpt + 1\n", " c = (a + b) / 2\n", " if f(a) * f(c) <= 0: b = c\n", " else: a = c\n", " return (c, cpt)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Testons sur $f:x\\mapsto x^2-2$ entre 1 et 2. Ce sera notre exemple fétiche jusqu'à la fin du notebook mais rien ne vous empêche d'essayer plein d'autres fonctions." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "racine(lambda x: x ** 2 - 2, 1, 2, 1e-15)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pourquoi 50 itérations ? Eh bien, à la $k$-ième itération, la longueur de l'intervalle $[a_n,b_n]$ est $\\frac{b-a}{2^k}$. Le nombre d'itérations est donc le plus petit entier $k$ tel que $\\frac{b-a}{2^k}\\le \\varepsilon$, ou encore $k\\ge \\lg \\frac{b-a}{\\varepsilon}$, $\\lg$ étant le logarithme en base 2." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "log(1 / 10 ** (-15)) / log(2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Effectivement, 50 est la bonne valeur." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Quelques remarques pour terminer :\n", "\n", "1. La fonction renvoie __une__ racine de $f$, en fin une approximation. Quelle racine ? On ne sait pas ...\n", "\n", "2. La méthode par dichotomie possède une propriété infiniment précieuse : elle renvoie un __encadrement__ d'une racine (à condition de renvoyer $(a, b)$ et pas $c$, évidemment). Ainsi, on maîtrise l'erreur commise." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Exercice__ : Soit $f:[a,b]\\to\\mathbb R$ continue ayant une racine. Soit $E=\\{x\\in[a,b], f(x)=0\\}$. Montrer que $E$ a un plus petit élément." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Points fixes d'une fonction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.1 Introduction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $f : I \\to I$ où $I$ est un intervalle de $\\mathbb R$. Un __point fixe__ de $f$ est un réel $x\\in I$ tel que $f(x)=x$. Donnons nous $x_0\\in I$ et considérons la suite $(x_n)_{n\\ge 0}$ définie par récurrence par $x_{n+1}=f(x_n)$. Sous certaines conditions, lorsque $n$ tend vers l'infini $x_n$ tend vers un point fixe de $f$. Nous expliciterons quelques conditions suffisantes dans ce qui va suivre. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2 Recherche d'un point fixe" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction `point_fixe` prend en paramètres une fonction $f$, un réel $x_0$ et un réel $\\varepsilon>0$. Tant que $|f(x)-x|>\\varepsilon$, elle remplace $x$ par $f(x)$. Toutefois, la fonction abandonne après 100 itérations, même si la condition d'arrêt n'est pas remplie.\n", "\n", "Ce test d'arrêt est plus que discutable, voire absurde. Mais faute de mieux on s'en contentera. Il est bon toutefois de noter que, sous certaines conditions, ce test est cohérent." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def point_fixe(f, x0, epsilon):\n", " x = x0\n", " cpt = 0\n", " while abs(f(x) - x) > epsilon:\n", " x = f(x)\n", " cpt = cpt + 1\n", " if cpt >= 100: break\n", " return (x, cpt)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "point_fixe(lambda x: cos(x), 0, 1e-15)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Un point fixe de la fonction cosinus est atteint \"à $10^{-15}$ près\" en 87 itérations." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.3 Visualiser les itérations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Petite fonction renvoyant une subdivision régulière du segment $[a,b]$. Bien pratique." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def subdivision(a, b, n):\n", " d = (b - a) / n\n", " return [a + k * d for k in range(n)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction `escalier` dessine la courbe d'une fonction $f$ sur le segment $[a,b]$, ainsi que la droite d'équation $y=x$. Puis elle visualise les différents termes de la suite de premier terme $x_0$ définie par $x_{n+1}=f(x_n)$. Vite écrite, mal écrite. On peut clairement améliorer la présentation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def escalier(f, a, b, c, d, x0, n):\n", " xs = subdivision(a, b, 200)\n", " ys = [f(x) for x in xs]\n", " pad = 0.\n", " plt.axis([a - pad, b + pad, c - pad, d + pad])\n", " plt.plot(xs, ys, 'b')\n", " plt.plot(xs, xs, 'k')\n", " x = x0\n", " coul = 'r'\n", " plt.plot([x0, x0], [0, f(x0)], coul)\n", " for k in range(n):\n", " plt.plot([x, f(x)], [f(x), f(x)], coul)\n", " plt.plot([f(x), f(x)], [f(x), f(f(x))], coul)\n", " x = f(x)\n", " plt.grid()\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "escalier(cos, 0, 1, 0, 1, 0.2, 10)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "escalier(lambda x: sin(x), 0, 3, 0, 1.5, 1.8, 10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.4 Une condition suffisante" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "À quelles conditions (suffisantes mais non triviales) une fonction possède-t-elle un point fixe ? L'ensemble $I$ ci-dessous est un intervalle fermé de $\\mathbb R$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Proposition__ : Soit $f:I\\to I$ une application $k$-lispchitzienne où $k\\in[0,1[$. Alors $f$ admet un unique point fixe $\\ell$. Mieux, pour tout $x_0\\in I$, la suite définie par $x_{n+1}=f(x_n)$ converge vers le point fixe $\\ell$. \n", "\n", "__Démonstration__ : Soient $\\ell_1$ et $\\ell_2$ deux points fixes de $f$. On a $|\\ell_1-\\ell_2|=|f(\\ell_1)-f(\\ell_2)|\\le k|\\ell_1-\\ell_2|$ d'où $(1-k)|\\ell_1-\\ell_2|\\le 0$. Mais $1-k>0$, donc $|\\ell_1-\\ell_2|\\le 0$. D'où $\\ell_1=\\ell_2$ et l'unicité.\n", "\n", "Passons à l'existence. On a pour tout entier $n\\ge 1$, $|x_{n+1}-x_n|=|f(x_n)-f(x_{n-1})|\\le k |x_n-x_{n-1}|$. On déduit par une récurrence facile que pour tout $n\\ge 0$, $|x_{n+1}-x_n|\\le k^n|x_1-x_0|$. Le terme général de la série $\\sum(x_{n+1}-x_n)$ est donc majoré en valeur absolue par le terme général d'une série géométrique convergente. Ainsi, la série $\\sum(x_{n+1}-x_n)$ est absolument convergente, donc convergente. Mais la $(n-1)$ième somme partielle de cette série est $S_{n-1}=\\sum_{k=0}^{n-1}(x_{k+1}-x_k)=x_n-x_0$ par téléscopage. On en déduit que $x_n$ converge vers un réel $\\ell$ lorsque $n$ tend vers l'infini. Comme l'intervalle $I$ est fermé, $\\ell\\in I$. Prenons l'égalité $x_{n+1} = f(x_n)$, valable pour tout entier $n$. Un passage à la limite ($f$ est lipschitzienne sur $I$, donc continue en $\\ell$) donne $\\ell = f(\\ell)$. On a donc l'existence et la limite adéquate. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Remarque__ : On a pour tout entier $n\\ge 1$, $|x_{n+1}-\\ell|=|f(x_n)-f(\\ell)|\\le k |x_n-\\ell|$. On déduit par une récurrence facile que pour tout $n\\ge 0$, $|x_{n}-\\ell|\\le k^n|x_0-\\ell|$. On a donc une majoration explicite de l'écart entre le $n$ième terme de la suite, $x_n$, et sa limite $\\ell$. On constate également que plus le coefficient de Lischitz $k$ est petit, plus la suite converge vite vers sa limite." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.5 Application à la recherche de racines" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $f:I\\to \\mathbb R$. Pour trouver une racine de $f$, une idée consiste à considérer la fonction $g:x\\mapsto x - f(x)$. En effet, pour tout $x\\in I$, on a $f(x)=0$ si et seulement si $g(x)=x$. Pour des raisons qui apparaîtront plus tard, on va plutôt prendre un réel $\\mu$ et s'intéresser à la fonction $g : x\\mapsto x - \\mu f(x)$. En effet, le paragraphe précédent nous dit qu'en abaissant le coefficient de Lispchitz de $g$ la suite $(x_n)$ tendra plus vite vers sa limite. En glissant un paramètre dans la fonciton $g$, on peut peut-être maîtriser le coefficient de Lipschitz.\n", "\n", "Regardons un peu ce qui se passe avec différentes valeurs de $\\mu$. On va s'apercevoir que pour certains $\\mu$ il y a divergence de l'algorithme, et pour d'autres convergence. De plus, selon le choix de $\\mu$, la convergence est plus ou moins rapide. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def racine2(f, x0, mu, epsilon):\n", " g = lambda x: x - mu * f(x)\n", " return point_fixe(g, x0, epsilon)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "racine2(lambda x: x ** 2 - 2, 2, 0.2, 1e-10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Prenons des valeurs de $\\mu$ entre -1 et 1 et regardons le nombre d'itérations nécessaires à la terminaison de notre algorithme." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "s = subdivision(0, 1, 1000)\n", "cs = [racine2(lambda x: x ** 2 - 2, 2, mu, 1e-10)[1] for mu in s]\n", "plt.plot(s, cs, 'k')\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Trop beau. Faisons un zoom dans la zone ou le nombre d'itérations est minimal ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "s = subdivision(0.31, 0.4, 1000)\n", "cs = [racine2(lambda x: x ** 2 - 2, 2, mu, 1e-10)[1] for mu in s]\n", "plt.plot(s, cs, 'k')\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Une valeur de $\\mu$ quelque part entre 0.35 et 0.36 a l'air optimale. Il s'avère que cette valeur est $\\frac{1}{2\\sqrt 2}$. Mais pourquoi ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "racine2(lambda x: x ** 2 - 2, 2, 1 / (2 * sqrt(2)), 1e-10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "L'inégalité des accroissements finis nous dit que si $|g'|$ est majorée par un réel $k$ alors $g$ est $k$-lipschitzienne. Pour accélérer la convergence de nos suites il suffit donc de minimiser $g'$, au moins au voisinage de $\\ell$. Si on prend $g'(\\ell)=0$, $|g'|$ sera $k$-lipschitzienne avec $k$ aussi petit qu'on veut pourvu qu'on se restreingne à un voisinage suffisamment petit de $\\ell$. \n", " \n", "Pour tout réel $x$, $g'(x)=1-\\mu f'(x)$. Le réel $\\mu =\\frac 1 {f'(\\ell)}$ paraît donc très prometteur." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dans l'exemple que nous avons regardé, $f(x)=x^2-2$, $\\ell=\\sqrt 2$ et, justement, $f'(\\ell)=\\frac 1 {2\\sqrt 2}$ !" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Il semblerait que l'on approche du Graal : la suite à considérer serait donc $x_{n+1}=g(x_n)$ où $g(x)=x-\\frac{f(x)}{f'(\\ell)}$, $\\ell$ étant la limite de la suite. Euh oui, mais on a un gros souci : on ne connaît pas $\\ell$ ! C'est toujours comme ça avec le Graal, chaque fois qu'on croit le tenir il nous échappe :-).\n", "\n", "Soyons rusés, en prenant un paramètre $\\mu$ variable qui tend vers $\\frac 1 {f'(\\ell)}$ lorsque $n$ tend vers l'infini. Lequel ? eh bien $\\mu = \\frac 1 {f'(x_n)}$. LA suite à étudier est donc la suite\n", "\n", "$$x_{n+1}= x_n - \\frac{f(x_n)}{f'(x_n)}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Il s'avère que cette suite est justement à la base d'un algorithme connu depuis des siècles, l'algorithme de Newton-Raphson. L'idée a été ébauchée par Newton et formalisée par Raphson. C'est le but du prochain paragraphe." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. L'algorithme de Newton-Raphson" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.1 Dérivation approchée" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici une fonction renvoyant la valeur approchée de la dérivée d'une fonction $f$ en un point $x$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def deriv(f, x, h=1e-5):\n", " return (f(x + h) - f(x - h)) / (2 * h)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "deriv(lambda x: exp(x), 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.2 La méthode de Newton" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $f : I\\to\\mathbb R$ une fonction dérivable en point $x\\in I$. L'équation de la tangente à la courbe représentative de $f$ au point $A=(x, f(x))$ est\n", "\n", "$$(T) Y = f(x) + (X-x)f'(x)$$\n", "\n", "Quelle est l'intersection de $T$ et de l'axe $Ox$ ? Il suffit d'annuler $Y$ dans l'équation de $T$. Il vient\n", "\n", "$$ X = x - \\frac {f(x)}{f'(x)}$$\n", "\n", "Surprise ... La suite récurrente que nous nous proposons d'étudier a donc une interprétation géométrique simple. Pour passer de $x_n$ à $x_{n+1}$, on trace la tangente à la courbe au point d'abscisse $x_n$. L'intersection de cette tangente avec l'axe $Ox$ a alors pour abscisse $x_{n+1}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici un dessin simpliste." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xs = subdivision(1.25, 2, 1000)\n", "f = lambda x: x ** 2 - 2\n", "ys = [f(x) for x in xs]\n", "plt.plot(xs, ys, 'k')\n", "plt.plot(xs, [0 for x in xs], 'b')\n", "x = 1.8\n", "x1 = x - f(x) / deriv(f, x) \n", "plt.plot([x, x1], [f(x), 0], 'r')\n", "plt.plot([x],[f(x)],'ko')\n", "plt.plot([x1],[0],'ko')\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Newton se contentait de prendre pour valeur approchée de la racine de $f$ l'intersection du segment bleu et du segment rouge. Nous devons à Raphson l'idée que l'on pouvait recommencer avec cette nouvelle valeur. Enfin, oui, j'avoue, Newton a quand même un peu itéré aussi :-) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.3 L'algorithme" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Il est immédiat. Il n'y a qu'à utiliser ce que nous avons déjà fait." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def newton(f, x0, epsilon):\n", " g = lambda x: x - f(x) / deriv(f, x)\n", " return point_fixe(g, x0, epsilon)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici $\\sqrt 2$ en 4 itérations avec 10 chiffres après la virgule." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "newton(lambda x: x ** 2 - 2, 2, 1e-10 )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Et voici $\\pi$ en 3 itérations avec 15 chiffres." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "newton(lambda x:sin(x), 3, 1e-15)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour bien mesurer à quel point la méthode de Newton est efficace, travailler avec 10 chiffres après la virgule ne suffit pas." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sympy import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculons $\\sqrt 2$ ... avec 2000 chiffres après la virugule ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prec = 2000\n", "newton(lambda x: x ** 2 - 2, N(2, prec + 1), S(10) ** (- prec))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Eh oui. En 12 itérations. L'algorithme de Newton-Raphson est TRÈS performant. Vous n'avez pas confiance, vous n'y croyez pas ? Comparons le résultat $r$ obtenu par Newton à une valeur approchée de $\\sqrt 2$, $x$, demandée à `sympy`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "prec = 2000\n", "r, c = newton(lambda x: x ** 2 - 2, N(2, prec + 1), S(10) ** (- prec))\n", "x = N(sqrt(2), prec)\n", "print(r - x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Demandez 10000 chiffres si cela vous tente ... Il faudra peut-être attendre 1 ou 2 secondes mais c'est normal. Faire une simple division avec autant de chiffres prend quand même un certain temps, et `sympy` n'a pas été construit pour cela. Je l'ai lancé avec 100000, mais je n'ai pas osé le million. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.4 Et sur $\\mathbb C$ ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Est-ce que la méthode de Newton permet de trouver une racine __non réelle__ de $x^3-1$ ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Exécutons `newton` avec $x_0=i$ ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "newton(lambda x: x ** 3 - 1, 0+1J, 1e-15)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Super, ça converge vers $j$, la racine cubique de 1.\n", "\n", "Maintenant ouvrons la boîte de Pandore. Il y a TROIS racines cubiques de 1. Quelles sont les valeurs de $x_0$ pour lesquelles la suite converge vers 1 ? Vers $j$ ? Vers $j^2$ ? Ci-dessous, on affiche en blanc les points $x_0\\in\\mathbb C$ du carré $[-1,1]\\times[-1,1]$ tels que la suite converge vers 1. Remarquer que pour $f(x)=x^3$, on a $x-\\frac{f(x)}{f'(x)}=\\frac{2x^3+1}{3x^2}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ne vous inquiétez pas pour la ligne ci-dessous. Vous me remercierez si vous faites l'exercice posé un peu plus bas ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from math import sqrt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def attracteur():\n", " N = 300\n", " n = 3\n", " f = lambda w : ((n - 1) * w ** n + 1) / (n * w ** (n - 1))\n", " #f = lambda w: ((2 * w ** 3 + 1) / (3 * w ** 2))\n", " A = N * [None]\n", " for i in range(N): A[i] = N * [0]\n", " for i in range(N):\n", " for j in range(N):\n", " x = -1 + 2 * i / N\n", " y = -1 + 2 * j / N\n", " z = x + y * 1j\n", " for k in range(32):\n", " if abs(z) > 1e-10: z = f(z)\n", " if abs(z - 1) < 1e-10: A[j][i] = 1\n", " plt.imshow(A, interpolation='bicubic', cmap='hot', origin='lower')\n", " plt.savefig('newton.png')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "attracteur()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Exercice__ : Modifier 1 ligne dans le code ci-dessus pour afficher en blanc les points pour lesquels la suite converge vers $j$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "L'ensemble cherché est __TRÈS__ compliqué, donc très intéressant. Imitons Pandore et refermons la boîte en y laissant l'espérance ... qu'un jour nous pourrons comprendre ce dessin." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4.5 La théorie" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 4.5.1 Hypothèses" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plaçons nous dans un cas particulier : $f:I\\to\\mathbb R$ est de classe $\\mathcal C^2$. $\\ell\\in I$ est une racine de $f$ et n'est pas la borne supérieure de $I$. On suppose $f''>0$, $f'(\\ell)>0$. Et on prend $x_0>\\ell$. Ça a l'air très particulier et ça l'est ... mais l'exposé sera plus clair.\n", "\n", "On considère la suite définie par la relation de récurrence $x_{n=1}=g(x_n)$ où $g(x)=x-\\frac{f(x)}{f'x)}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 4.5.2 Convergence" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Proposition__ : la suite $(x_n)$ est décroissante et minorée par $\\ell$.\n", "\n", "__Démonstration__ : Soit $n\\in\\mathbb N$. Supposons $\\ell\\ell$ et $x_{n+1} < x_n$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Proposition__ : $x_n \\to \\ell$ lorsque $n$ tend vers l'infini.\n", "\n", "__Démonstration__ : $(x_n)$ est décroissante et minorée, donc converge vers $\\ell'\\in\\mathbb R$. Comme les $x_n$ appartiennent à $[\\ell, x_0]$, la limite $\\ell'$ de la suite aussi et cette limite vérifie $g(\\ell')=\\ell'$, c'est à dire $f(\\ell')=0$. Mais seul $\\ell$ vérifie cela sur $[\\ell, x_0]$. Donc $\\ell'=\\ell$. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### 4.5.3 Vitesse de convergence" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Appliquons maintenant le théorème des accroissements finis à $g$ entre $\\ell$ et $x_{n+1}$. Il existe $c_n$ (pas le même que ci-dessus )), $\\ell