{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Un joli résultat sur les polynômes\n",
    "\n",
    "Marc Lorenzi\n",
    "\n",
    "26 mars 2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "import numpy.polynomial.polynomial as poly\n",
    "import cmath, math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (10, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 Qu'allons-nous faire ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous allons dans ce petit notebook prouver le théorème suivant :\n",
    "\n",
    "__Théorème (Gauss-Lucas)__ : Soit $P\\in\\mathbb C[X]$ un polynôme non constant. Les racines de $P'$ sont dans l'enveloppe convexe de l'ensemble des racines de $P$.\n",
    "\n",
    "Nous illustrerons ce théorème avec quelques expériences pythonesques sur des polynômes \"aléatoires\".\n",
    "\n",
    "__Convention__ : Dans toute la suite nous identifierons les nombres complexes et les __points__ du plan, ainsi que les nombres complexes et les __vecteurs__ du plan."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Droites, demi-plans"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $a\\in\\mathbb C$, vu comme un point du plan. Soit $b\\in\\mathbb C^*$, vu comme un vecteur non nul. Les points de la droite $\\mathcal D$ passant par $a$ et dirigée par $b$ sont les $z\\in\\mathbb C$ de la forme $z=a+bt$ où $t\\in\\mathbb R$. Dit autrement, $z\\in\\mathcal D$ si et seulement si\n",
    "\n",
    "$${\\rm Im\\ } \\frac{z-a}{b}=0$$\n",
    "\n",
    "Orientons la droite $\\mathcal D$ par le vecteur $b$. $\\mathcal D$ partage le plan en deux demi-plans, l'un, $H_g$, à gauche de $\\mathcal D$, l'autre, $H_d$, à droite. À droite ? À gauche ? Oui, grâce à l'orientation que nous avons imposée à $\\mathcal D$. En nous plaçant sur la droite $\\mathcal D$ et en regardant dans le sens de $b$, $H_g$ est à notre gauche et $H_d$ à notre droite. Ces deux demi-plans sont caractérisés par des inégalités :\n",
    "\n",
    "$$(H_g) \\quad {\\rm Im\\ } \\frac{z-a}{b}>0$$\n",
    "\n",
    "et\n",
    "\n",
    "$$(H_d) \\quad {\\rm Im\\ } \\frac{z-a}{b}<0$$\n",
    "\n",
    "On définit également les demi-plans __fermés__ $\\overline H_g=H_g\\cup \\mathcal D$ et $\\overline H_d=H_d\\cup \\mathcal D$. Il suffit de remplacer les inégalités strictes ci-dessus par des inégalités larges."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exemple__ : Prenons $a=0$ et $b=1$. La droite $\\mathcal D$ est alors l'axe $Ox$, orienté dans le sens usuel. $H_g$ est donné par ${\\rm Im\\ }z>0$, c'est le demi-plan supérieur de $\\mathbb C$. Et $H_d$ est le demi-plan inférieur.\n",
    "\n",
    "__Exercice__ : Prenez $a=0$ et $b=i$. Posez $z=x+iy$, où $x,y\\in\\mathbb R$. À quelle condition sur $x$ et $y$ a-t-on $z\\in H_g$ ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercice__ :\n",
    "1. Montrer que pour tous $z,z'\\in H_g$, le segment $[z,z']$ est inclus dans $H_g$. Ce qui montre que $H_g$ est __convexe__.\n",
    "2. Montrer que pour tout $z\\in H_g$ et tout $z'\\in H_d$, le segment $[z,z']$ contient un point de  $\\mathcal D$.\n",
    "\n",
    "Avec un peu de dextérité mathématique on pourrait en fait montrer le résultat numéro 2 ci-dessus en remplaçant \"segment $[z,z']$\" par \"courbe d'origine $z$ et d'extrémité $z'$\". Pour les adeptes de topologie, cela prouve que $\\mathbb C\\setminus\\mathcal D$ possède deux composantes connexes, qui sont $H_g$ et $H_d$.\n",
    "\n",
    "__Remarque__ : On peut aussi faire un dessin et se dire que c'est évident :-)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Une décomposition en éléments simples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soi $P$ un polynôme de degré $n\\ge 1$. Soient $a_1, a_2,\\ldots,a_n$ les racines (éventuellement égales) de $P$. On a\n",
    "\n",
    "$$P=c\\prod_{k=1}^n(X-a_k)$$\n",
    "\n",
    "où $c\\in\\mathbb C^*$ est le coefficient dominant de $P$. Dérivons $P$. On obtient une somme de $n$ termes, chaque terme étant obtenu en dérivant l'un des facteurs de $P$.\n",
    "\n",
    "$$P'=c\\sum_{j=1}^n\\prod_{k\\ne j}(X-a_k)=\\sum_{j=1}^n\\frac{P}{X-a_j}$$\n",
    "\n",
    "Ainsi,\n",
    "\n",
    "$$\\frac{P'}{P}=\\sum_{j=1}^n\\frac 1 {X-a_j}$$\n",
    "\n",
    "__Avant de voir le rapport avec notre théorème, faisons quelques expériences en Python.__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Racines d'un polynôme et de sa dérivée"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 `numpy` et les polynômes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le module `numpy.polynomial.polynomial` contient un certain nombre de fonctions permettant de manipuler des polynômes. Un polynôme $P=\\sum_{k=0}^n a_k X^k$ y est représenté par la liste $[a_0,\\ldots,a_n]$ (ou, de façon équivalente, un tableau `numpy`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Polynômes \"aléatoires\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fabriquons une fonction `random_poly`. Elle prend en paramètre un entier $n$. Elle calcule $n$ nombres complexes \"aléatoires\" $z_1,\\ldots, z_n$ de parties réelle et imaginaire entre $-1$ et $1$ et renvoie le polynôme\n",
    "\n",
    "$$P=\\prod_{k=1}^n (X-z_k)$$\n",
    "\n",
    "La fonction `poly.polyfromroots` fait le travail à notre place. Elle prend la liste des racines en paramètre et renvoie la liste des coefficients du polynôme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_poly(n):\n",
    "    s = [random.uniform(-1, 1) + 1j * random.uniform(-1, 1) for k in range(n)]\n",
    "    # print(s)\n",
    "    P = poly.polyfromroots(s)\n",
    "    return P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un petit exemple ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = random_poly(10)\n",
    "print(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remarquez le coefficient dominant égal à 1. C'est normal ! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Représentation graphique d'un polynôme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comment représenter graphiquement une fonction polynôme $P:\\mathbb C\\to\\mathbb C$ ? Comme nous nous intéressons ici aux racines des polynômes, nous allons dessiner $|P|$. La racines de $P$ donnent lieu au minimum, nul, de son module.\n",
    "\n",
    "Une première possibilité serait une représentation sous forme de surface, mais les résultats obtenus pour les polynômes sont assez décevants. Orientons-nous vers un système de couleurs. Chaque réel positif sera associé à une couleur : il nous suffira donc avec cette méthode de remplir une matrice dont les coefficients sont les valeurs du module de $P$. Plus précisément, le coefficient ligne $i$, colonne $j$ de la matrice est $|P(z)|$ où $z$ est le nombre complexe qui \"correspond\" au couple $(i, j)$. Nos polynômes ayant leurs racines dans le disque de centre $O$ et de rayon 1, nous tracerons $P$ pour des nombres complexes de parties réelle et imaginaire entre $-1.2$ et $1.2$. Il n'est alors pas bien difficile d'établir une correspondance entre le couple $(i, j)$ et $z$.\n",
    "\n",
    "La fonction `plot_poly` prend en paramètres\n",
    "\n",
    "- Un polynôme $P$ donné par la liste de ses coefficients.\n",
    "- Un entier $N$ qui est le nombre de lignes (et de colonnes) de la matrice à calculer.\n",
    "- Un réel positif `m`.\n",
    "\n",
    "La fonction renvoie les valeurs de $z \\mapsto |P(z)|^{1/n}\\bmod m$, où $n$ est le degré de $P$, dans une matrice. Pourquoi ce paramètre $m$ ? Plus $m$ est petit, plus les couleurs varieront vite en fonction de $z$. Une petite valeur de $m$ permet donc d'augmenter le contraste au voisinage des points où le polynôme varie peu. Pour ce que nous voulons observer, $m=0.1$ est un bon compromis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_poly(P, N=200, m=0.1):\n",
    "    e = 1 / (len(P) - 1)\n",
    "    z = [[0 for i in range(N)] for j in range(N)]\n",
    "    for i in range(N):\n",
    "        for j in range(N):\n",
    "            x = -1.2 + j * 2.4 / N\n",
    "            y = 1.2 - i * 2.4 / N\n",
    "            w = abs(poly.polyval(x + 1j * y, P)) ** e\n",
    "            if w % (2 * m) <= m: z[i][j] = (w % m) / m\n",
    "            else: z[i][j] = 1 - (w % m) / m\n",
    "    plt.imshow(z, cmap='viridis', interpolation='bicubic', extent=[-1.2, 1.2, -1.2, 1.2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici un exemple sur un polynôme dont on connaît bien les racines : $P=X(X^3-1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "j = cmath.exp(2 * 1j * math.pi / 3)\n",
    "P = poly.polyfromroots([0, 1, j, j ** 2])\n",
    "plot_poly(P)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et un test sur un polynôme aléatoire."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = random_poly(15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_poly(P, m=0.1)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercice__ : Trouvez sur le dessin les racines de $P$.\n",
    "\n",
    "__Exercice__ : Changez ci-dessus la valeur de $m$ pour voir ce qui se passe."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Racines de $P$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculons idiotement les racines de $P$. La fonction `poly.polyroots` est là pour ça.\n",
    "\n",
    "__Exercice__ : Pourquoi idiotement ? _Parce que nous avons calculé $P$ à partir de ses racines !_ Décommentez la ligne `print(s)` dans `random_poly` pour constater que `polyroots` fait bien son travail. Puis re-commentez-la."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rs = poly.polyroots(P)\n",
    "print(rs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Et dessinons les racines par dessus le tracé du polynôme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = random_poly(20)\n",
    "plot_poly(P, m=0.1)\n",
    "rs = poly.polyroots(P)\n",
    "xs = [complex(z).real for z in rs]\n",
    "ys = [complex(z).imag for z in rs]\n",
    "plt.plot(xs, ys, 'or', markersize=10)\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercice__ : Le tracé des cercles sur les racines est presque superflu. Commentez la ligne correspondante, vous verrez les racines tout aussi bien !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.5 $P'$ et ses racines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dérivons maintenant $P$ grâce à la fonction `numpy.polyder`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P1 = poly.polyder(P)\n",
    "print(P1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remarquez le coefficient dominant égal au degré de $P$ :-).\n",
    "\n",
    "Calculons illico les racines de $P'$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rs1 = poly.polyroots(P1)\n",
    "print(rs1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enfin, affichons le tout ! En rouge, les racines de $P'$. En noir, les racines de $P$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs1 = [complex(z).real for z in rs1]\n",
    "ys1 = [complex(z).imag for z in rs1]\n",
    "plt.plot(xs1, ys1, 'or')\n",
    "plt.plot(xs, ys, 'ok')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.6 Combinons le tout"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Regroupons le tout dans une fonction `test_theoreme`. Les racines de $P$ ne sont pas affichées. Les racines de $P'$ sont marquées par des cercles rouges et blancs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_theoreme(P, **opt):\n",
    "    plot_poly(P, **opt)\n",
    "    #rs = poly.polyroots(P)\n",
    "    #xs = [complex(z).real for z in rs]\n",
    "    #ys = [complex(z).imag for z in rs]\n",
    "    #plt.plot(xs, ys, 'or', markersize=10)\n",
    "    P1 = poly.polyder(P)\n",
    "    rs1 = poly.polyroots(P1)\n",
    "    xs1 = [complex(z).real for z in rs1]\n",
    "    ys1 = [complex(z).imag for z in rs1]\n",
    "    plt.plot(xs1, ys1, 'ow', markersize=10)\n",
    "    plt.plot(xs1, ys1, 'or', markersize=5)\n",
    "    plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vous pouvez maintenant tester. N'abusez pas sur le degré du polynôme : pour cause d'erreurs d'arrondi et de dépassement de capacité des flottants, `numpy` renvoie des résultats erronés si le degré du polynôme est trop grand (essayez 200, juste pour rire). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = random_poly(10)\n",
    "#P = poly.polyfromroots([0, 1, -0.5 + math.sqrt(3)/2 * 1j, -0.5 - math.sqrt(3)/2 * 1j, math.sqrt(2) / 2 * (1 + 1j)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_theoreme(P, m=0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il est intéressant de se poser la question : Et alors ? Que voit-on ?\n",
    "\n",
    "En fait on voit bien des choses mais nous allons nous fixer sur l'une d'entre-elles ...\n",
    "\n",
    "Si vous tracez une droite, celle qui vous fait plaisir, et que tous les points noirs sont d'un côté de la droite, alors tous les points rouges sont aussi du même côté. En d'autres termes, tout demi-plan contenant les racines de $P$ contient aussi les racines de $P'$.\n",
    "\n",
    "Si nous notons $\\mathcal C(P)$ l'intersection de __tous__ les demi-plans contenant les racines de $P$, alors les racines de $P'$ appartiennent à $\\mathcal C(P)$.\n",
    "\n",
    "On peut montrer que $\\mathcal C(P)$ est une partie convexe du plan (ça c'est facile), et que c'est la plus petite partie convexe du plan contenant les racines de $P$ (c'est un peu moins facile). C'est ce que l'on appelle __l'enveloppe convexe__ de l'ensemble des racines. Encore mieux, $\\mathcal C(P)$ est un polygone convexe dont les sommets sont des racines de $P$.\n",
    "\n",
    "__NB__ : Si vous voulez plus de détails sur le sujet des enveloppes convexes, consultez le notebook qui lui est consacré. Un algorithme permettant le calcul efficace de l'enveloppe convexe s'appelle __l'algorithme de Graham__."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercice__ : Réévaluez la cellule ci-dessus avec un polynôme unitaire de degré 2 ayant deux racines distinctes. Où est l'unique racine de $P'$ ? Prouvez-le, c'est facile.\n",
    "\n",
    "__Exercice__ : Réévaluez la cellule ci-dessus avec un polynôme unitaire de degré 3 ayant trois racines distinctes et tel que $P'$ n'ait qu'une racine, double. Où est l'unique racine de $P'$ ? Prouvez-le aussi. Indication : $P'=3(X-a)^2$, cela ne devrait pas être difficile de trouver $P$ et ses racines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Fin de la démonstration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $P$ un polynôme de degré $n\\ge 1$ dont les racines (éventuellement égales) sont $a_1,a_2,\\ldots,a_n$.\n",
    "\n",
    "Soient $a\\in\\mathbb C$ et $b\\in\\mathbb C^*$. Soit $\\overline H$ le demi-plan fermé défini par\n",
    "\n",
    "$$(\\overline H) \\quad {\\rm Im\\ } \\frac{z-a}{b}\\ge 0$$\n",
    "\n",
    " Soit $H'$ son complémentaire. C'est le demi-plan ouvert défini par\n",
    "\n",
    "$$(H') \\quad {\\rm Im\\ } \\frac{z-a}{b}<0$$\n",
    "\n",
    "Supposons que $\\overline H$ contient toutes les racines du polynôme $P$. On a donc pour tout $k$ entre 1 et $n$\n",
    "\n",
    "$${\\rm Im\\ } \\frac{a_k-a}{b}\\ge 0$$\n",
    "\n",
    "Soit $z\\in H'$. On a donc ${\\rm Im\\ } \\frac{z-a}{b}<0$. Pour tout $k$ entre 1 et $n$, on a\n",
    "\n",
    "$${\\rm Im\\ } \\frac{z-a_k}{b}={\\rm Im\\ } (\\frac{z-a}{b}-\\frac{a_k-a}{b})={\\rm Im\\ } \\frac{z-a}{b}-{\\rm Im\\ }\\frac{a_k-a}{b}<0$$\n",
    "\n",
    "Le passage à l'inverse change le signe d'un nombre complexe. On a donc\n",
    "\n",
    "$${\\rm Im\\ } \\frac b{z-a_k}> 0$$\n",
    "\n",
    "Rappelons-nous maintenant notre décomposition en éléments simples :\n",
    "\n",
    "$$\\frac{P'}{P}=\\sum_{k=1}^n\\frac 1 {X-a_k}$$\n",
    "\n",
    "On en déduit\n",
    "\n",
    "$$b\\frac{P'(z)}{P(z)}=\\sum_{k=1}^n\\frac {b} {z-a_k}$$\n",
    "\n",
    "D'où, en prenant la partie imaginaire :\n",
    "\n",
    "$${\\rm Im\\ } \\left(b\\frac{P'(z)}{P(z)}\\right)=\\sum_{k=1}^n{\\rm Im\\ } \\frac {b} {z-a_k}> 0$$\n",
    "\n",
    "La partie imaginaire de $b\\frac{P'(z)}{P(z)}$ est non nulle, donc $P'(z)\\ne 0$. Contraposons :\n",
    "\n",
    "Si $P'(z)=0$, alors $z\\not\\in H'$, donc $z\\in \\overline H$. Et ceci est vrai pour __tout__ demi-plan (fermé) contenant les racines de $P$.\n",
    "\n",
    "__Notre théorème est démontré.__\n",
    "\n"
   ]
  },
  {
   "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
}
