6_interpolation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Interpolation\n",
    "\n",
    "**Definition:**  Given a discrete set of values $y_i$ at locations $x_i$, an *interpolant* is a (piece-wise) continuous function $f(x)$ that passes exactly through the data (*i.e.* $f(x_i) = y_i$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Applications\n",
    "\n",
    " - Data filling\n",
    " - Function approximation\n",
    " - Fundamental component of other algorithms\n",
    "   - Root finding (secant method)\n",
    "   - Optimization, minima/maxima (successive parabolic interpolation)\n",
    "   - Numerical integration and differentiation\n",
    "   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Polynomial Interpolation (1-D)\n",
    "\n",
    "**Theorem:**  There is a *unique* polynomial of degree $N$, $P_N(x)$, that passes exactly through $N + 1$ values $y_1, y_2, \\ldots, y_N, y_{N+1}$ at *distinct* points $x_1, x_2, \\ldots, x_N, x_{N+1}$.\n",
    "\n",
    "Consequence of the number of unknowns in $P_N(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example 1: 2 Points\n",
    "\n",
    "Given points are $(x_0, y_0)$ and $(x_1, y_1)$ which will lead to a line:\n",
    "\n",
    "Define $P_1(x) = p_1 x + p_0$ and use the two points to find $p_0$ and $p_1$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$y_0 = p_1 x_0 + p_0 ~~~~ \\Rightarrow ~~~~ p_0 = y_0 - p_1 x_0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$y_1 = p_1 x_1 + p_0 ~~~~ \\Rightarrow ~~~~ y_1 = p_1 x_1 + y_0 - p_1 x_0 ~~~~ \\Rightarrow ~~~~ p_1 = \\frac{y_1 - y_0}{x_1 - x_0} ~~~~ \\Rightarrow ~~~~ p_0 = y_0 - \\frac{y_1 - y_0}{x_1 - x_0} x_0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$P_1(x) = \\frac{y_1 - y_0}{x_1 - x_0} x + y_0 - \\frac{y_1 - y_0}{x_1 - x_0} x_0 = \\frac{y_1 - y_0}{x_1 - x_0} (x - x_0) + y_0$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example 2:  3 Points\n",
    "\n",
    "Given points are $(x_0, y_0)$, $(x_1, y_1)$, and $(x_2, y_2)$ which will lead to quadratic polynomial:\n",
    "\n",
    "Define $P_2(x) = p_0 x^2 + p_1 x + p_2$ leading to the equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$y_0 = p_2 x_0^2 + p_1 x_0 + p_0$$\n",
    "\n",
    "$$y_1 = p_2 x_1^2 + p_1 x_1 + p_0$$\n",
    "\n",
    "$$y_2 = p_2 x_2^2 + p_1 x_2 + p_0$$\n",
    "\n",
    "This gets complicated quickly!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Proof\n",
    "\n",
    "Let \n",
    "\n",
    "$$P_N(x) = \\sum^N_{n=0} p_n x^n $$\n",
    "or\n",
    "$$P_N(x) = p_0 + p_1 x + \\cdots + p_{N - 1} x^{N - 1} + p_{N} x^N$$\n",
    "\n",
    "and require $P_N(x_i) = y_i$ for $i=0,1,\\ldots,N$ and $x_i \\neq x_j ~~~ \\forall i,j$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Assume there exists another polynomial \n",
    "\n",
    "$$Q_N(x) = \\sum^N_{n=0} q_n x^n$$\n",
    "\n",
    "that passes through the same set of points such that $Q_N(x_i) = y_i$.  Now compute $T_N(x) = P_N(x) - Q_N(x)$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We know that by assumption that $T_N(x_i) = 0$ but what about for all $x$?\n",
    "\n",
    "$$T_N(x) = P_N(x) - Q_N(x) = \\sum^N_{n=0} p_n x^n - q_n x^n = \\sum^N_{n=0} (p_n - q_n) x^n$$\n",
    "\n",
    "But if $T_N(x_i) = 0$ implies that $p_n - q_n = 0$ individually and therefore $P_N(x) = Q_N(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Monomial Basis\n",
    "\n",
    "Let $P_N(x) = \\sum^N_{n=0} p_n x^n$, $P_N(x)$ can be represented by a *linear combination* of the monomials\n",
    "\n",
    "$$1, x, x^2, x^3, \\ldots, x^{N-1}, x^N$$ \n",
    "\n",
    "with weights \n",
    "\n",
    "$$p_0, p_1, p_2, p_3, \\ldots, p_{N-1},~~ \\text{and}~~ p_N$$\n",
    "\n",
    "respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example 3:  Monomial Basis\n",
    "\n",
    "Consider $P_3(x) = p_0 + p_1 x + p_2 x^2 + p_3 x^3$ with the four data points $(x_i, y_i), ~~ i = 0,1,2,3$.  We have four equations and four unknowns as expected:\n",
    "\n",
    "$$P_3(x_0) = p_0 + p_1 x_0 + p_2 x_0^2 + p_3 x_0^3 = y_0$$\n",
    "\n",
    "$$P_3(x_1) = p_0 + p_1 x_1 + p_2 x_1^2 + p_3 x_1^3 = y_1$$\n",
    "\n",
    "$$P_3(x_2) = p_0 + p_1 x_2 + p_2 x_2^2 + p_3 x_2^3 = y_2$$\n",
    "\n",
    "$$P_3(x_3) = p_0 + p_1 x_3 + p_2 x_3^2 + p_3 x_3^3 = y_3$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Lets rewrite these as a matrix equation:\n",
    "\n",
    "$$\\vec{x} = \\begin{bmatrix} x_0 \\\\ x_1 \\\\ x_2 \\\\ x_3 \\end{bmatrix} ~~~~ \\vec{y} = \\begin{bmatrix} y_0 \\\\ y_1 \\\\ y_2 \\\\ y_3 \\end{bmatrix} ~~~~ \\vec{p} = \\begin{bmatrix} p_0 \\\\ p_1 \\\\ p_2 \\\\ p_3 \\end{bmatrix}$$\n",
    "\n",
    "Define the *Vandermonde* matrix as\n",
    "\n",
    "$$\n",
    "V = \\begin{bmatrix} \n",
    "     1 & x_0 & x_0^2 & x_0^3 \\\\\n",
    "     1 & x_1 & x_1^2 & x_1^3 \\\\\n",
    "     1 & x_2 & x_2^2 & x_2^3 \\\\\n",
    "     1 & x_3 & x_3^2 & x_3^3\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "which allows us to write the system of linear equations as $V \\vec{p} = \\vec{y}$:\n",
    "\n",
    "$$\\begin{bmatrix} \n",
    "     1 & x_0 & x_0^2 & x_0^3 \\\\\n",
    "     1 & x_1 & x_1^2 & x_1^3 \\\\\n",
    "     1 & x_2 & x_2^2 & x_2^3 \\\\\n",
    "     1 & x_3 & x_3^2 & x_3^3\n",
    "\\end{bmatrix} \\begin{bmatrix} p_0 \\\\ p_1 \\\\ p_2 \\\\ p_3 \\end{bmatrix} = \\begin{bmatrix} y_0 \\\\ y_1 \\\\ y_2 \\\\ y_3 \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Vandermonde matrices in general are defined as\n",
    "\n",
    "$$V = \\begin{bmatrix} \n",
    "    1      & x_0    & x_0^2  & \\cdots & x_0^N  \\\\\n",
    "    1      & x_1    & x_1^2  & \\cdots & x_1^N  \\\\\n",
    "    \\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "    1      & x_m    & x_m^2  & \\cdots & x_m^N  \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "where $V$ is a $m \\times n$ matrix with points $(x_i, y_i)$ for $i = 0, 1, 2, 3, \\ldots m$ and for an order $N$ polynomial $P_N(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Finding $p_i$\n",
    "\n",
    "Finding the coefficients of $P_N(x)$ can be done by solving the system outlined above but there are other approaches:\n",
    "\n",
    " - Use `numpy.polyfit(x, y, x.shape[0] - 1)`\n",
    " - Use `numpy.vander(x, N=None)` to construct the matrix and use a linear solver routine.\n",
    " - Use a Lagrangian basis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Lagrangian Basis\n",
    "\n",
    "Given $N+1$ points $(x_0,y_0), (x_1,y_1), \\ldots, (x_{N},y_{N})$ again assuming the $x_i$ are all unique, the interpolating polynomial $P_N(x)$ can be written as\n",
    "\n",
    "$$P_N(x) = \\sum^{N}_{i=0} y_i \\ell_i(x)$$\n",
    "\n",
    "where \n",
    "\n",
    "$$\\ell_i(x) = \\prod^{N}_{j=0, j \\neq i} \\frac{x - x_i}{x_j - x_i} = \\frac{x - x_0}{x_i - x_0} \\frac{x - x_1}{x_i - x_1} \\cdots \\frac{x - x_{i-1}}{x_i - x_{i-1}}\\frac{x - x_{i+1}}{x_i - x_{i+1}} \\cdots \\frac{x - x_{N}}{x_i - x_{N}}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example 4:  $N = 1$ Lagrange Polynomial\n",
    "\n",
    "Given 2 points $(x_0, y_0)$ and $(x_1, y_1)$ the Lagrange form of $P_N(x)$ is given by\n",
    "\n",
    "$$\\ell_0(x) = \\frac{x - x_1}{x_0 - x_1}$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\\ell_1(x) = \\frac{x - x_0}{x_1 - x_0}$$\n",
    "\n",
    "so that\n",
    "\n",
    "$$P_1(x) = \\ell_0(x) \\cdot y_0 + \\ell_1(x) \\cdot y_1 = \\frac{x - x_1}{x_0 - x_1} \\cdot y_0 + \\frac{x - x_0}{x_1 - x_0} \\cdot y_1$$\n",
    "\n",
    "One important aspect of Lagrange polynomials to note is that the $\\ell_i(x)$ functions are exactly 1 when $x  = x_i$ and that every other $\\ell_j(x)$ where $j \\neq i$ is 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "data = numpy.array([[-1.5, -0.5], [0.0, 0.0]])\n",
    "N = data.shape[0] - 1\n",
    "M = data.shape[0]\n",
    "x = numpy.linspace(-2.0, 2.0, 100)\n",
    "\n",
    "# ====================================================\n",
    "# Compute the Lagrange basis (\\ell_i(x))\n",
    "lagrange_basis = numpy.ones((N + 1, x.shape[0]))\n",
    "for i in xrange(N + 1):\n",
    "    for j in xrange(N + 1):\n",
    "        if i != j:\n",
    "            print i, data[i, 0], data[j, 0]\n",
    "            lagrange_basis[i, :] *= (x - data[j, 0]) / (data[i, 0] - data[j, 0])\n",
    "\n",
    "# Calculate full polynomial\n",
    "P = numpy.zeros(x.shape[0])\n",
    "for n in xrange(N + 1):\n",
    "    P += lagrange_basis[n, :] * data[n, 1]\n",
    "\n",
    "# ====================================================\n",
    "# Plot individual basis functions\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "for i in xrange(N + 1):\n",
    "    axes.plot(x, lagrange_basis[i, :], label=\"$\\ell_{%s}(x)$\" % i)\n",
    "#     axes.plot(data[i, 0], data[i, 1], 'ko')\n",
    "\n",
    "axes.set_title(\"Lagrange Basis $\\ell_i(x)$\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$\\ell_i(x)$\")\n",
    "axes.legend(loc=8)\n",
    "\n",
    "# Plot full polynomial P_N(x)\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, P, label=\"$P_{%s}(x)$\" % N)\n",
    "for point in data:\n",
    "    axes.plot(point[0], point[1], 'ko')\n",
    "axes.set_title(\"$P_N(x)$\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$P_N(x)$\")\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example 5:  Interpolate six points from $sin(\\pi x)$\n",
    "\n",
    "Use six points to approximate $\\sin$ on the interval $x \\in [-1, 1]$.  What is the behavior as $N \\rightarrow \\infty$?  Also plot the error between $f(x)$ and the interpolant $P_N(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example 6: Runge's Function\n",
    "\n",
    "Interpolate $f(x) = \\frac{1}{1 + 25 x^2}$ using 6 points of your choosing on $x \\in [-1, 1]$.\n",
    "\n",
    "Try it with 11 points.\n",
    "\n",
    "Keep increasing the number of points and see what happens."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return 1.0 / (1.0 + 25.0 * x**2)\n",
    "\n",
    "x = numpy.linspace(-1, 1, 100)\n",
    "\n",
    "data = numpy.empty((6, 2))\n",
    "data[:, 0] = numpy.linspace(-1, 1, 6)\n",
    "data[:, 1] = f(data[:, 0])\n",
    "N = data.shape[0] - 1\n",
    "\n",
    "# Calculate interpolant\n",
    "lagrange_basis = numpy.ones((N + 1, x.shape[0]))\n",
    "for i in xrange(N + 1):\n",
    "    for j in xrange(N + 1):\n",
    "        if i != j:\n",
    "            lagrange_basis[i, :] *= (x - data[j, 0]) / (data[i, 0] - data[j, 0])\n",
    "\n",
    "# Calculate full polynomial\n",
    "P = numpy.zeros(x.shape[0])\n",
    "for n in xrange(N + 1):\n",
    "    P += lagrange_basis[n, :] * data[n, 1]\n",
    "    \n",
    "# Plot the results\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(x, P, 'b', label=\"$P_6(x)$\")\n",
    "axes.plot(x, f(x), 'k', label=\"True $f(x)$\")\n",
    "axes.plot(data[:, 0], data[:, 1], 'ro', label=\"data\")\n",
    "          \n",
    "axes.set_title(\"Interpolation of Runge's function\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"y\")\n",
    "axes.legend(loc=1)\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Rules of Thumb\n",
    "\n",
    " - Avoid high-order interpolants when possible!  Keep increasing the number of points and see what happens.\n",
    " - Avoid extrapolation - Increase the range of $x$ in the above example and check how good the approximation is beyond our sampling interval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Error Analysis\n",
    "\n",
    "**Theorem:** Lagrange Remainder Theorem - Let $f(x) \\in C^{N+1}[-1, 1], then\n",
    "\n",
    "$$f(x) = P_N(x) + R_N(x)$\n",
    "\n",
    "where $P_N(x)$ is the interpolating polynomial and\n",
    "\n",
    "$$R_N(x) = Q(x) \\frac{f^{(N+1)}(c)}{(N+1)!} ~~~~ \\text{with}~~~~ c \\in [-1,1]$$\n",
    "\n",
    "and\n",
    "\n",
    "$$Q(x) = (x-x_0)(x-x_1)\\cdots(x-x_N).$$\n",
    "\n",
    "A few things to note:\n",
    " - For Taylor's theorem note that $Q(x) = (x - x_0)^{N+1}$ and the error only vanishes at $x_0$.\n",
    " - For Lagrange's theorem the error vanishes at all $x_i$.\n",
    " - To minimize $R_N(x)$ requires minimizing $|Q(x)|$ for $x \\in [-1, 1]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Chebyshev Polynomials\n",
    "\n",
    "*Chebyshev polynomials* $T_N(x)$ are another basis that can be used for interpolation.\n",
    "\n",
    "First 5 polynomials\n",
    "$$T_0(x) = 1$$\n",
    "\n",
    "$$T_1(x) = x$$\n",
    "\n",
    "$$T_2(x) = 2 x^2 - 1$$\n",
    "\n",
    "$$T_3(x) = 4 x^3 - 3 x$$\n",
    "\n",
    "$$T_4(x) = 8x^4 - 8x^2 + 1$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(-1, 1, 100)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(x, numpy.ones(x.shape), label=\"$T_0$\")\n",
    "axes.plot(x, x, label=\"$T_1$\")\n",
    "axes.plot(x, 2.0 * x**2 - 1.0, label=\"$T_2$\")\n",
    "axes.plot(x, 4.0 * x**3 - 3.0 * x, label=\"$T_3$\")\n",
    "axes.plot(x, 8.0 * x**4 - 8.0 * x**2 + 1, label=\"$T_4$\")\n",
    "\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "axes.set_title(\"Chebyshev Polynomials\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$T_N(x)$\")\n",
    "axes.legend(loc=1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Properties:\n",
    "1. Defined by a recurrence relation \n",
    "\n",
    "   $$T_k(x) = 2 x T_{k-1}(x) - T_{k-2}(x)$$\n",
    "   \n",
    "1. Leading coefficient of $x^N$ in $T_N(x)$ is $2^{N-1}$ for $N \\geq 1$\n",
    "1. Extreme values:\n",
    "\n",
    "   $$|T_N(x)| \\leq 1 ~~~~\\text{for}~~~~~ -1 \\leq x \\leq 1$$\n",
    "   \n",
    "1. Minimax principle:  The polynomial\n",
    "\n",
    "   $$T(x) = \\frac{T_{N+1}(x)}{2^N}$$\n",
    "   \n",
    "   is a monic polynomial with the property that\n",
    "   \n",
    "   $$\\max |T(x)| \\leq \\max |Q(X)| ~~~ \\text{for} ~~~ x \\in [-1, 1], ~~~~ \\text{and}$$\n",
    "   \n",
    "   $$\\max |T(x)| = \\frac{1}{2^N}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Error Analysis Redux\n",
    "\n",
    "Recall that in order to minimize the error with Lagrange's Theorem we needed to minimize $|Q(x)|$ for $x \\in [-1, 1]$.  This implies that in order to minimize $R_N(x)$ with Chebyshev polynomials set $Q(x) = T(x)$.  We only know $Q(x)$ through its roots $x_0, x_1, \\ldots, x_N$ so we will require that the points used (nodes) in the interpolation by the zeros of $T_{N+1}(x)$.\n",
    "\n",
    "The zeros of $T_N(x)$ in the interval $[-1, 1]$ can be shown to satisfy\n",
    "\n",
    "$$x_k = \\cos\\left( \\frac{(2k - 1) \\pi}{2 N} \\right ) ~~~~~ \\text{for}~~~~ k=0, 1, \\ldots, N-1$$\n",
    "\n",
    "These nodal points (sampling the function at these points) can be shown to minimize interpolation error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(0, numpy.pi, 100)\n",
    "N = 15\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(numpy.cos(x), numpy.sin(x), 'r--')\n",
    "axes.plot(numpy.linspace(-1.1, 1.1, 100), numpy.zeros(x.shape), 'r')\n",
    "for k in xrange(1, N + 1):\n",
    "    location = [numpy.cos((2.0 * k - 1.0) * numpy.pi / (2.0 * N)),\n",
    "                numpy.sin((2.0 * k - 1.0) * numpy.pi / (2.0 * N))]\n",
    "    axes.plot(location[0], location[1], 'ko')\n",
    "    axes.plot(location[0], 0.0, 'ko')\n",
    "    axes.plot([location[0], location[0]], [0.0, location[1]], 'k--')\n",
    "    \n",
    "axes.set_xlim((-1.1, 1.1))\n",
    "axes.set_ylim((-0.1, 1.1))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return 1.0 / (1.0 + 25.0 * x**2)\n",
    "\n",
    "# Parameters\n",
    "x = numpy.linspace(-1, 1, 100)\n",
    "num_points = 6\n",
    "\n",
    "# ============================================================\n",
    "#  Equidistant nodes\n",
    "equidistant_data = numpy.empty((num_points, 2))\n",
    "equidistant_data[:, 0] = numpy.linspace(-1, 1, num_points)\n",
    "equidistant_data[:, 1] = f(equidistant_data[:, 0])\n",
    "N = equidistant_data.shape[0] - 1\n",
    "\n",
    "lagrange_basis = numpy.ones((N + 1, x.shape[0]))\n",
    "for i in xrange(N + 1):\n",
    "    for j in xrange(N + 1):\n",
    "        if i != j:\n",
    "            lagrange_basis[i, :] *= (x - equidistant_data[j, 0]) / (equidistant_data[i, 0] - equidistant_data[j, 0])\n",
    "\n",
    "# Calculate full polynomial\n",
    "P_lagrange = numpy.zeros(x.shape[0])\n",
    "for n in xrange(N + 1):\n",
    "    P_lagrange += lagrange_basis[n, :] * equidistant_data[n, 1]\n",
    "\n",
    "# ============================================================\n",
    "#  Chebyshev nodes\n",
    "chebyshev_data = numpy.empty((num_points, 2))\n",
    "chebyshev_data[:, 0] = numpy.cos((2.0 * numpy.arange(N + 1) - 1.0) * numpy.pi / (2.0 * N))\n",
    "chebyshev_data[:, 1] = f(chebyshev_data[:, 0])\n",
    "\n",
    "# Project data onto Chebyshev polynomials\n",
    "coeff = numpy.polynomial.chebyshev.chebfit(chebyshev_data[:, 0], chebyshev_data[:, 1], N)\n",
    "\n",
    "# Find needed Chebyshev polynomials\n",
    "T = numpy.empty((N, x.shape[0]))\n",
    "T[0, :] = numpy.ones(x.shape)\n",
    "T[1, :] = x\n",
    "for k in xrange(2, N):\n",
    "    T[k, :] = 2.0 * x * T[k-1, :]- T[k-2, :]\n",
    "    \n",
    "# Could plot them here but actually just want to evaluate the polynomial at x\n",
    "P_cheby = numpy.zeros(x.shape)\n",
    "for k in xrange(N):\n",
    "    P_cheby += T[k, :] * coeff[k]\n",
    "\n",
    "# ============================================================\n",
    "# Plot the results\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2.0)\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(x, P_lagrange, 'b', label=\"$P_%s(x)$\" % N)\n",
    "axes.plot(x, f(x), 'k', label=\"True $f(x)$\")\n",
    "axes.plot(equidistant_data[:, 0], equidistant_data[:, 1], 'ro', label=\"data\")\n",
    "axes.set_title(\"Interpolation of Runge's function at Equispaced Points\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"y\")\n",
    "axes.legend(loc=1)\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.plot(x, P_cheby, 'b', label=\"$P_%s(x)$\" % N)\n",
    "axes.plot(x, f(x), 'k', label=\"True $f(x)$\")\n",
    "axes.plot(chebyshev_data[:, 0], chebyshev_data[:, 1], 'ro', label=\"data\")\n",
    "          \n",
    "axes.set_title(\"Interpolation of Runge's function at Chebyshev Points\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"y\")\n",
    "axes.legend(loc=1)\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}