{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# $x^{x^{x^{x^{\\ldots}}}}=$ ?\n", "\n", "Marc Lorenzi\n", "\n", "15 décembre 2018" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import math\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams['figure.figsize'] = (12, 6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Exponentielle infinie" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.1 Introduction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $x\\in\\mathbb R_+^*$. Quel sens donner à $\\varphi(x)=x^{x^{x^{x^{\\ldots}}}}$ ? Des points de suspensions avec rien au bout nous suggèrent une limite de suite.\n", "\n", "Posons $x^{[0]}=1$ et, pour tout $n\\in\\mathbb N$, $x^{[n+1]}=x^{x^{[n]}}$. On a $x^{[1]}=x$, $x^{[2]}=x^x$, $x^{[3]}=x^{x^x}$. Notons $x^{[\\infty]}$ la limite, si elle existe, de $x^{[n]}$ lorsque $n$ tend vers l'infini.\n", "\n", "On a $x^{[n]}\\to x^{[\\infty]}$. Mais $x^{x^{[n]}}= x^{[n+1]}\\to x^{[\\infty]}$. On en déduit que \n", "\n", "$$x^{x^{[\\infty]}}= x^{[\\infty]}$$\n", "\n", "Posons dorénavant $\\varphi(x)=x^{[\\infty]}$. Nous venons donc de \"définir\" une fonction $\\varphi : \\mathcal D \\to \\mathbb R$ par\n", "\n", "$$x^{\\varphi(x)}=\\varphi(x)$$\n", "\n", "où $\\mathcal D$ est l'ensemble des réels $x>0$ tel que notre suite converge. Pourvu que $\\mathcal D$ soit non vide, sinon ce notebook risque de tourner court :-)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Exercice__ : Montrer que $1\\in\\mathcal D$ et que $\\varphi(1)=1$.\n", "\n", "__Exercice__ : Montrer que pour tout $n\\in\\mathbb N$, $0^{[2n]}=1$ et $0^{[2n+1]}=0$. En déduire que $0\\not\\in\\mathcal D$.\n", "\n", "__Exercice__ : Montrer que pour tout $x>2$, $x\\not\\in\\mathcal D$. Ainsi, $\\mathcal D\\subset ]0, 2[$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici quelques questions auxquelles nous allons répondre :\n", "\n", "- Que vaut précisément $\\mathcal D$ ?\n", "- Comment calculer efficacement pour $x\\in\\mathcal D$ une valeur approchée de $\\varphi(x)$ ? \n", "- Y-a-t-il des réels $x\\in\\mathcal D$ pour lesquels on peut calculer explicitement $\\varphi(x)$ ?\n", "\n", "Nous verrons que la fonction $\\varphi$ est étroitement liée à une autre fonction, la fonction $W$ de Lambert." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.2 Une valeur approchée" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $x>0$. Le réel $\\varphi(x)$ (s'il existe) est un point fixe de la fonction $t\\mapsto x^t$. Considérer une suite $(u_n)$ définie par $u_{n+1}=x^{u_n}$ devrait peut-être permettre d'approcher ce point fixe. Une simple boucle `while` ferait-elle l'affaire ?\n", "\n", "La fonction `inf_exp` ci-dessous prend un réel $x$ en paramètre. Elle renvoie une valeur que l'on espère approchée de $\\varphi(x)$. Le second paramètre, `dbg` permet l'affichage du nombre d'itérations lorsqu'il est mis à `True`. Le troisième paramètre, `tol`, permet de régler la précision du résultat." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def inf_exp(x, dbg=False, tol=1e-15):\n", " cnt = 0\n", " u, v = 1, x\n", " while abs(v - u) > tol and cnt <= 1e6:\n", " cnt += 1\n", " u, v = v, x ** v\n", " if dbg: print(cnt)\n", " return v" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "inf_exp(1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "y = inf_exp(1.2, True)\n", "print(y)\n", "print(1.2 ** y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.3 Quelques tentatives plus ou moins heureuses" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tentative 1 ... On a $\\sqrt 2^2=2$. Se pourrait-il donc que $\\sqrt 2^{[\\infty]}=2$ ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = math.sqrt(2)\n", "y = inf_exp(x, True)\n", "print(y)\n", "print(x ** y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Cela est très prometteur. Mais bien entendu l'approximation précédente ne le __montre__ pas. Elle nous le __suggère__." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tentative 2 ... On a $\\left(e^{\\frac 1 e}\\right)^e=e$. Se pourrait-il que $\\left(e^{\\frac 1 e}\\right)^{[\\infty]}=e$ ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = math.e ** (1 / math.e)\n", "y = inf_exp(x, True, tol=1e-10)\n", "print(y)\n", "print(x ** y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Cela se pourrait. Remarquez la convergence terriblement lente. Il va falloir améliorer cela." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tentative 3 ... On a $\\left(e^{-e}\\right)^{\\frac 1 e}=e^{-1}=\\frac 1 e$. Se pourrait-il que $\\left(e^{-e}\\right)^{[\\infty]}=\\frac 1 e$ ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = math.e ** (- math.e)\n", "y = inf_exp(x, True, tol=1e-5)\n", "print(y)\n", "print(x ** y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Argh, le million d'itérations a été dépassé ! Mais les valeurs numériques renvoyées sont encourageantes. Relançons avec une tolérance de $10^{-2}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = math.e ** (- math.e)\n", "y = inf_exp(x, dbg=True, tol=1e-2)\n", "print(y)\n", "print(x ** y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ultime tentative ... Et si nous prenions $x$ __complexe__ ? Que vaut $i^{i^{i^{\\ldots}}}$ ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "inf_exp(1j, True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Bigre, ça a l'air de converger. Notre tâche va donc être complexe :-). Quels sont les __nombres complexes__ $z$ pour lesquels $z^{[\\infty]}$ existe ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1.4 Les résultats sur $\\mathbb R$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Proposition__ : L'ensemble de définition de $\\varphi$ est $\\mathcal D=[e^{-e},e^{1/e}]$.\n", "\n", "__Démonstration__ : La démonstration est un peu longue, quoique sans difficulté particulière. Nous ne la ferons pas ici. Pour être honnêtes, $D=[e^{-e},e^{1/e}]\\cup\\{-1\\}$.\n", "\n", "Contentons nous de tracer la courbe de la fonction $\\varphi$." ] }, { "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 + 1)]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xs = subdivision(math.e ** (-math.e), math.e ** (1 / math.e), 200)\n", "ys = [inf_exp(x, False, 1e-2) for x in xs]\n", "plt.plot(xs, ys, 'k')\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Remarquez les zigzags au voisinage de $e^{-e}$ ... ceci est dû à la mauvaise convergence des suites au voisinage de ce point. Il va nous falloir remédier à cela.\n", "\n", "Nous allons voir maintenant que $\\varphi$ est liée à une autre fonction, la __fonction de Lambert__." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. La fonction $W$ de Lambert" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.1 C'est quoi ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La fonction de Lambert est la fonction $W$ \"définie\" par $W(xe^x)=x$. Disons pour faire vite que $W$ est la \"réciproque\" de la fonction $f :x\\mapsto xe^x$. Euh, oui, sauf que $f$ n'est pas bijective. Voici le graphe de $f$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xs = subdivision(-5, 1, 200)\n", "ys = [x * math.exp(x) for x in xs]\n", "zs = [0 for x in xs]\n", "plt.plot(xs, ys, 'k')\n", "plt.plot(xs, zs, 'r')\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Une étude rapide de $f$ montre que $f$ est une bijection de $]-\\infty,-1]$ sur $[-\\frac 1 e, 0[$, et aussi une bijection de $[-1,+\\infty[$ sur $[-\\frac 1 e,+\\infty[$.\n", "\n", "Notons $f_1:]-\\infty,-1]\\to[-\\frac 1 e, 0[$ et $f_2:[-1,+\\infty[\\to[-\\frac 1 e,+\\infty[$ les bijections induites par $f$. Nous aurons donc deux fonctions de lambert, réciproques de $f_1$ et $f_2$ :\n", "\n", "__Proposition__ :\n", "- $W_1:[-\\frac 1 e, 0[\\to]-\\infty,-1]$ est une bijection continue strictement décroissante. La fonction $W_1$ est de classe $\\mathcal C^{\\infty}$ sur $[-\\frac 1 e, 0[$. \n", "- $W_2:[-\\frac 1 e, +\\infty[\\to[-1,+\\infty[$ est une bijection continue strictement croissante. La fonction $W_2$ est de classe $\\mathcal C^{\\infty}$ sur $[-\\frac 1 e, +\\infty[$.\n", "\n", "__Démonstration__ : Il suffit d'appliquer les théorèmes généraux sur les fonctions strictement monotones. Que vaut la dérivée de $W_2$ ? Pour tout $x\\in[-\\frac 1 e, +\\infty[$, on a $W'_2(x)=(f_2^{-1})'(x)=\\frac 1 {f_2'\\circ f_2^{-1}(x)}$. Or, $f'_2(t)=(t+1)e^t=e^t+f_2(t)$. Donc, \n", "\n", "$$W'_2(x)=\\frac 1{x+\\exp(W_2(x))}$$\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Quelques valeurs intéressantes prises par la fonction $W_2$ :\n", "\n", "- $W_2(0) = 0$.\n", "- $W_2(-\\frac 1 e) = -1$.\n", "- $W_2(e) = 1$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Exercice__ :\n", "\n", "- Que valent $W'_2(0)$ ? $W'_2(e)$ ?\n", "- Quelle est la limite de $W'_2(x)$ lorsque $x$ tend vers $-\\frac 1 e$ ?\n", "- Quelle est la limite de $W'_2(x)$ lorsque $x$ tend vers $+\\infty$ ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici le graphe de $W_1$ obtenu en \"trichant\". Il n'y a qu'à tracer $f_1$ en échangeant abscisses et ordonnées." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xs = subdivision(-4, -1, 200)\n", "ys = [x * math.exp(x) for x in xs]\n", "plt.plot(ys, xs, 'k')\n", "plt.xlabel('x')\n", "plt.ylabel('W1(x)')\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Et voici le graphe de $W_2$, obtenu par la même technique." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xs = subdivision(-1, 2, 200)\n", "ys = [x * math.exp(x) for x in xs]\n", "plt.plot(ys, xs, 'k')\n", "plt.xlabel('x')\n", "plt.ylabel('W2(x)')\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Évidemment, nous aimerions bien pouvoir calculer $W_i(x)$, $i=1,2$, ou du moins une approximation. Patience, cela va venir." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.2 Le rapport avec les puissances infinies" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Quel est le rapport entre les puissances infinies et la fonctiion $W$ ?\n", "\n", "Rappelons la notation $\\varphi(x)=x^{[\\infty]}$. Soit $x\\in[e^{-e},e^{1/e}]$. On a $x^{\\varphi(x)}=\\varphi(x)$. Passons aux logarithmes pour en déduire $\\varphi(x)\\ln x=\\ln\\varphi(x)$. Si nous posons $\\alpha =\\ln\\varphi(x)$, il vient $e^\\alpha \\ln x = \\alpha$, ou encore $-\\ln x= -\\alpha e^{-\\alpha}$. Ainsi, $\\alpha=-W(-\\ln x)$, d'où\n", "\n", "$$\\varphi(x)=e^{-W(-\\ln x)}$$\n", "\n", "Certes, mais de quel $W$ parlons-nous ? S'agit-il de $W_1$ ou de $W_2$ ? Comme $e^{-e}\\le x \\le e^{1/e}$, on a $-e\\le \\ln x\\le \\frac 1 e$, donc $-\\frac 1 e \\le -\\ln x \\le e$. Si $-\\ln x \\ge 0$, c'est à dire $x\\le 1$, aucun doute nous parlons de $W_2$. Et sinon ? Dans ce cas, $-\\ln x$ appartient aux ensembles de définition de $W_1$ et $W_2$ et nous voilà bien embêtés. Nous admettrons que la \"bonne\" fonction est encore $W_2$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Proposition__ :\n", "\n", "- $\\varphi(e^{-e})=\\frac 1 e$\n", "- $\\varphi(e^{1/e})=e$\n", "- $\\varphi(\\sqrt 2)=2$\n", "\n", "__Démonstration__ : On a $W_2(-\\ln e^{-e})=W_2(e)=1$. Donc $\\varphi(e^{-e})=\\exp(W_2(-\\ln e^{-e}))=e^{-1}=\\frac 1 e$. Je vous laisse calculer les deux autres valeurs." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Proposition__ : la fonction $\\varphi$ est une bijection continue strictement croissante de $[e^{-e},e^{1/e}]$ sur $[\\frac 1 e, e]$. Elle est de classe $\\mathcal C^\\infty$ sur $[e^{-e},e^{1/e}[$ et dans cet intervalle on a\n", "\n", "$$\\varphi'(x)=\\frac{\\varphi(x)^2}{x(1-\\varphi(x)\\ln x)}$$\n", "\n", "__Démonstration__ : On a $\\varphi(x)=e^{-W_2(-\\ln x)}$. Il n'y a qu'à utiliser les résultats montrés pour $W_2$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Exemples__ :\n", "\n", "- $\\varphi'(e^{-e})=\\frac1 2 e^{e-2}$.\n", "- $\\varphi'(1)=1$.\n", "- Quand $x\\to e^{1/e}$, $\\varphi(x)\\to \\varphi(e^{1/e})=e$. On voit donc que $\\varphi'(x)$ tend vers $+\\infty$. Ainsi, $\\varphi$ n'est pas dérivable en $e^{1/e}$ et la courbe de $\\varphi$ a une tangente verticale en ce point. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "math.e ** (math.e - 2) / 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Si nous disposons d'un algorithme efficace du calcul de $W_i$, alors nous pourrons calculer efficacement $x^{[\\infty]}$ pour tout $x$. Nous allons tout d'abord implémenter la méthode de Newton, qui donne d'excellents résultats pour le calcul de $W$, sauf au voisinage de $-\\frac 1 e$. Puis nous verrons une méthode encore plus efficace, la méthode de Halley, une sorte de méthode de Newton de Newton." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.3 Calcul de $W$ par la méthode de Newton" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $f:I\\to \\mathbb R$ une fonction dérivable définie sur un intervalle de $\\mathbb R$. Il s'agit de trouver une racine de $f$. La méthode de Newton consiste à considérer la suite $(u_n)$ définie par $u_0\\in I$ et, pour tout entier $n$, $u_{n+1}=u_n-\\frac{f(u_n)}{f'(u_n)}$. Sous certaines conditions que je n'expliciterai pas ici (__voir le notebook sur les racines des équations !__), la suite $(u_n)$ converge vers une racine de $f$ sur l'intervalle $I$.\n", "\n", "Soit $x\\in \\mathbb R$. Notre problème ici est de trouver un réel $w$ tel que $we^w=x$. On considère donc la fonction $f$ définie par $f(w)=we^w-x$. On a $w-\\frac{f(w)}{f'(w)}=\\frac{w^2e^w+x}{(w+1)e^w}$, après quelques simplifications. \n", "\n", "La fonction `flamb` ci-dessous s'impose donc." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def flamb(w, x):\n", " e = math.exp(w)\n", " return (w ** 2 * e + x) / ((w + 1) * e)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "En itérant la fonction `flamb` nous devrions donc avoir des suites convergeant vers $W(x)$. Euh, oui, mais quel $W$ ? $W_1$ ou $W_2$ ? Sans démonstration, l'idée est de choisir le premier terme de la suite dans l'image de la fonction $W$ qui nous intéresse. Ainsi, si nous voulons calculer $W_1(x)$, nous choisissons $u_0=-2\\in W_1([-\\frac 1 e, 0[)$ et si nous voulons calculer $W_2(x)$, nous choisissons $u_0=2\\in W_1([-\\frac 1 e, +\\infty[)$. Ces choix sont arbitraires et risqués. On pourrait faire mieux (et on le fera plus loin) mais cela nous suffira pour l'instant.\n", "\n", "Voici donc la fonction `lambertw`. Elle prend en paramètres un réel $x$, un numéro $i$ de branche ($i=$ 1 ou 2), et un paramètre optionnel `dbg` qui permet d'afficher le nombre d'itérations effectuées. Elle renvoie $W_i(x)$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def lambertw(x, branche, dbg=False, tol=1e-14):\n", " cnt = 0\n", " if branche == 1: u = -2\n", " else: u = 2\n", " #while abs(u - v) > 1e-15:\n", " while abs(u * math.exp(u) - x) > tol:\n", " cnt += 1\n", " u = flamb(u, x)\n", " if dbg: print(cnt)\n", " return u" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "lambertw(-0.25, 1, True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "lambertw(-0.25, 2, True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "lambertw(-1 / math.e, 2, True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici la courbe de la fonction $W_2$. Et sans tricher cette fois-ci :-)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xs = subdivision(-1 / math.e, 3, 500)\n", "ys = [lambertw(x, 2) for x in xs]\n", "plt.plot(xs, ys, 'k')\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Et voici la courbe de la fonction $W_1$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xs = subdivision(-1 / math.e, -0.001, 500)\n", "ys = [lambertw(x, 1) for x in xs]\n", "plt.plot(xs, ys, 'k')\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Parfait. Nous pouvons calculer assez rapidement $W_i(x)$, $i=1,2$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.3' Petit aparté pour Romain Blanchard" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On désire résoudre l'équation d'inconnues $x,y>0$, $(E)\\ x^y=y^x$. On vérifie facilement que ceci équivaut à $\\frac{\\ln x}{x}=\\frac{\\ln y}{y}$. L'étude de la fonction $t\\mapsto \\frac{\\ln t} t$ permet de voir que les couples solutions vérifient $x,y>1$, l'un des deux réels $x,y$ étant inférieur à $e$ et l'autre supérieur à $e$.\n", "\n", "Posons $x=e^{-u}$ et $y=e^{-v}$. Le couple $(x,y)$ vérifie $(E)$ si et seulement si le couple $(u,v)$ vérifie $ue^u=ve^v$ (exercice). On commence à voir le rapport avec la fonction de Lambert ...\n", "\n", "Regardez la courbe de la fonction $f:x\\mapsto xe^x$, vue plus haut. Pour tout réel $a\\in]-\\frac 1 e,0[$ il existe exactement deux réels $u$ et $v$ tels que $f(u)=f(v)=a$. Ces deux réels sont bien sûr $W_1(a)$ et $W_2(a)$.\n", "\n", "Pour revenir à l'équation qui nous intéresse, on prend $x\\in]1,+\\infty[$ arbitraire. On a alors $u=-\\ln x< 0$ et $\\alpha=ue^u\\in]-\\frac 1 e, 0[$. On a $u=W_1(\\alpha)$ ou $u=W_2(\\alpha)$. En tout cas, le couple $(u',v')=(W_1(\\alpha), W_2(\\alpha))$, vérifie que $f(u')=f(v')$ et que l'un des deux réels $u'$ et $v'$ est $u$. Le couple $(x',y')=(e^{-u'},e^{-v'})$ est alors une solution de $(E)$ telle que l'un des deux réels $x'$ et $y'$ soit $x$.\n", "\n", "On en déduit la fonction `solxyyx` qui prend en paramètre un réel $x>1$ et renvoie le couple $(x,y)$ (ou $(y, x)$) tel que $x^y=y^x$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def solxyyx(x):\n", " u = - math.log(x)\n", " a = u * math.exp(u)\n", " u = lambertw(a, 1)\n", " v = lambertw(a, 2)\n", " x = math.exp(-u)\n", " y = math.exp(-v)\n", " return(x, y)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x, y = solxyyx(2)\n", "print('x = ', x)\n", "print('y = ', y)\n", "print('x^y = ', x ** y)\n", "print('y^x = ', y ** x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.4 $x^{[\\infty]}$, plus efficacement" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rappelons-nous l'inefficacité de notre fonction `inf_exp` au voisinage des points $e^{-e}$ et $e^{1/e}$. Écrivons donc une nouvelle fonction. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def inf_exp2(x, dbg=False, tol=1e-14):\n", " return math.exp(-lambertw(-math.log(x), 2, dbg, tol))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "inf_exp2(math.sqrt(2), True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = math.exp(-math.e)\n", "y = inf_exp2(x, True)\n", "print(y)\n", "print(1 / math.e)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "6 itérations seulement, au lieu de 233158 !" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = math.exp(1 / math.e)\n", "y = inf_exp2(x, True)\n", "print(y)\n", "print(math.e)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Et ici, 24 itérations. Remarquez tout de même que la précision obtenue n'est pas terrible. La méthode de Newton n'aime pas les points à tangente horizontale, ce qui est le cas ici pour la fonction $x\\mapsto x e^x$ en $-1/e$.\n", "\n", "Voici enfin la courbe tanta-tendue de $x\\mapsto x^{[\\infty]}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams['figure.figsize'] = (6, 10)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "xs = subdivision(math.exp(-math.e), math.exp(1 / math.e), 1000)\n", "ys = [inf_exp2(x) for x in xs]\n", "plt.plot(xs, ys, 'k')\n", "plt.xlabel('x')\n", "plt.ylabel('phi(x)')\n", "plt.axis('equal')\n", "plt.grid()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams['figure.figsize'] = (12, 6)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Rappelez-vous la remarque concernant le cafouillis au voisinage de $e^{-e}$. Nous n'avons plus ce problème. Notre nouvel alorithme se comporte bien sur tout l'ensemble de définition de $x\\mapsto x^{[\\infty]}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2.4 Calcul encore plus efficace de $W$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici un extrait de la documentation de la fonction `lambertw` du module `scipy` :\n", "\n", "...................................................................................................................\n", "\n", "All branches are supported by lambertw:\n", "\n", "`lambertw(z)` gives the principal solution (branch 0)\n", "`lambertw(z, k)` gives the solution on branch $k$\n", "The Lambert $W$ function has two partially real branches: the principal branch ($k = 0$) is real for real $z > -1/e$, and the $k = -1$ branch is real for $-1/e < z < 0$. All branches except $k = 0$ have a logarithmic singularity at $z = 0$.\n", "\n", "__Possible issues__\n", "\n", "The evaluation can become inaccurate very close to the branch point at $-1/e$. In some corner cases, `lambertw` might currently fail to converge, or can end up on the wrong branch.\n", "\n", "__Algorithm__\n", "\n", "Halley’s iteration is used to invert $w * exp(w)$, using a first-order asymptotic approximation ($O(log(w))$ or $O(w)$) as the initial estimate.\n", "\n", "The definition, implementation and choice of branches is based on [R207].\n", "\n", "..................................................................................................................." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Comment fonctionne la méthode de Halley ? Soit à résoudre l'équation $f(x)=0$. On considère la suite définie par son premier terme $x_0$ et la récurrence\n", "\n", "$$x_{n+1}=x_n-\\frac{2f(x_n)f'(x_n)}{2f'(x_n)^2-f(x_n)f''(x_n)}$$\n", "\n", "__Proposition__ : On suppose que $f$ est de classe $\\mathcal C^3$ ur un intervalle $I$, que $a\\in I$ vérifie $f(a)=0$ et $f'(a)\\ne 0$. Alors, pour tout $x_0$ suffisamment proche de $a$, la suite $(x_n)$ converge vers $a$. Il existe un réel $K>0$ tel que, pour tout $n$,\n", "\n", "$$|x_{n+1}-a|\\le K|x_n-a|^3$$\n", "\n", "La présence de la puissance 3 indique une vitesse foudroyante de convergence. La méthode de Halley est une __méthode d'ordre 3__. À titre de comparaison, la méthode de Newton est une méthode d'ordre 2." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Exercice__ : Prenez $f(w)=we^w-x$. Que devient la méthode de Halley pour $f$ ? Du coup, vous avez compris la fonction ci-dessous." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def flamb2(w, x):\n", " e = math.exp(w)\n", " wx = w * e - x\n", " return w - wx / (e * (w + 1) - (w + 2) * wx / (2 * (w + 1)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici notre nouvelle fonction de Lambert. J'ai effectué quelques modifications :\n", "\n", "\n", "- La considération à part du cas où $x=-\\frac 1 e$. Dans ce cas $f'(x)=0$ et la méthode de Halley échoue.\n", "- Le choix mystérieux du premier terme de la suite, $-1\\pm\\sqrt{2(ex+1)}$. Je n'entre pas ici dans les détails. Ce choix est dicté par un développement en série de la fonction de Lambert au voisinage de $x$.\n", "- Lorsque `dbg` vaut `True`, on renvoie le couple $(W_i(x), n)$ où $n$ est le nombre d'itérations effectuées." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def lambertw2(x, branch, dbg=False, tol=1e-14):\n", " cnt = 0\n", " p = math.sqrt(2 * (math.e * x + 1))\n", " if p == 0:\n", " if dbg: return (-1, 0)\n", " else: return -1\n", " if branch == 1: u = -1 - p\n", " else: u = -1 + p\n", " while abs(u * math.exp(u) - x) > tol:\n", " cnt += 1\n", " u = flamb2(u, x)\n", " if dbg: return (u, cnt)\n", " else: return u" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vérifions que tout va bien ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(lambertw(-0.25, 1, True))\n", "print(lambertw2(-0.25, 1, True))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(lambertw(-1 / math.e + 0.001, 2, True))\n", "print(lambertw2(-1 / math.e + 0.001, 2, True))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dans l'exemple ci-dessous, nous sommes obligés de diminuer la tolérance pour `lambertw`, sinon celle-ci échoue lamentablement. Pour `lambertw2`, en revanche, aucun problème." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(lambertw(100, 2, True, tol=1e-12))\n", "print(lambertw2(100, 2, True))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xs = subdivision(-1 / math.e, 5, 200)\n", "ys = [lambertw2(x, 2) for x in xs]\n", "plt.plot(xs, ys, 'k')\n", "plt.grid()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Et voici notre nouvelle exponentielle infinie." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def inf_exp3(x, dbg=False, tol=1e-14):\n", " if dbg:\n", " y, v = lambertw2(-math.log(x), 2, dbg, tol)\n", " return (math.exp(-y), v)\n", " else:\n", " return math.exp(-lambertw2(-math.log(x), 2, dbg, tol))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = math.exp(- math.e)\n", "print(inf_exp3(x, True))\n", "print(1 / math.e)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x = math.exp(1 / math.e)\n", "y = inf_exp3(x, True)\n", "print(y)\n", "print(math.e)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notre nouvelle fonction est-elle vraiment efficace ? Traçons, en fonction de $x$, le nombre d'itérations nécessaires pour calculer $x^{[\\infty]}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xs = subdivision(math.exp(-math.e), math.exp(1/math.e), 500)\n", "ys = []\n", "for x in xs:\n", " (_, n) = inf_exp3(x, True, tol=1e-15)\n", " ys.append(n)\n", "plt.plot(xs, ys, 'ko', markersize=2)\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Maximum 5 ! C'est __très__ efficace." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. $z^{z^{z^{z^\\ldots}}}$ dans $\\mathbb C$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "L'exposé qui suit reste à un niveau élémentaire. Je me contenterai de quelques explications et de dessins. Pas de preuves ...\n", "\n", "__Remarque__ : pour tout nombre complexe $z\\ne 0$ nous noterons $\\arg z$ l'unique argument de $z$ dans l'intervalle $]-\\pi,\\pi]$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.1 Logarithmes complexes" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import cmath" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Soit $z\\in\\mathbb C^*$. Quels sont les nombres complexes $w$ tels que $e^w=z$ ? Posons $z=re^{i\\theta}$ où $r>0$ est le module de $z$ et $-\\pi<\\theta\\le\\pi$ est son argument. Posons $w=x+iy$, où $x,y\\in\\mathbb R$. On a $e^w=z$ si et seulement si $e^x e^{iy}=re^{i\\theta}$, ou encore $x=\\ln r$ et $\\exists k\\in\\mathbb Z, y=\\theta+2k\\pi$. Ainsi, $w=\\ln |z| + i\\arg z + 2ik\\pi$.\n", "\n", "Posons, pour tout $k\\in\\mathbb Z$, $\\log_k(z)=\\ln |z| + i\\arg z + 2ik\\pi$. Nous venons de définir pour tout entier $k$ la $k$ième branche du _logarithme complexe_. Pour $k=0$, nous noterons plus simplement $\\log z = \\log_0 z$. La fonction $\\log$ est appelée la _branche principale_ du logarithme.\n", "\n", "Que dire de la fonction $\\log_k$ ? Eh bien, par définition même, $\\log_k : \\mathbb C^* \\to \\mathbb C$. Et pour tout $z\\ne 0$, nous avons $e^{\\log_k z}=z$.\n", "\n", "Soit maintenant $z=x+iy\\in\\mathbb C$. Que vaut $\\log_k e^z$ ? On a $|e^z| = e^x$ et $\\arg e^z=iy-2i\\ell\\pi$ où $\\ell\\in\\mathbb Z$ est l'unique entier tel que $-\\pi0$ est un réel \"petit\", on a $\\log(-1+i\\varepsilon)\\simeq i\\pi$ alors que $\\log(-1-i\\varepsilon)\\simeq -i\\pi$. les nombres complexes $-1+i\\varepsilon$ et $-1-i\\varepsilon$ sont très proches mais leurs logarithmes diffèrent d'environ $2i\\pi$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "z1 = cmath.log(-1+0.001j)\n", "z2 = cmath.log(-1-0.001j)\n", "print(z1)\n", "print(z2)\n", "print(z1 - z2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Exercice__ : Montrer que pour tous $z,z'\\in\\mathbb C$, on a $\\log(zz')\\equiv \\log z + \\log z'[2i\\pi]$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voici à titre d'illustration le graphe du logarithme. Mais c'est quoi le graphe d'une fonction de $\\mathbb C$ vers $\\mathbb C$ ? Il y a plusieurs façons de voir les choses. L'idée que nous allons retenir ici est la suivante. Prenons des courbes du plan et regardons leurs images par la fonction $\\log$. Nous allons obtenir de nouvelles courbes. En traçant ce nouvelles courbes nous voyons comment le logarithme _déforme_ le plan complexe.\n", "\n", "Quelles courbes choisir ? Voici quelques choix possibles :\n", "\n", "- Les droites parallèles aux axes de coordonnées. C'est le choix \"cartésien\", que nous allons prendre ici.\n", "- Les cercles centrés en $O$ et les demi-droites issues de $O$. C'est le choix \"polaire\".\n", "- D'autres familles de courbes ...\n", "\n", "Remarquez que nous prenons à chaque fois __DEUX__ familles de courbes qui se coupent à angles droits. Nous y reviendrons plus loin.\n", "\n", "Pour éviter le problème de discontinuité du logarithmes aux réels négatifs, nous ne traçons ici que pour $y>0$. Pour $y<0$ on obtiendrait une figure symétrique. Pourquoi ? Parce que :\n", "\n", "__Exercice__ : Montrer que pour tout $z\\in\\mathbb C^*$ on a $\\log \\overline z = \\overline{\\log z}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Images des droites x = cte\n", "xmin, xmax = -10, 10\n", "ymin, ymax = 0.1, 10\n", "for x in subdivision (xmin, xmax, 100):\n", " xs = []\n", " ys = []\n", " for y in subdivision(ymin, ymax, 200):\n", " if x != 0 or y != 0:\n", " w = cmath.log(x + 1j * y)\n", " xs.append(w.real)\n", " ys.append(w.imag)\n", " plt.plot(xs, ys, 'k')\n", "# Images des droites y = cte\n", "for y in subdivision (ymin, ymax, 20):\n", " xs = []\n", " ys = []\n", " for x in subdivision(xmin, xmax, 1000):\n", " if x != 0 or y != 0:\n", " w = cmath.log(x + 1j * y)\n", " xs.append(w.real)\n", " ys.append(w.imag)\n", " plt.plot(xs, ys, 'k')\n", "\n", "\n", "plt.grid()\n", "plt.axis('equal')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Que voyons-nous ? Les droites horizontales et verticales se sont transformées en des courbes dont nous pourrions trouver éventuellement l'équation. Je vous laisse faire cet exercice.\n", "\n", "__Remarque__ : Les courbes du dessin ont l'air de se couper à angles droits. Est-ce une coïncidence ? Non ! La fonction $\\log$ \"conserve les angles\". C'est ce que l'on appelle une __transformation conforme__." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.2 Puissances complexes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour $a,b\\in\\mathbb C$, $a\\ne 0$, posons $a^b=e^{b\\log a}$.\n", "\n", "__Exemples__ :\n", "\n", "- Si $b$ est entier, $a^b=e^{b\\log a}=(e^{\\log a})^b=\\ldots a^b$, au sens usuel de l'expression. Tout va bien !\n", "- $i^i=e^{i\\log i}=e^{-\\frac \\pi 2}\\in\\mathbb R$.\n", "\n", "Ci-dessous le graphe de $z\\mapsto 2^z$.\n", "\n", "__Exercice__ : Pourquoi ces choix pour `ymin` et `ymax` ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Images des droites x = cte\n", "xmin, xmax = -3, 3\n", "ymin, ymax = -math.pi / math.log(2), math.pi / math.log(2)\n", "for x in subdivision (xmin, xmax, 30):\n", " xs = []\n", " ys = []\n", " for y in subdivision(ymin, ymax, 200):\n", " w = 2 ** (x + 1j * y)\n", " xs.append(w.real)\n", " ys.append(w.imag)\n", " plt.plot(xs, ys, 'k')\n", "# Images des droites y = cte\n", "for y in subdivision (ymin, ymax, 30):\n", " xs = []\n", " ys = []\n", " for x in subdivision(xmin, xmax, 200):\n", " w = 2 ** (x + 1j * y)\n", " xs.append(w.real)\n", " ys.append(w.imag)\n", " plt.plot(xs, ys, 'k')\n", "\n", "plt.grid()\n", "plt.axis('equal')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Exercice__ : Commenter ce dessin. Demi-droites, cercles, angles droits." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Munis de puissances complexes quelconques, nous pouvons nous attaquer aux puissances infinies ..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.3 $z^{[\\infty]}$ dans $\\mathbb C$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous admettrons le résultat suivant :" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Proposition__ : Soit $z\\in\\mathbb C$. Le nombre complexe $z^{[\\infty]}=z^{z^{z^\\ldots}}$ est défini si et seulement si $z\\in V$ où\n", "\n", "$V = \\{e^{te^{-t}}, |t|<1$ ou $\\exists n\\ge 1$, $t^n=1\\}$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dessinons l'ensemble $V$. Plus précisément, dessinons $V_0 = \\{e^{te^{-t}}, |t|<1\\}$ et, pour tout $n\\ge 1$, les ensembles $V_n = \\{e^{te^{-t}}, t^n=1\\}$. Si nous posons $\\Phi(z)=e^{ze^{-z}}$, nous avons alors\n", "\n", "- $V_0=\\Phi(D)$ où $D=\\{z\\in\\mathbb C, |z|<1\\}$ est le disque ouvert de centre $O$ et de rayon 1\n", "- $V_n=\\Phi(\\mathcal U_n)$ où $\\mathcal U_n=\\{z\\in\\mathbb C, z^n=1\\}$ est l'ensemble des racines $n$ièmes de l'unité." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Commençons par $V_0$. Nous adoptons la méthode déjà vue plus haut, sauf que $D$ étant rond nous travaillons en coordonnées polaires. Nous traçons les images par $\\Phi$ des cercles de centre $O$ et des demi-droites issues de $O$.\n", "\n", "Dans le dessin ci-dessous, la ligne extérieure rouge correspond à $r=1$. Elle n'est donc pas incluse dans $V_0$. En fait, $V_0$ est la partie du plan située strictement à l'intérieur de cette ligne." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "for r in subdivision(0, 1, 20):\n", " xs = []\n", " ys = []\n", " for theta in subdivision (0, 2 * math.pi, 300):\n", " t = r * cmath.exp(1j * theta)\n", " w = cmath.exp(t * cmath.exp(-t))\n", " xs.append(w.real)\n", " ys.append(w.imag)\n", " if r < 1: col = 'k'\n", " else: col = 'r'\n", " plt.plot(xs, ys, col)\n", "for theta in subdivision(0, 2 * math.pi, 50):\n", " xs = []\n", " ys = []\n", " for r in subdivision (0, 1, 200):\n", " t = r * cmath.exp(1j * theta)\n", " w = cmath.exp(t * cmath.exp(-t))\n", " xs.append(w.real)\n", " ys.append(w.imag)\n", " plt.plot(xs, ys, 'k')\n", "plt.grid()\n", "plt.axis('equal')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ensemble compliqué, certes, mais abordable. Regardez en particulier les deux \"pincements\", ils se trouvent pile aux réels $e^{-e}$ et $e^{1/e}$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Et les ensembles $V_{n}$ ? Ils forment une partie dense de la ligne rouge. En effet, l'ensemble de toutes les racines $n$ièmes de l'unité, pour toutes les valeurs de $n$, est dense dans le cercle unité $\\mathcal U$. Il en résulte par continuité de $\\Phi$ que $\\cup_{n\\ge 1} V_n$ est dense dans la ligne extérieure, qui n'est autre que $\\Phi(\\mathcal U)$ !" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous avons donc bien dessiné l'ensemble $V$. Nous retrouvons l'intervalle $[e^{-e},e^{1/e}]$ pour les réels de cet ensemble. Et nous voyons que $i\\in V$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.4 $W:\\mathbb C\\to\\mathbb C$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Réécrivons la fonction de Lambert, ainsi que les puissances infinies pour que tout cela fonctionne avec un paramètre complexe.\n", "\n", "__Remarque__ : La fonction de Lambert complexe possède en fait une __infinité__ de branches. Nous nous cantonnons dans ce qui suit aux deux branches déjà évoquées pour les réels." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def flambc(w, x):\n", " e = cmath.exp(w)\n", " wx = w * e - x\n", " return w - wx / (e * (w + 1) - (w + 2) * wx / (2 * (w + 1)))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def lambertwc(x, branch, tol=1e-14):\n", " p = cmath.sqrt(2 * (math.e * x + 1))\n", " if p == 0: return -1\n", " if branch == 1: u = -1 - p\n", " else: u = -1 + p\n", " while abs(u * cmath.exp(u) - x) > tol:\n", " u = flambc(u, x)\n", " return u" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def inf_expc(x, tol=1e-14):\n", " return cmath.exp(-lambertwc(-cmath.log(x), 2, tol))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "inf_expc(1j)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "inf_expc(math.sqrt(2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Et maintenant, traçons le graphe de $z\\mapsto z^{[\\infty]}$ ! Pour être honnête, l'ensemble de départ de la fonction étant très compliqué, traçons le graphe de $\\Phi(z)^{[\\infty]}$. L'ensemble de départ de cette fonction étant un disque, le choix de coordonnées polaires s'impose." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "for r in subdivision(0, 1, 20):\n", " xs = []\n", " ys = []\n", " for theta in subdivision (0, 2 * math.pi, 300):\n", " t = r * cmath.exp(1j * theta)\n", " w = inf_expc(cmath.exp(t * cmath.exp(-t)))\n", " xs.append(w.real)\n", " ys.append(w.imag)\n", " if r < 1:\n", " plt.plot(xs, ys, 'k')\n", " else:\n", " plt.plot(xs, ys, 'r')\n", "for theta in subdivision(0, 2 * math.pi, 50):\n", " xs = []\n", " ys = []\n", " for r in subdivision (0, 1, 200):\n", " t = r * cmath.exp(1j * theta)\n", " w = inf_expc(cmath.exp(t * cmath.exp(-t)))\n", " xs.append(w.real)\n", " ys.append(w.imag)\n", " plt.plot(xs, ys, 'k')\n", "plt.grid()\n", "plt.axis('equal')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Exercice__ : Comprenez le dessin ci-dessus." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3.5 Que se passe-t-il __en dehors__ de l'ensemble de définition ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Que se passe-t-il __en dehors__ de l'ensemble de définition de $z\\mapsto z^{[\\infty]}$ ? Eh bien, la suite $(z^{[n]})$ diverge, me direz-vous. Oui, mais il y a plusieurs façons de diverger. La suite peut tendre (en module) vers l'infini, ou alors les termes d'indices pairs converger vers une certaine limite, et ceux d'indice impair vers une autre. Ou alors il faut peut-être regarder modulo 3 ou 4. Ou alors la suite peut faire n'importe quoi.\n", "\n", "Bref, tentons une expérience. Pour $z\\in\\mathbb C$, calculons $z^{[0]}, z^{[1]}, z^{[2]}$, etc, jusqu'à ce que $|z^{[n]}|$ devienne \"très grand\" (disons plus grand que 200), ou que $n$ dépasse une certaine valeur `niter`. Et renvoyons le nombre d'itérations effectuées. La fonction `iterer` fait le travail. Elle renvoie également 0 si la suite à l'air de converger." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def iterer(z, niter):\n", " k = 0\n", " w = 1\n", " w1 = z ** w\n", " while k < niter and abs(w1) < 200:\n", " w, w1 = w1, z ** w1\n", " k = k + 1\n", " if abs(w - w1) < 1e-5: k = 0\n", " return k" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Par exemple, pour $z=2$ la suite tend très vite vers l'infini." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "iterer(2, 256)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour $z=i$, la suite converge." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "iterer(1j, 256)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pour $z=3+2i$, on a divergence, mais la suite diverge \"moins vite\" que pour \"z=2\"." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "iterer(3+2j, 256)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Maintenant, itérons pour tous les $z$ appartenant à un certain rectangle et mettons les valeurs retournées par `iterer` dans une matrice. Puis affichons. Nous obtenons une carte qui nous donne la \"vitesse de divergence\" en fonction de $z$. Je ne commenterai pas. Contentez-vous de regarder.\n", "\n", "Soyez patients, l'affichage prend un certain temps (environ 7 secondes sur ma machine)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams['figure.figsize'] = (12, 12)\n", "a = 3\n", "xmin, xmax = -a, a\n", "ymin, ymax = -a, a\n", "nx, ny = 400, 400\n", "niter = 64\n", "M = (ny + 1) * [None]\n", "for i in range(ny + 1):\n", " M[i] = (nx + 1) * [0]\n", "for i in range(nx + 1):\n", " for j in range(ny + 1):\n", " x = xmin + i * (xmax - xmin) / nx\n", " y = ymin + j * (ymax - ymin) / ny\n", " z = x + 1j * y\n", " M[j][i] = 4 * iterer(z, niter)\n", "plt.imshow(M, origin='lower', cmap='hot', interpolation='bicubic', extent=(xmin, xmax, ymin, ymax))\n", "plt.grid()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Nous voyons au centre l'ensemble déjà vu plus haut (une espèce de néphroïde diront certains), là où la suite converge. Et tout autour, des tas de choses intéressantes. Une structure très compliquée, et certainement digne d'être étudiée. Ce que je ne ferai pas ici.\n", "\n", "Quelques suggestions :\n", "\n", "- Augmentez les valeurs de `nx` et `ny`, à 600 par exemple. Attention, multiplier ces valeurs par 2 multiplie le temps de calcul par 4.\n", "- Faites un zoom autour des points $-1$, $e^{-e}$ et $e^{\\frac 1 e}$ en ajustant les valeurs de `xmin`, `xmax`, `ymin` et `ymax`.\n", "- Cherchez d'autres points intéressants." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Petite bibliographie" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Outre les sources standard d'information, les deux articles ci-dessous m'ont été d'une aide précieuse. Ils se trouvent facilement sur Internet.\n", "\n", "- Corless, Gonnet, Hare, Jeffrey, Knuth - _On the Lambert $W$ Function_, Advances in Computational Mathematics, __5__, 329-359.\n", "- Galidakis - _On an Application of Lambert's $W$ Function to Infinite Exponentials_, Complex Variables, Vol 49, N$^o$ 11, 15 September 2004, 759-780.\n", "- [Documentation de scipy : fonction lambertw](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.lambertw.html#scipy.special.lambertw)\n", "- [Documentation de mpmath : fonction lambertw](http://mpmath.org/doc/current/functions/powers.html)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.4" } }, "nbformat": 4, "nbformat_minor": 2 }