{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Interpolation polynomiale - Étude du phénomène de Runge" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Importation de packages pour Python" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import numpy.linalg as npl\n", "from scipy import interpolate" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Retour sur le phénomène de Runge\n", "\n", "On a vu dans le cours le phénomène de Runge qui se traduit par une mauvaise interpolation, lorsque l'on augmente le degré du polynôme d'interpolation de Lagrange, de la fonction $f$ définie sur $\\mathbb{R}$ par\n", "$$\n", "f(x) = \\dfrac{1}{1+x^2}.\n", "$$\n", "Le but du TP est d'observer ce phénomène mais aussi de mettre en oeuvre une meilleure répartition des points d'interpolation à l'aide des racines des polynômes de Chebyshev et de constater l'atténuation du phénomène de Runge. Plus particulièrement, vous implémenterez des fonctions permettant de calculer le polynôme d'interpolation de Lagrange par la méthode directe. L'interpolation par la méthode de Lagrange et celle de Newton vous est ensuite donnée. Enfin, à partir de la méthode de Newton, vous ferez l'interpolation en utilisant les points de Chebyshev." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fonction de Runge" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ">**À faire :** Implémenter une fonction **Runge** qui à un vecteur $x$ sous format numpy retourne le vecteur $f(x)$. Puis, tracer de cette fonction sur l'intervalle $[-3,3]$." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def Runge(x):\n", " return" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Construction de points d'interpolation équirépartis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dans un premier temps, nous allons considérer des points d'interpolation uniformément répartis sur l'intervalle d'interpolation $[a,b]$, $a**À faire :** Implémenter une fonction **Interp_Equi** qui prend en arguments d'entrée les valeurs $a$ et $b$ (qui définissent l'intervalle d'interpolation) ainsi que $m$ de manière à retourner un vecteur $x$ de $m$ points d'interpolations équirépartis sur $[a,b]$. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "def Interp_Equi(a,b,m):\n", " return" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Méthode directe de construction d'un polynôme d'interpolation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On rappelle que la méthode directe de construction consiste à simplement poser le système d'équations suivant\n", "$$ p(x_j) = y_j, \\quad 0\\leq j\\leq n,$$\n", "sous la forme d'un système linéaire dont la solution correspond au vecteur $(a_j)_{0\\leq j\\leq n}$ de coefficients du polynôme d'interpolation. Plus précisément, le système s'écrit\n", "$$ \\left[ \n", "\\begin{array}{ccccc}\n", "1 & x_0 & \\cdots & x_0^{n-1} & x_0^n \\\\ \n", "1 & x_1 & \\cdots & x_1^{n-1} & x_1^n \\\\ \n", "\\vdots & \\vdots & & & \\vdots \\\\ \n", "1 & x_n & \\cdots & x_{n}^{n-1} & x_{n}^n\n", "\\end{array}\n", "\\right]\\left[\\begin{array}{c}a_0\\\\ a_1\\\\ \\vdots \\\\ a_n \\end{array}\\right] = \\left[\\begin{array}{c}y_0\\\\ y_1\\\\ \\vdots \\\\ y_n \\end{array}\\right].\n", " $$\n", " \n", " >**À faire :** Implémenter une fonction **Vandermonde** qui prend en arguments d'entrée un vecteur $x$ de points d'interpolation et qui rend, en sortie, la matrice de Vandermonde associée au système linéaire précédent. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def Vandermonde(x):\n", " return" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " >**À faire :** À l'aide de la fonction **npl.solve**, implémenter une fonction **Methode_directe** qui prend en argument les points d'interpolations $x$ ainsi qu'une fonction $f$ à interpoler et qui rend en sortie un vecteur $a$ qui contient les coefficients du polynôme d'interpolation de $f$ aux points $x$." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def Methode_directe(x,f):\n", " return" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " >**À faire :** Grâce aux fonctions **PtsInterp_Equi** et **Methode_directe**, tracer sur une même figure la fonction **Runge** sur l'intervalle $[-3,3]$ ainsi que son polynôme d'interpolation pour $5$, $10$ et $15$ points d'interpolation." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ " >**À faire :** Utiliser à présent une boite noire de python pour retrouver ces résultats. Pour cela, s'appuyer sur la méthode interpolate.KroghInterpolator de scipy." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Méthode de Lagrange" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "Nous avons vu dans le cours que le polynôme d'interpolation de Lagrange peut être calculé par la méthode de Lagrange. Cette méthode consiste à trouver une base de polynômes, qui sont les polynômes de Lagrange, dans laquelle on exprime le polynôme d'interpolation. Étant donné une fonction $f$ et un ensemble $(x_j)_{0\\leq j\\leq n}$ de points d'interpolation, le polynôme d'interpolation $p$ s'exprime comme\n", "$$\n", "p(x) = \\sum_{j = 0}^{n} f(x_j) L_{j,n}(x), \\quad \\forall x \\in \\mathbb{R},\n", "$$\n", "avec $(L_{j,n})_{0\\leq j\\leq n}$ la base de polynôme de Lagrange donnée par\n", "$$\n", "L_{j,n}(x) = \\prod_{\\substack{k = 0\\\\k\\neq j}}^{n} \\dfrac{x-x_k}{x_j-x_k}.\n", "$$\n", "\n", " >**À faire :** Écrire une fonction **Methode_Lagrange** qui implémente la méthode précédente avec en arguments d'entrée un vecteur $x\\_interp$ de points d'interpolation et un vecteur $x$ de points d'abscisse. Cette fonction returnera la matrice $L$ dont chaque colonne correspond au vecteur $\\ell_i = L_{i,n}(x)$. \n", " \n", "On obtient alors le vecteur $y = p(x)$ grâce au produit matrice-vecteur $y = Lz$ où $z = f(x\\_interp)$. Le produit matrice-vecteur peut se faire grâce à la fonction **np.dot**." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def Methode_Lagrange(x_interp,x):\n", " m = len(x_interp)\n", " n = len(x)\n", " L = np.ones((n,m))\n", " for i in range(m):\n", " for j in range(m):\n", " if j != i:\n", " L[:,i] = L[:,i]*(x - x_interp[j])/(x_interp[i] - x_interp[j])\n", " return L" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " >**À faire :** Tracer ensuite sur une même figure la fonction **Runge** sur l'intervalle $[-3,3]$ ainsi que son polynôme d'interpolation obtenue avec la méthode de Lagrange pour $5$, $10$ et $15$ points d'interpolation." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "ename": "NameError", "evalue": "name 'np' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mRunge\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mn\u001b[0m \u001b[0;32min\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m15\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mx_interp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mInterp_Equi\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mz\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mRunge\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx_interp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'np' is not defined" ] } ], "source": [ "x = np.linspace(-3,3,1000)\n", "plt.plot(x,Runge(x))\n", "for n in [5,10,15]:\n", " x_interp = Interp_Equi(-3,3,n)\n", " z = Runge(x_interp)\n", " L = Methode_Lagrange(x_interp,x)\n", " P = L.dot(z)\n", " plt.plot(x,P)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "### Méthode de Newton\n", "\n", "La méthode de Newton, comme pour celle de Lagrange, permet d'exprimer le polynôme d'interpolation de Lagrange dans une base de polynômes, qui sont les polynômes de Newton. Étant donné une fonction $f$ et un ensemble $(x_j)_{0\\leq j\\leq n}$ de points d'interpolation, le polynôme d'interpolation $p$ s'exprime comme\n", "$$\n", "p(x) = \\sum_{j = 0}^{n} z_j N_{j,n}(x), \\quad \\forall x \\in \\mathbb{R},\n", "$$\n", "avec $(N_{j,n})_{0\\leq j\\leq n}$ la base de polynôme de Newton donnée par\n", "$$\n", "N_{j,n}(x) = \\prod_{k = 0}^{j} (x-x_k),\n", "$$\n", "et $(z_j)_{0\\leq j\\leq n}$ les coefficients qui sont solution du système triangulaire\n", "$$\n", "\\left[ \n", "\\begin{array}{cccccc}\n", "1 & 0 & 0 & \\cdots & 0 \\\\ \n", "1 & (x_1-x_0) & 0 & \\cdots & 0 \\\\ \n", "1 & (x_2-x_0) & (x_2-x_0)(x_2-x_1) & \\ddots & \\vdots \\\\ \n", "\\vdots & \\vdots & \\vdots & \\ddots & 0 \\\\ \n", "1 & (x_n-x_0) & (x_n-x_0)(x_n-x_1) & \\cdots & \\prod_{j = 0}^{n-1}(x_n-x_j)\n", "\\end{array}\n", "\\right]\\left[\\begin{array}{c}z_0\\\\ z_1\\\\ \\vdots \\\\ z_n \\end{array}\\right] = \\left[\\begin{array}{c}y_0\\\\ y_1\\\\ \\vdots \\\\ y_n \\end{array}\\right].\n", "$$\n", "\n", " >**À faire :** Écrire une fonction **Methode_Newton** qui implémente la méthode précédente avec en arguments d'entrée un vecteur $x\\_interp$ de points d'interpolation et un vecteur $x$ de points d'abscisse. Cette fonction returnera la matrice $N$ dont chaque colonne correspond au vecteur $n_i = N_{i,n}(x)$. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def Methode_Newton(x_interp,x):\n", " m = len(x_interp)\n", " n = len(x)\n", " N = np.ones((n,m))\n", " for i in range(m):\n", " for j in range(i):\n", " N[:,i] = N[:,i]*(x - x_interp[j])\n", " return N" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ">**À faire :** Écrire une fonction **Coeff_Newton** qui permet de calculer les coefficients $(z_j)_{0\\leq j\\leq n}$ en résolvant le système triangulaire inférieur par un algorithme de descente. Cette fonction prendra en entrée le vecteur $x\\_interp$ des points d'interpolation ainsi que la fonction $f$ à interpoler et rendra le vecteur $z$ des coefficients." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def Coeff_Newton(x_interp,f):\n", " m = len(x_interp)\n", " T = np.zeros((m,m))\n", " for i in range(m):\n", " for j in range(m):\n", " if j == 0:\n", " T[i,j] = 1\n", " elif j<=i:\n", " T[i,j] = T[i,j-1]*(x_interp[i] - x_interp[j-1])\n", " z = npl.solve(T,f(x_interp))\n", " return z" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On obtient alors le vecteur $y = p(x)$ grâce au produit matrice-vecteur $y = Lz$.\n", "\n", ">**À faire :** Tracer ensuite sur une même figure la fonction **Runge** sur l'intervalle $[-3,3]$ ainsi que son polynôme d'interpolation obtenue avec la méthode de Newton pour $5$, $10$ et $15$ points d'interpolation." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "ename": "NameError", "evalue": "name 'np' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mRunge\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mn\u001b[0m \u001b[0;32min\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m15\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mx_interp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mInterp_Equi\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mz\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mCoeff_Newton\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx_interp\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mRunge\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'np' is not defined" ] } ], "source": [ "x = np.linspace(-3,3,1000)\n", "plt.plot(x,Runge(x))\n", "for n in [5,10,15]:\n", " x_interp = Interp_Equi(-3,3,n)\n", " z = Coeff_Newton(x_interp,Runge)\n", " N = Methode_Newton(x_interp,x)\n", " P = N.dot(z)\n", " plt.plot(x,P)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Points d'interpolation de Chebyshev" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Les points de Chebyshev dans l'intervalle $[-1,1]$ sont donnés par la formule suivante\n", "$$\n", "y_j = \\cos\\left(\\dfrac{2i+1}{2(n+1)}\\pi\\right), \n", "\\quad 0\\leq j\\leq n.\n", "$$\n", "Afin d'adapter ces points à un intervalle $[a,b]$, on se basera sur la formule suivante\n", "$$\n", "x_j = \\frac{a+b}2 + \\frac{a-b}2 y_j,\n", "\\quad 0\\leq j\\leq n.\n", "$$\n", "\n", ">**À faire :** Implémenter une fonction **Interp_Chebyshev** qui prend en arguments d'entrée les valeurs $a$ et $b$ (qui définissent l'intervalle d'interpolation) ainsi que $m$ de manière à retourner un vecteur $x$ de $m$ points d'interpolations répartis selon les points de Chebyshev sur $[a,b]$. " ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def Interp_Chebyshev(a,b,m):\n", " return" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ ">**À faire :** Tracer ensuite sur une même figure la fonction **Runge** sur l'intervalle $[-3,3]$ ainsi que son polynôme d'interpolation obtenue avec la méthode de Newton pour $5$, $10$ et $15$ points d'interpolation." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "### Comparaison d'erreur selon la répartition des points d'interpolation\n", "\n", "Nous sommes à présent en mesure de quantifier numériquement la réduction de l'erreur qu'apporte la répartition suivant les points de Chebyshev en comparaison avec les points équirépartis. On introduit pour cela l'erreur\n", "$$\n", "e(n) = \\sup_{x\\in[-3,3]} |f(x) - p(x)|,\n", "$$\n", "où $f$ est la fonction à interpoler et $p$ est le polynôme d'interpolation de Lagrange de degré $n$. On veut comparer l'évolution de l'erreur en fonction de $n$ et du choix de la répartition des points d'interpolation.\n", "\n", ">**À faire :** Écrire une fonction **Erreur_Interp** permettant de calculer l'erreur $e(n)$ dans le cas de points d'interpolation équirépartis et le cas des points de Chebyshev. Cette fonction prendra en entrée l'entier $n$ ainsi qu'un entier $p$ et donnera en sortie un scalaire correspondant à l'erreur $e(n)$ pour les points équirépartis si $p = 0$ ou pour les points de Chebyshev si $p = 1$. L'interpolation se fera à l'aide de la méthode de Newton sur l'intervalle $[-3,3]$." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def Erreur_Interp(n,p):\n", " return" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ ">**À faire :** Tracer sur une même figure l'évolution de l'erreur $\\log(e(n))$, en fonction de $n$, pour les points d'interpolation équirépartis et pour les points de Chebyshev de $1$ jusqu'à $100$ points." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": true }, "outputs": [], "source": [ "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.11" } }, "nbformat": 4, "nbformat_minor": 1 }