{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.4 Interpolation par intervalles\n",
    "\n",
    "### Interpolation linéaire par morceaux\n",
    "\n",
    "Soit $f: \\mathbf [a,b] \\rightarrow R$ continue et $a=x_0<\\ldots<x_n=b$. \n",
    "\n",
    "On choisit une partition de $[a,b]$ en $N$ sous-intervalles de la forme $[x_i, x_{i+1}], i=0,...,N-1]$.\n",
    "Sur chaque sous-intervalle, on fait une interpolation de degré 1 avec les 2 noeuds $x_i, x_{i+1}$.\n",
    "Sur chaque sous-intervalle $I_i=[x_i, x_{i+1}]$, on interpole $f_{\\mid I_i}$ par\n",
    "un polynôme de degré 1. Le polynôme par morceaux (polynôme\n",
    "  composite) qu'on obtient est noté $\\Pi_1^H f(x)$ et on a:\n",
    "$$\\Pi_1^H f(x) = f(x_i) + \\dfrac{f(x_{i+1}) - f(x_i)}{x_{i+1} - x_i}\n",
    "(x - x_i) \\quad \\text{pour } x\\in [x_i, x_{i+1}]$$\n",
    "\n",
    "Dans le cas de données $y$, cela s'écrit\n",
    "$$\\Pi_1^H (x) = y_i + \\dfrac{y_{i+1} - y_i}{x_{i+1} - x_i}\n",
    "(x - x_i) \\quad \\text{pour } x\\in [x_i, x_{i+1}]$$\n",
    "\n",
    "Le choix le plus simple est le suivant :\n",
    "\n",
    "* on pose $H = \\frac{b-a}{N}$\n",
    "* ensuite $x_i = a + iH$ pour $i=0,...,N$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# importing libraries used in this book\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from InterpolationLib import PiecewiseLinearInterpolation as  PiH1\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fonction `PiH1` implémente l'interpolation linéaire par morceaux pour des points équidistribués\n",
    "\n",
    "```python\n",
    "def PiecewiseLinearInterpolation(a,b,N,f,z):\n",
    "    # the input variables are:\n",
    "    # a,b : x[0] = a, x[n] = b\n",
    "    # f : the corresponding data at the points x\n",
    "    # z : where to evaluate the function\n",
    "```\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercice\n",
    "\n",
    "Soient $x_k, k=0,...,4$ des points équidistribués sur l'intervalle $[1,5]$ et \n",
    "$y_k = (3.38, 3.86, 3.85, 3.59, 3.49)$\n",
    "les valeurs d'une fonction en ces points. \n",
    "\n",
    "* Dessinez le graphe de l'interpolateur par morceaux de cette fonction\n",
    "* Calculez numériquement (sur papier) la valeur de $\\Pi^H_1(4.5)$ et vérifiez le résultat sur le graphique\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# intervalle d'interpolation\n",
    "a = 1; b = 5\n",
    "\n",
    "# define the vector of values at equidistributed points x\n",
    "# N the number of intervals\n",
    "# y the array of values\n",
    "# then use PiH1(a,b,N,y,z) to compute and plot the piecewise linear interpolant\n",
    "\n",
    "# >>> COMPLETE HERE <<<\n",
    "\n",
    "\n",
    "\n",
    "plt.xlabel('x'); plt.ylabel('y');\n",
    "plt.legend(['data','p(x)'])\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercice\n",
    "\n",
    "Dessinez le graphe de la fonction de Runge et l'interpolateur linéaire par morceaux \n",
    "sur l'intervalle $[-5,5]$ pour $N=3,5,10$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# interval and function\n",
    "a = -5; b = 5\n",
    "f = lambda x : 1/(1+x**2)\n",
    "\n",
    "# >>> COMPLETE HERE <<<\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Erreur d'interpolation linéaire par morceaux \n",
    "\n",
    "#### Théorème \n",
    "\n",
    "Soient $f\\in C^{2}([a,b])$, $H = \\frac{b-a}{N}$, $x_i = a + iH$ pour $i=0,...,N$.\n",
    "  \n",
    "Soit $E^H_1f(x) = f(x) - \\Pi_1^{H} f(x)$, alors\n",
    "$$\\max_{x\\in I}\\mid E^H_1f(x) \\mid\\leq \\dfrac{H^2}{8}\\max_{x \\in I} | f''(x)|.$$\n",
    "\n",
    "**Preuve**\n",
    "D'après la formule, sur chaque intervalle $I_i=[x_i,x_{i+1}]$, on a\n",
    "  \n",
    "$$\\max_{x\\in[x_i,x_{i+1}]}\\mid E^H_1f(x)\\mid \\leq\n",
    "      \\dfrac{H^2}{4(1+1)}\\max_{x\\in I_i}\\mid f''(x)\\mid.$$\n",
    "\n",
    "**Remarque**\n",
    "On peut aussi montrer que, si l'on utilise un polynôme\n",
    "de degré $n$ ($\\geq 1$) et si l'on dénote $E^H_nf(x) = f(x) - \\Pi_n^{H} f(x)$,\n",
    " dans chaque sous-intervalle $I_i$, on\n",
    "trouve\n",
    "\n",
    "$$\\max_{x\\in I}\\mid E^H_nf(x) \\mid \\leq \\dfrac{H^{n+1}}{4(n+1)} \\max_{x \\in I} |f^{(n+1)} (x)| \\, .$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpolation quadratique par morceaux\n",
    "\n",
    "Soit $f: \\mathbf [a,b] \\rightarrow R$ continue et $a=x_0<\\ldots<x_n=b$. \n",
    "Sur chaque sous-intervalle $[x_i,x_{i+1}]$, on fait une interpolation de **degré 2** avec \n",
    "les **3 noeuds** $x_i, x_{i+\\frac12}, x_{i+1}$, où $x_{i+\\frac12}$ est le milieu de $[x_i,x_{i+1}]$.\n",
    "Sur chaque sous-intervalle $I_i=[x_i, x_{i+1}]$, on interpole $f_{\\mid I_i}$ par\n",
    "un polynôme de degré 2. Le polynôme par morceaux (polynôme\n",
    "  composite) qu'on obtient est noté $\\Pi_2^H f(x)$ et on a:\n",
    "\n",
    "$$\\Pi_2^H f(x) = f(x_i)\\varphi^{(i)}_0(x) + f(x_{i+\\frac12})\\varphi^{(i)}_1(x) \n",
    "    + f(x_{i+1})\\varphi^{(i)}_2(x) \\quad \\text{pour } x\\in [x_i, x_{i+1}]$$\n",
    "\n",
    "où $\\varphi^{(i)}_0(x), \\varphi^{(i)}_1(x), \\varphi^{(i)}_0(x)$ sont les polynômes de la base de Lagrange\n",
    "associés aux noeuds $(x_i,x_{i+\\frac12}, x_{i+1})$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from InterpolationLib import PiecewiseQuadraticInterpolation as  PiH2\n",
    "\n",
    "# intervalle et fonction\n",
    "a = -5; b = 5\n",
    "f = lambda x : 1/(1+x**2)\n",
    "\n",
    "# Values of N to use\n",
    "Nrange = [3,5,10]\n",
    "# plotting points\n",
    "z = np.linspace(-5, 5, 100)\n",
    "\n",
    "for N in Nrange :\n",
    "    plt.plot(z, PiH2(a,b,N,f,z), ':')\n",
    "\n",
    "plt.plot(z, f(z), 'b')\n",
    "plt.xlabel('x'); plt.ylabel('y'); #plt.title('data')\n",
    "plt.legend(Nrange)\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercice\n",
    "\n",
    "Calculez l'erreur d'interpolation de la fonction de Runge dans l'intervalle $[-5,5]$ quand on utilise l'interpolation linéaire par morceaux\n",
    "\n",
    "1. Théoriquement\n",
    "2. Faites un graphique qui montre la convergence en fonction de $H=\\frac{b-a}N=\\frac{10}{N}$\n",
    "3. Refaites le même exercice pour l'interpolation quadratique par morceaux"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Remarque\n",
    "\n",
    "On s'attend à ce que l'erreur soit quadratique en H. Pour le voir il faut utiliser \n",
    "des axes logarithmiques dans les deux directions.\n",
    "\n",
    "Admettons que l'erreur converge quadratiquement vers 0 par rapport à $H$, e.g\n",
    "```Python\n",
    "err = lambda H : 3*H**2 + 2*H**3;\n",
    "\n",
    "```\n",
    "Comment est le graphique de cette fonction dans l'intervalle $[10^{-6},1]$ ?\n",
    "Que se passe-t-il si on change les axes avec `plt.xscale('log')` et  `plt.yscale('log')` ?\n",
    "Quelle est la pente de `err` dans ce système ?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Assume the error is quadratic\n",
    "err = lambda H : 10*H**2 + 20*H**3;\n",
    "\n",
    "# take h has powers of 2: 2^(-20),2^(-19), ..., 2^(-1) \n",
    "h = np.power(2,np.linspace(-20, -1, 20, endpoint=True) )\n",
    "\n",
    "plt.plot(h, err(h), 'b:.')\n",
    "\n",
    "plt.xlabel('h'); plt.ylabel('err');\n",
    "\n",
    "# plt.xscale('log') \n",
    "# plt.yscale('log')\n",
    "\n",
    "plt.grid(True)\n",
    "# try with and wothout equal axis\n",
    "plt.axis('equal')\n",
    "plt.show()\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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}