{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpolation de Lagrange"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Marc Lorenzi\n",
    "\n",
    "Février 2019"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "init_printing()\n",
    "x = Symbol('x')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (10, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. La théorie"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dans tout ce notebook, \n",
    "\n",
    "- $n$ désigne un entier naturel. \n",
    "- $I$ est un intervalle de $\\mathbb R$. \n",
    "- $x_0<x_1<\\ldots<x_n$ sont $n+1$ points distincts de $I$.\n",
    "- $f$ est une fonction de $I$ vers $\\mathbb R$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.1 Le problème de l'interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Problème__ : trouver un polynôme $P$ de degré inférieur ou égal à $n$ tel que\n",
    "\n",
    "$$\\forall i\\in[0,n], P(x_i)=f(x_i)$$\n",
    "\n",
    "Nous allons voir que ce problème admet une unique solution. Le polynôme correspondant est appelé le polynôme d'interpolation de $f$ aux points $x_i$.\n",
    "\n",
    "Nous allons dans ce notebook :\n",
    "\n",
    "1. Montrer l'existence et l'unicité d'un tel polynôme.\n",
    "2. Calculer efficacement la valeur du polynôme d'interpolation en un point $x$.\n",
    "3. Majorer $|f(x)-P(x)|$ pour $x\\in I$.\n",
    "4. Voir comment rendre ce majorant aussi petit que possible en choisissant astucieusement les $x_i$.\n",
    "5. Apprendre à marcher sur les mains.\n",
    "6. Et bien d'autres choses ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Les polynômes élémentaires de Lagrange"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour $k=0,\\ldots,n$ le $k$ième polynôme de Lagrange est le polynôme \n",
    "\n",
    "$$L_k=\\frac{\\prod_{j\\ne k}(X-x_j)}{\\prod_{j\\ne k}(x_k-x_j)}$$\n",
    "\n",
    "C'est l'unique polynôme de degré inférieur ou égal à $n$ qui s'annule en tous les $x_j$, sauf en $x_k$ où il prend la valeur 1.\n",
    "\n",
    "__Proposition__ : $\\mathcal B=(L_0,L_1,\\ldots,L_n)$ est une base de l'espace $\\mathbb R_n[X]$ des polynômes de degré inférieur ou égal à $n$.\n",
    "\n",
    "__Démonstration__ : Comme ${\\rm Card\\ } \\mathcal B=\\dim\\mathbb R_n[X]$ il suffit de montrer que $\\mathcal B$ est, par exemple, libre. Soient $\\lambda_k, k=0,1,\\ldots,n$ $n+1$ réels. Supposons que\n",
    "\n",
    "$$\\sum_{k=0}^n\\lambda_kL_k=0$$\n",
    "\n",
    "On a pour tous $i,k$ $L_k(x_i)=\\delta_{ki}$ où le symbole de Kronecker $\\delta_{ki}$ vaut 1 si $k=i$ et $0$ sinon. Évaluons l'égalité ci-dessus en $x_i$. Il vient\n",
    "\n",
    "$$0=\\sum_{k=0}^n\\lambda_kL_k(x_i)=\\sum_{k=0}^n\\lambda_k\\delta_{ki}=\\lambda_i$$\n",
    "\n",
    "ceci pour tout $i$. D'où la liberté. Tout polynôme $P$ de degré inférieur ou égal à $n$ est donc combinaison linéaire des $L_k$, d'où leur nom \"élémentaires\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `lagrange` prend en paramètre un entier $k$, un réel $x$ et une liste $xs=[x_0,\\ldots,x_n]$ de réels distincts. Elle renvoie $L_k(x)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lagrange(k, x, xs):\n",
    "    p = 1\n",
    "    n = len(xs) - 1\n",
    "    for j in range(n + 1):\n",
    "        if j != k:\n",
    "            p *= (x - xs[j]) / (xs[k] - xs[j])\n",
    "    return p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ayant défini au début du notebook $x$ comme un symbole, nous pouvons donc obtenir une experssion explicite des $L_k$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expand(lagrange(1, x, [-2, -1, 0, 1, 2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "expand(lagrange(2, x, [-2, -1, 0, 1, 2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Subdivisions régulières"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soient $a, b\\in\\mathbb R$, $a<b$. Une subdivision du segment $[a,b]$ est une suite $x_0=a<x_1<\\ldots<x_n=b$. Pour tout entier $n\\ge 1$, la subdivision régulière à $n+1$ points de $[a, b]$ est définie par \n",
    "\n",
    "$$x_k=a+k\\frac{b-a}n, k=0,\\ldots,n$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def subdi(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 = subdi(-1, 1, S(10))\n",
    "xs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.4 Graphes des polynômes élémentaires"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous pouvons maintenant effectuer quelques petits affichages. Commençons par le graphe des polynômes élémentaires."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = subdi(-1, 1, S(5))\n",
    "ps = [expand(lagrange(k, x, xs)) for k in range(len(xs))]\n",
    "ps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `tracer_poly` trace la courbe de $p$ pour $x$ entre `xmin` et `xmax`. Les paramètres optionnels de la fonction sont passés à `plt.plot`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tracer_poly(p, xmin, xmax, *opt1, **opt2):\n",
    "    xs = subdi(xmin, xmax, 300)\n",
    "    ys = [p.subs(x, t) for t in xs]\n",
    "    plt.plot(xs, ys, *opt1, **opt2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici les polynômes élémentaires de Lagrange associés à une subdivision régulière."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = subdi(-1, 1, 10)\n",
    "ps = [expand(lagrange(k, x, xs)) for k in range(len(xs))]\n",
    "for p in ps:\n",
    "    tracer_poly(p, -1., 1., 'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercice__ : Parmi les courbes ci-dessus, trouver la courbe de $L_3$. Facile, $L_3$ est le seul à ne pas s'annuler en $x_3$. Une fois que vous l'avez trouvée, surlignez-la sur l'écran de votre ordinateur avec un feutre rouge :-)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.5 Le polynôme d'interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Proposition__ : Il existe un unique polynôme $P$ de degré inférieur ou égal à $n$ tel que $P(x_k)=f(x_k)$ pour $k=0,\\ldots,n$. Le polynôme $P$ est donné par \n",
    "\n",
    "$$P=\\sum_{k=0}^n f(x_k)L_k$$\n",
    "\n",
    "où les $L_k$ sont les polynômes élémentaires de Lagrange.\n",
    "\n",
    "__Démonstration__ : Soit $P$ le polynôme défini ci-dessus. Tout d'abord, $P$ est de degré inférieur ou égal à $n$. De plus, pour $i=0,1,\\ldots,n$,\n",
    "\n",
    "$$P(x_i)=\\sum_{k=0}^n f(x_k)L_k(x_i)=\\sum_{k=0}^n f(x_k)\\delta_{ki}=f(x_i)$$\n",
    "\n",
    "Ce polynôme est donc solution du problème. Par ailleurs, supposons que $P$ et $Q$ sont solutions. Le polynôme $P-Q$ est de degré inférieur ou égal à $n$ et s'annule en les $n+1$ points $x_0, x_1,\\ldots, x_n$. C'est donc le polynôme nul. D'où l'unicité."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `interpoler` prend en paramètres une fonction $f$ et une liste $xs$ de points. Elle renvoie le polynôme d'interpolation de $f$ aux points de la liste $xs$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def interpoler(f, xs):\n",
    "    ps = [lagrange(k, x, xs) for k in range(len(xs))]\n",
    "    p = 0\n",
    "    n = len(xs)\n",
    "    for k in range(n):\n",
    "        p = p + f(xs[k]) * ps[k]\n",
    "    return p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Testons sur un petit exemple. Prenons $f:x\\mapsto \\frac 1 {x^2+1}$ et interpolons sur une subdivision de $[-2,2]$. Nous allons fréquemment utiliser cette fonction dans toute la suite, alors appelons-la _l'exemple_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exemple(x):\n",
    "    return 1 / (x ** 2 + 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = expand(interpoler(exemple, subdi(-2, 2, S(5))))\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = subdi(-2, 2, 300)\n",
    "ys = [exemple(t) for t in xs]\n",
    "plt.plot(xs, ys, 'r')\n",
    "tracer_poly(p, -2, 2, 'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Que se passe-t-il si on augmente le nombre de points de la subdivision ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = expand(interpoler(exemple, subdi(-2, 2, 10)))\n",
    "xs = subdi(-2, 2, 300)\n",
    "ys = [exemple(t) for t in xs]\n",
    "plt.plot(xs, ys, 'r')\n",
    "tracer_poly(p, -2, 2, 'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sur cet exemple, augmenter le nombre de points d'interpolation conduit à deux phénomènes : si l'on excepte un voisinage des bords de l'intervalle, on augmente la qualité de l'interpolation. En revanche, au bord de l'intervalle, on assiste à une dégradation de l'interpolation. C'est le phénomène de Runge ... "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un deuxième exemple avant de poursuivre : la fonction sinus sur l'intervalle $[0,2\\pi]$. Nous l'appellerons ... $\\sin$ :-)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = expand(interpoler(math.sin, subdi(0, 2 * math.pi, 5)))\n",
    "xs = subdi(0, 2 * math.pi, 300)\n",
    "ys = [math.sin(t) for t in xs]\n",
    "plt.plot(xs, ys, 'r')\n",
    "tracer_poly(p, 0, 2 * math.pi, 'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cette fois-ci, la concordance entre la fonction et le polynôme d'interpolation est tout à fait remarquable. Pourquoi ? Nous verrons cela plus loin."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Arrêtons les calculs formels et penchons-nous sur des problèmes numériques. Comment calculer efficacement $P(x)$ pour une valeur donnée de $x$ ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_printing(pretty_print=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. L'algorithme de Neville"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Il existe un certain nombre d'algorithmes permettant le calcul de $P(x)$, où $P$ est le polynôme d'interpolation de Lagagrange. Nous allons en examiner un dans cette partie."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Relations de récurrence sur les polynômes d'interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Proposition__ : Soit $n\\ge 1$. Soit $P$ le polynôme d'interpolation de $f$. Pour $i=0,1,\\ldots,n$ soit $P_{i}$ le polynôme d'interpolation de $f$ aux points $x_0,\\ldots,x_{i-1},x_{i+1},\\ldots,x_n$. Soient $0\\le i,j\\le n$, $i\\ne j$. Soit $x\\in I$. Alors\n",
    "\n",
    "$$P(x)=\\frac{(x-x_j)P_{j}(x)-(x-x_i)P_{i}(x)}{x_i-x_j}$$\n",
    "\n",
    "__Démonstration__ : Posons \n",
    "\n",
    "$$Q=\\frac{(X-x_j)P_{j}-(X-x_i)P_{i}}{x_i-x_j}$$\n",
    "\n",
    "Soit $k\\in[0, n]$. Trois cas se présentent :\n",
    "\n",
    "- $k\\ne i, j$. Alors \n",
    "\n",
    "$$Q(x_k)=\\frac{(x_k-x_j)P_{j}(x_k)-(x_k-x_i)P_{i}(x_k)}{x_i-x_j}=\\frac{(x_k-x_j)f(x_k)-(x_k-x_i)f(x_k)}{x_i-x_j}=f(x_k)$$\n",
    "\n",
    "- $k=i$. On a \n",
    "\n",
    "$$Q(x_i)=\\frac{(x_i-x_j)P_{j}(x_i)-(x_i-x_i)P_{i}(x_i)}{x_i-x_j}=P_{j}(x_i)=f(x_i)$$\n",
    "\n",
    "- $k=j$. Identique au cas précédent.\n",
    "\n",
    "Remarquons enfin que le degré de $Q$ est inférieur ou égal à $n$. Par l'unicité du polynôme d'interpolation, on a $Q=P$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 L'algorithme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Changeons un peu nos notations par rapport au paragraphe précédent et notons bien lourdement $P_{i_1\\ldots i_k}$ le polynôme d'interpolation de $f$ aux points $x_{i_1},\\ldots,x_{i_k}$. Ce qu nous voulons, c'est calculer $P_{01\\ldots n}$.\n",
    "\n",
    "On connaît évidemment $P_j$, polynôme d'interpolation de $f$ en l'unique point $x_j$, qui est constant égal à $f(x_j)$. Dans la fonction `neville` ci-dessous, $Q$ est une liste de $n+1$ flottants. Pour $j=0,\\ldots, n$, on initialise $Q[j]$ à la valeur $f(x_j)$. Suit une boucle indexée par $j$, allant de 1 à $n$ :\n",
    "\n",
    "- Première itération : pour $k = n, n-1, \\ldots, 1$, on affecte à $Q[k]$ la valeur\n",
    "\n",
    "$$\\frac{(x-x_k)Q[k-1]-(x-x_{k-1})Q[k]}{x_{k-1}-x_k}$$\n",
    "\n",
    "qui n'est autre que $P_{(k-1)k}(x)$.\n",
    "\n",
    "- Deuxième itération : pour $k = n, n-1, \\ldots, 2$, on affecte à $Q[k]$ la valeur\n",
    "\n",
    "$$\\frac{(x-x_k)Q[k-1]-(x-x_{k-2})Q[k]}{x_{k-2}-x_k}$$\n",
    "\n",
    "c'est à dire\n",
    "\n",
    "$$\\frac{(x-x_k)P_{(k-2)(k-1)}(x)-(x-x_{k-2})P_{(k-1)k}(x)}{x_{k-2}-x_k}$$\n",
    "\n",
    "qui n'est autre que $P_{(k-2)(k-1)k}(x)$.\n",
    "\n",
    "- Troisième itération : pour $k = n, n-1, \\ldots, 3$, on affecte à $Q[k]$ la valeur $P_{(k-3)(k-2)(k-1)k}(x)$.\n",
    "\n",
    "- Et ainsi de suite.\n",
    "\n",
    "À la sortie de la boucle, $Q[n]$ contient $P_{0\\ldots n}(x)=P(x)$. On le renvoie.\n",
    "\n",
    "__Remarque__ : La complexité de la fonction `neville` est clairement en $O(n^2)$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def neville(f, xs, x):\n",
    "    n = len(xs) - 1\n",
    "    Q = (n + 1) * [0]\n",
    "    for j in range(n + 1): Q[j] = f(xs[j])\n",
    "    for j in range(1, n + 1):\n",
    "        for k in range(n, j - 1, -1):\n",
    "            Q[k] = ((x - xs[k]) * Q[k - 1] - (x - xs[k - j]) * Q[k]) / (xs[k - j] - xs[k])\n",
    "    return Q[n]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous pouvons maintenant tenter d'interpoler en un nombre de points plus important qu'au paragraphe précédent. Soyons toutefois raisonnables, une complexité en $O(n^2)$ nous invite à la modestie :-)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = subdi(-2, 2, 100)\n",
    "xs = subdi(-2, 2, 500)\n",
    "ys = [exemple(t) for t in xs]\n",
    "zs = [neville(exemple, ts, x) for x in xs]\n",
    "plt.ylim(-3, 3)\n",
    "plt.plot(xs, ys, 'r')\n",
    "plt.plot(xs, zs,'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aux bords de l'intervalle, les problèmes s'aggravent lorsqu'on augmente le nombre de points. Si vous n'êtes pas convaincus, enlevez dans la cellule ci-dessus la ligne `plt.ylim(-3, 3)`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. L'erreur commise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nous supposons dans cette partie que l'intervalle $I$ de définition de $f$ est un __segment__, et que $f$ est de classe $\\mathcal C^{n+1}$ sur $I$. La fonction $|f^{(n+1)}|$ est continue sur le segment $I$, elle y possède donc un maximum que nous noterons $M_{n+1}$.\n",
    "\n",
    "\n",
    "Rappelons que $n\\in\\mathbb N$, $x_0<x_1<\\ldots<x_n$ sont $n+1$ points de $I$, et $P$ est le polynôme d'interpolation de $f$ aux points $x_i$.\n",
    "\n",
    "Notons enfin $Q=(X-x_0)(X-x_1)\\ldots(X-x_n)$. Le polynôme $Q$ est unitaire, de degré $n+1$ et s'ennule en tous les $x_i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Une majoration de l'erreur"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $x\\in I$ différent des $x_i$. Considérons la fonction $\\varphi:I\\to\\mathbb R$, de classe $\\mathcal C^{n+1}$, définie pour $t\\in I$ par\n",
    " \n",
    "$$\\varphi(t)=f(t)-P(t)-\\frac{f(x)-P(x)}{Q(x)}Q(t)$$\n",
    "\n",
    "On a pour tout $i\\in[0,n]$, $\\varphi(x_i)=0$, et également $\\varphi(x)=0$. La fonction $\\varphi$ s'annule donc au moins $n+2$ fois sur $I$. De plus, $\\varphi$ est dérivable sur $I$. Par le théorème de Rolle $\\varphi'$ s'annule au moins $n+1$ fois sur $I$. On recommence : $\\varphi''$ s'annule au moins $n$ fois sur $I$, $\\ldots$, $\\varphi^{(n+1)}$ s'annule au moins 1 fois sur $I$ en un point que nous noterons $\\xi$.\n",
    "\n",
    "Comme le polynôme $P$ est de degré inférieur ou égal à $n$, on a $P^{(n+1)}=0$. De là, \n",
    "\n",
    "$$0=\\varphi^{(n+1)}(\\xi)=f^{(n+1)}(\\xi)-\\frac{f(x)-P(x)}{Q(x)}Q^{(n+1)}(\\xi)=f^{(n+1)}(\\xi)-(n+1)!\\frac{f(x)-P(x)}{Q(x)}$$\n",
    "\n",
    "En effet, le polynôme $Q$ est unitaire, de degré $n+1$. Sa dérivée d'ordre $n+1$ est donc constante, égale à $(n+1)!$.\n",
    "\n",
    "On déduit de l'égalité précédente\n",
    "\n",
    "$$f(x)-P(x)=\\frac{f^{(n+1)}(\\xi)}{(n+1)!}Q(x)$$\n",
    "\n",
    "__Proposition__ : Pour tout $x\\in I$, \n",
    "\n",
    "$$|f(x)-P(x)|\\le\\frac{M_{n+1}}{(n+1)!}|Q(x)|$$\n",
    "\n",
    "__Démonstration__ : On vient de le faire si $x$ est différent des $x_i$, et c'est évident si $x$ est l'un des $x_i$ puisque dans ce cas les deux côtés de l'inégalité sont nuls."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La majoration que nous venons d'obtenir est intéressante. Le facteur $(n+1)!$ au dénominateur est prometteur, puisqu'il tend très vite vers $+\\infty$ lorsque $n$ augmente. Le facteur $M_{n+1}$ ne dépend que de $f$ et pas des $x_i$. Et le facteur $|Q(x)|$ ne dépend que des $x_i$ et de $x$, et pas de $f$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 La valeur de $M_{n+1}$ pour notre exemple préféré"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Revenons à $f:x\\mapsto \\frac 1 {x^2+1}$, définie sur le segment $[-2,2]$. Pour tout $x\\in I$,\n",
    "\n",
    "$$f(x)=\\frac{1}{2i}\\left(\\frac 1 {x-i}-\\frac 1 {x+i}\\right)$$\n",
    "\n",
    "On vérifie facilement par récurrence sur $n$ que pour tout entier naturel $n$,\n",
    "\n",
    "$$f^{(n)}(x)=\\frac{(-1)^nn!}{2i}\\left(\\frac 1 {(x-i)^{n+1}}-\\frac 1 {(x+i)^{n+1}}\\right)$$\n",
    "\n",
    "On a $|x-i|=\\sqrt{x^2+1}\\ge 1$ et de même $|x+i|\\ge 1$. On en déduit que le module de la grosse parenthèse est au plus égal à 2, et donc \n",
    "\n",
    "$$|f^{(n)}(x)|\\le n!$$\n",
    "\n",
    "__Conclusion peu optimiste__ : $M_{n+1}\\le (n+1)!$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Que vaut $f^{(n)}(0)$ ? Eh bien,\n",
    "\n",
    "$$f^{(n)}(0)=\\frac{(-1)^nn!}{2i}\\left(\\frac 1 {(-i)^{n+1}}-\\frac 1 {i^{n+1}}\\right)$$\n",
    "\n",
    "Lorsque $n=2p$ est impair, cela donne\n",
    "\n",
    "$$f^{(n)}(0)=(-1)^{p+1}n!$$\n",
    "\n",
    "et donc $M_n\\ge n!$.\n",
    "\n",
    "__Conclusion pessimiste__ : Si $n$ est impair, $M_{n+1}=(n+1)!$.\n",
    "\n",
    "Cette valeur de $M_{n+1}$ annihile évidemment tous les bienfaits du facteur $(n+1)!$ au dénominateur dans le majorant de l'erreur commise. Pour la fonction $f$ de notre exemple, nous avons $|f(x)-P(x)|\\le|Q(x)|$, et pas mieux. En tout cas pas mieux de la façon dont nous avons réalisé notre majoration !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Une majoration de $|Q(x)|$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comment majorer $|Q(x)|$ ? Soit $h$ le _pas_ de la subdivision $(x_0,\\ldots,x_n)$, c'est à dire la plus grande des différences $x_{i+1}-x_i$, $i=0,\\ldots, n-1$. Soit $j$ tel que $x_j\\le x \\le x_{j+1}$. On a \n",
    "\n",
    "$$Q(x)=(x-x_0)\\ldots(x-x_{j-1})(x-x_j)(x_{j+1}-x)(x_{j+2}-x)\\ldots(x_{n}-x)$$\n",
    "\n",
    "Maintenant, faites un petit dessin ... On a $x-x_0\\le(j+1)h$, $x-x_1\\le jh$, $\\ldots$, $x-x_{j-1}\\le 2 h$. \n",
    "\n",
    "De même $x_{j+2}-x\\le 2h$, $\\ldots$, $x_{n}-x\\le (n-j)h$. De là,\n",
    "\n",
    "$$Q(x)\\le h^{n-1}(x-x_j)(x_{j+1}-x)(j+1)!(n-j)!$$\n",
    "\n",
    "Petite astuce maintenant : $\\binom {n+1}{j+1}\\ge n+1$, donc $(j+1)!(n-j)!\\le n!$. Enfin, une petite étude de fonction montre que \n",
    "\n",
    "$$(x-x_j)(x_{j+1}-x)\\le \\frac 1 4 (x_{j+1}-x_j)^2\\le \\frac 1 4 h^2$$\n",
    "\n",
    "Ainsi,\n",
    "\n",
    "$$|Q(x)|\\le \\frac {h^{n+1}n!}{4}$$\n",
    "\n",
    "__Corollaire__ : $$|f(x)-P(x)|\\le \\frac {h^{n+1}M_{n+1}}{4(n+1)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4 Subdivisions régulières"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En prenant une subdivision régulière du segment $I=[a, b]$, on a $h=\\frac {b-a}n$, et donc \n",
    "\n",
    "$$|f(x)-P(x)|\\le \\frac {(b-a)^{n+1}M_{n+1}}{4n^{n+1}(n+1)}$$\n",
    "\n",
    "Pour notre fonction d'exemple $x\\mapsto\\frac 1 {x^2+1}$ sur le segment $[-2,2]$, avec $n$ impair, nous avons $b-a=4$, $M_{n+1}=(n+1)!$, et donc\n",
    "\n",
    "$$|f(x)-P(x)|\\le\\frac{4^{n}n!}{n^{n+1}}$$\n",
    "\n",
    "Une petite application de la formule de Stirling permet de trouver un équivalent de notre majorant :\n",
    "\n",
    "$$\\frac{4^{n}n!}{n^{n+1}}\\sim \\left(\\frac {4} e\\right)^n\\sqrt{\\frac{2\\pi}n}$$\n",
    "\n",
    "Lorsque $n$ tend vers l'infini, comme $e < 4$, le majorant tend vers ... $+\\infty$ :-(.\n",
    "\n",
    "Comment peut-on espérer faire mieux ? La fonction $f$ étant ce qu'elle est, nous ne pouvons rien faire pour la valeur de $M_{n+1}$. En revanche, nous pouvons agir sur le polynôme $Q$. Pouvons nous trouver un polynôme $Q$ unitaire de degré $n+1$ tel que $\\max_{[a,b]}|Q|$ soit minimal parmi tous les polynômes unitaires de degré $n+1$ ? Oui, et c'est ce que nous allons maintenant étudier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Tchebychev"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.1 Les polynômes de Tchebychev"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour tout $n\\in\\mathbb N$, considérons la fonction $T_n:[-1,1]\\to\\mathbb R$ définie par $T_n(x)=\\cos(n\\arccos x)$. Clairement, pour tout $x\\in[-1,1]$, $|T_n(x)|\\le 1$. \n",
    "\n",
    "__Exercice__ : Montrer que pour tout $x\\in[-1,1]$, $T_0(x)=1$, $T_1(x)=x$, et pour tout $n\\ge 0$, $T_{n+2}(x)=2xT_{n+1}(x)-T_n(x)$. Indication : poser $x=\\cos t$. Puis remarquer que $\\cos(n+2)t + \\cos(nt) = 2 \\cos t\\cos (n+ 1) t$.\n",
    "\n",
    "__Conséquence__ : $T_n$ est la restriction à $[-1,1]$ d'une fonction polynomiale que nous continuerons à appeler $T_n$. Les polynômes $T_n$ sont donc définis par la récurrence :\n",
    "\n",
    "- $T_0=0$, $T_1=X$.\n",
    "- Pour tout $n\\in\\mathbb N$, $T_{n+2}=2XT_{n+1}-T_n$.\n",
    "\n",
    "On les appelle les __polynômes de Tchebychev de première espèce__."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Exercice__ : N'y a-t-il pas une question \"d'unicité\" de $T_n$ à se poser ?\n",
    "\n",
    "__Exercice__ : Montrer :\n",
    "\n",
    "- Pour tout $n\\ge 0$, $T_n$ est un polynôme de degré $n$.\n",
    "- Pour tout $n\\ge 1$, le coefficient dominant de $T_n$ est $2^{n-1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `tchebychev` renvoie le $n$ième polynôme de Tchebychev."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tchebychev(n, x):\n",
    "    a = 1\n",
    "    b = x\n",
    "    for k in range(n):\n",
    "        a, b = b, 2 * x * b - a\n",
    "    return a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_printing(pretty_print=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tchebs = [expand(tchebychev(k, x)) for k in range(8)]\n",
    "tchebs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_printing(pretty_print=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Traçons quelques courbes ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (8, 8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k in range(len(tchebs)):\n",
    "    if k != 5: tracer_poly(tchebs[k], -1., 1., 'k')\n",
    "    else: tracer_poly(tchebs[k], -1., 1., 'r')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = (10, 6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Le polynôme $T_5$ est tracé en rouge sur la figure. Combien a-t-il de racines dans l'intervalle $[-1,1]$ ? Apparemment 5. Plus généralement, quelles sont les racines de $T_n$ ?\n",
    "\n",
    "__Proposition__ : Pour tout $n\\ge 1$, les racines de $T_n$ sont les réels $x_k=\\cos(\\frac\\pi{2n}+\\frac{k\\pi}{n})$, $k=0,\\ldots,n-1$.\n",
    "\n",
    "__Démonstration__ : Il est évident que ces réels sont des racines distinctes de $T_n$. Comme $T_n$ est de degré $n$, ce sont les seules racines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2 Une propriété de minimalité"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Soit $Q=\\frac 1{2^n} T_{n+1}$. Le polynôme $Q$ est unitaire, de degré $n+1$, et on a\n",
    "\n",
    "$$\\forall x\\in[-1,1], |Q(x)|\\le\\frac 1 {2^n}$$\n",
    "\n",
    "De plus, $Q(1)=\\frac 1{2^n} T_{n+1}(1)=\\frac 1{2^n}\\cos(n\\arccos 1)=\\frac 1{2^n}$. Ainsi,\n",
    "\n",
    "$$\\max_{[-1,1]}|Q|=\\frac 1 {2^n}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Proposition__ : Soit $P$ un polynôme unitaire de degré $n+1$. Supposons \n",
    "\n",
    "$$\\max_{[-1,1]}|P|\\le\\frac 1 {2^n}$$\n",
    "\n",
    "Alors $P=Q$.\n",
    "\n",
    "__Démonstration__ : Soit $P$ un polynôme unitaire de degré $n+1$. Supposons que pour tout $x\\in[-1,1]$ on ait $|P(x)|\\le \\frac 1 {2^n}$. Posons $R=P-Q$. Le polynôme $R$ est de degré inférieur ou égal à $n$, car $P$ et $Q$ sont unitaires et de degré $n+1$. \n",
    "\n",
    "Pour $k=0,1,\\ldots, n+1$, soit $x_k=\\cos \\frac{k\\pi}{n+1}$. On a $T_{n+1}(x_k)=\\cos(k\\pi)=(-1)^k$. De là, $R(x_k)=P(x_k)-\\frac{(-1)^k}{2^n}$. Comme $|P(x_k)|\\le \\frac 1 {2^n}$, il s'ensuit que les réels $R(x_k)$ sont positifs si $k$ est impair et négatifs si $k$ est pair. Par le théorème des valeurs intermédiaires, le polynôme $R$ s'annule au moins $n+1$ fois. Mais le degré de $R$ est inférieur ou égal à $n$, donc $R=0$ et ainsi $P=Q$.\n",
    "\n",
    "__Exercice__ : J'ai un peu triché. Il se pourrait que deux racines successives de $R$ trouvées par le TVI soient égales. Comment rétablir la démonstration ? Indication : racine multiple.\n",
    "\n",
    "__Remarque__ : ce résultat montre que le polynôme $Q$ est, parmi tous les polynômes unitaires de degré $n+1$, celui dont le maximum de la valeur absolue est minimal (si j'ose dire :-))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.3 Subdivisions de Tchebychev"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plaçons-nous sur un segment $I=[a,b]$. Soient $y_0,\\ldots,y_n$ les racines de $T_{n+1}$.\n",
    "\n",
    "Pour $i=0,\\ldots,n$, posons $x_i=m+\\delta y_i$ où $m=\\frac{a+b} 2$ et $\\delta=\\frac{b-a}2$. On vérifie facilement que les $x_i$ forment une subdivision de $[a,b]$. \n",
    "\n",
    "Posons $Q=\\frac {\\delta^{n+1}}{2^{n}}T_{n+1}(\\frac{X-m}\\delta)$. Le polynôme $Q$ est unitaire, et ses racines sont les $x_i$. Ainsi, $Q=(x-x_0)\\ldots(x-x_n)$. De plus, pour tout $x\\in[a,b]$,\n",
    "\n",
    "$$|Q(x)|\\le \\frac {\\delta^{n+1}} {2^{n}}=\\frac {(b-a)^{n+1}} {2^{2n+1}}$$\n",
    "\n",
    "et ainsi\n",
    "\n",
    "$$|f(x)-P(x)|\\le\\frac{(b-a)^{n+1}M_{n+1}}{2^{2n+1}(n+1)!}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def subdi_tchebychev(a, b, n):\n",
    "    m = (a + b) / 2\n",
    "    d = (b - a) / 2\n",
    "    return [m - d * math.cos(math.pi / (n + 1) * (k + 1 / 2)) for k in range(n + 1)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subdi_tchebychev(-2, 2, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voici les polynômes élémentaires de Lagrange pour une subdivision de Tchebychev."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = subdi_tchebychev(-1, 1, 10)\n",
    "ps = [expand(lagrange(k, x, xs)) for k in range(len(xs))]\n",
    "for p in ps:\n",
    "    tracer_poly(p, -1., 1., 'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vous souvenez-vous de ce que nous avions obtenu avec une subdivision régulière ? Le revoici. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = subdi(-1, 1, 10)\n",
    "ps = [expand(lagrange(k, x, xs)) for k in range(len(xs))]\n",
    "for p in ps:\n",
    "    tracer_poly(p, -1., 1., 'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Je ne commenterai pas plus avant mais les polynômes associés à une subdivision de Tchebychev ont l'air \"mieux\" que ceux associés à une subdivision régulière, surtout aux bords de l'intervalle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.5 Retour à l'exemple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Que devient la majoration sur notre exemple ?\n",
    "\n",
    "$$|f(x)-P(x)|\\le\\frac{4^{n+1}(n+1)!}{2^{2n+1}(n+1)!}=2$$\n",
    "\n",
    "Notre majorant, égal à 2, ne tend toujours pas vers 0 lorsque $n$ tend vers l'infini, mais il est au moins borné. \n",
    "\n",
    "__Remarque__ : Depuis le début de ce notebook, aucune majoration n'a été efficace pour notre fonction d'exemple :-).\n",
    "\n",
    "Remarquons toutefois que si nous nous plaçons sur un segment $I=[a,b]$ quelconque, on a\n",
    "\n",
    "$$|f(x)-P(x)|\\le2\\left(\\frac{b-a}4\\right)^{n+1}$$\n",
    "\n",
    "Ce majorant tend vers 0 dès que $b-a<4$. Nous avons donc montré la convergence uniforme des polynômes d'interpolation de notre exemple $f$ vers la fonction $f$ lorsque l'intervalle est de longueur strictement inférieure à 4. Par exemple, sur $[-1.9999, 1.9999]$. Mais pas sur $[-2, 2]$.\n",
    "\n",
    "Peut-on faire mieux ? Oui, en fait il y a convergence uniforme sur $[-2,2]$, mais ceci est une autre histoire.\n",
    "\n",
    "Petite illustration pour terminer ce notebook, notre fonction `exemple` et son polynôme d'interpolation en des points de Tchebychev ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = subdi_tchebychev(-2, 2, 8)\n",
    "xs = subdi(-2, 2, 300)\n",
    "ys = [exemple(t) for t in xs]\n",
    "zs = [neville(exemple, ts, x) for x in xs]\n",
    "plt.plot(xs, ys, 'r')\n",
    "plt.plot(xs, zs, 'k')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "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": 1
}