06_interpolation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"./images/CC-BY.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "from __future__ import absolute_import\n",
    "\n",
    "%matplotlib inline\n",
    "import numpy\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Interpolation\n",
    "\n",
    "There are times when you have estimates for the values of a function for specific inputs. The values of the function may be obtained in a variety of ways either through experiment or through the use of other approximation techniques.  Our goal in this chapter is to explore techniques that allow us to determine a new function whose values match the known observations at a set of predetermined input values. We first formally define the term we will use to describe the process."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**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$).\n",
    "\n",
    "**Example 0** The linear polynomial \n",
    "\n",
    "$$\n",
    "    P_1(x) = 2(x-1)+3\n",
    "$$\n",
    "\n",
    "interpolates the coordinates $(1,3)$ and $(3,7)$.\n",
    "\n",
    "In general a polynomial of degree $N$ can be used to interpolate $N+1$ data points. There are many different kinds of functions to use to interpolate values, but here we focus on polynomials."
   ]
  },
  {
   "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",
    "   - The Finite Element Method\n",
    "   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Polynomial Interpolation\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": [
    "We first note that we have two equations and two unknowns. The two equations can be found by assuming the function $P_1(x)$ interpolates the two data points\n",
    "$$\n",
    "    \\begin{align}\n",
    "       y_0 &= p_1 x_0 + p_0, \\\\\n",
    "       y_1 &= p_1 x_1 + p_0.\n",
    "    \\end{align}\n",
    "$$\n",
    "In this example we will solve the first equation for $p_0$, substitute the result into the second equation, and then solve for $p_1$.\n",
    "\n",
    "$$y_0 = p_1 x_0 + p_0 \\quad \\Rightarrow \\quad p_0 = y_0 - p_1 x_0$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\begin{aligned}\n",
    "    y_1 &= p_1 x_1 + p_0 & \\Rightarrow \\\\\n",
    "    y_1 &= p_1 x_1 + y_0 - p_1 x_0 & \\Rightarrow \\\\\n",
    "    p_1 &= \\frac{y_1 - y_0}{x_1 - x_0} & \\Rightarrow \\\\\n",
    "    p_0 &= y_0 - \\frac{y_1 - y_0}{x_1 - x_0} x_0 &\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "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! Note, we have three equations and three unknowns, and the previous system is a linear system of three equations.\n",
    "and in general, the problem will reduce to a linear system\n",
    "$$ \n",
    "    A(\\mathbf{x})\\mathbf{p} = \\mathbf{y}\n",
    "$$\n",
    "A more general approach to solving the system will be explored later, but first it is important to determine whether or not the system even has a solution.    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Proof - Uniqueness of Polynomial Interpolants\n",
    "\n",
    "Let \n",
    "\n",
    "$$\\mathcal{P}_N(x) = \\sum^N_{n=0} p_n x^n $$\n",
    "or\n",
    "$$\\mathcal{P}_N(x) = p_0 + p_1 x + \\cdots + p_{N - 1} x^{N - 1} + p_{N} x^N$$\n",
    "\n",
    "and require $\\mathcal{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": [
    "### Preliminaries: Monomial Basis\n",
    "\n",
    "We can think of  $\\mathcal{P}_N(x) = \\sum^N_{n=0} p_n x^n$ as a polynomial,  or more fundamentally as a *linear combination* of a set of simpler functions, 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": [
    "### Linear independence of the Monomials\n",
    "\n",
    "The monomials, form a *linearly independent* set of functions such that no monomial $x^n$ can be written as a linear combination of any other monomial.  We can see this graphically, for the first few monomials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(-1,1,100)\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "for n in range(4):\n",
    "    axes.plot(x,x**n,label='$x^{}$'.format(n))\n",
    "axes.set_xlabel('x')\n",
    "axes.grid()\n",
    "axes.legend(loc='best')\n",
    "axes.set_title('The First 4 Monomials')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "But more fundamentally. A set of functions is linearly independent if the only linear combination that add to form the zero function, e.g.\n",
    "\n",
    "$$\n",
    "    P_N(x) = p_0 1 + p_1 x + p_2 x^2 + \\ldots + p_n x^n = 0\n",
    "$$\n",
    "\n",
    "is if all the coefficients $p_i = 0$, $\\forall i=0,\\ldots N$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Theorem**:  The monomials $x^0,\\ldots, x^n$ are linear independent.\n",
    "\n",
    "**Proof**: consider $P_N(x) = 0$ for all $x$.  Since the polynomials (and monomials) are differentiable at least $n$ times,  differentiate $n$ times to yield \n",
    "$$\n",
    "    P^{(n)}_N(x) = n!p_n = 0\n",
    "$$ \n",
    "which implies $p_n=0$.  \n",
    "\n",
    "Using  this result and differentiating $n-1$ times shows $p_{n-1}=0$, which by induction gives all $p_i = 0$.\n",
    "\n",
    "Put another way,  the only $n$th degree polynomial that is zero everywhere is if all coefficients are zero."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### The Fundamental theorem of algebra\n",
    "\n",
    "Every $n$th degree polynomial has exactly  $n$ complex roots, i.e. \n",
    "\n",
    "$$\n",
    "    P_N(x) = (x - a_1)(x - a_2)\\ldots(x - a_n)\n",
    "$$\n",
    "for $a_i\\in \\mathbb{C}$.  Therefore, a non-trivial $n$th order polynomial can only be zero at $n$ points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Proof - Uniqueness of Polynomial Interpolants\n",
    "\n",
    "Let \n",
    "\n",
    "$$\\mathcal{P}_N(x) = \\sum^N_{n=0} p_n x^n $$\n",
    "\n",
    "**interpolate** the $N+1$ points $y_i$ at $x_i$.\n",
    "\n",
    "i.e. \n",
    "$$\\mathcal{P}_N(x_i) = y_i,\\quad \\mathrm{for}\\quad i=0,1,\\ldots,N$\n",
    "$$\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) = \\mathcal{P}_N(x) - Q_N(x)$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now, by construction,  $T_N(x_i) = 0$ which implies that it is equal to zero at $n+1$ points.  However, \n",
    "\n",
    "$$T_N(x) = \\mathcal{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",
    "is a $n$th order polynomial which has at most $n$ real roots.  The only way to reconcile this is if T_n(x) = 0, for all $x$, and therefore $p_n - q_n = 0$ individually and therefore $\\mathcal{P}_N(x) = Q_N(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example 3:  Monomial Basis\n",
    "\n",
    "Consider $\\mathcal{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",
    "$$\\mathcal{P}_3(x_0) = p_0 + p_1 x_0 + p_2 x_0^2 + p_3 x_0^3 = y_0$$\n",
    "\n",
    "$$\\mathcal{P}_3(x_1) = p_0 + p_1 x_1 + p_2 x_1^2 + p_3 x_1^3 = y_1$$\n",
    "\n",
    "$$\\mathcal{P}_3(x_2) = p_0 + p_1 x_2 + p_2 x_2^2 + p_3 x_2^3 = y_2$$\n",
    "\n",
    "$$\\mathcal{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",
    "$$\\mathbf{x} = \\begin{bmatrix} x_0 \\\\ x_1 \\\\ x_2 \\\\ x_3 \\end{bmatrix} \\quad \\mathbf{y} = \\begin{bmatrix} y_0 \\\\ y_1 \\\\ y_2 \\\\ y_3 \\end{bmatrix} \\quad \\mathbf{p} = \\begin{bmatrix} p_0 \\\\ p_1 \\\\ p_2 \\\\ p_3 \\end{bmatrix}$$\n",
    "\n",
    "When we write the system in matrix/vector form the matrix that arises is called *Vandermonde* matrix:\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",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can now write the system of linear equations as $V \\mathbf{p} = \\mathbf{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": "fragment"
    }
   },
   "source": [
    "**Note**: the columns  of $V$ are simply the monomial functions sampled at the discrete points $x_i$.  Because the monomials are linearly independent, so are the columns of $V$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\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}$$\n",
    "\n",
    "- What happens if we have redundant data?  Either $(x_i, y_i)$ is repeated or for one $i$ we have two values of $y$.\n",
    "- What if we have more points then the order of polynomial we want?\n",
    "- How does this relate to solving the above linear system of equations?"
   ]
  },
  {
   "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 $\\mathcal{P}_N(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Finding $p_i$\n",
    "\n",
    "Finding the coefficients of $\\mathcal{P}_N(x)$ can be done by solving the system outlined above.  There are functions in `numpy` that can do this for us such as:\n",
    " - `numpy.polyfit(x, y, x.shape[0] - 1)`\n",
    " - `numpy.vander(x, N=None)` to construct the matrix and use a linear solver routine.\n",
    "\n",
    "We can also use a different **basis** that might be easier to use."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Basis\n",
    "\n",
    "**Def:**  A basis for a $N$ dimensional vector space is a set of linearly independent vectors that span the space.\n",
    "\n",
    "The monomials,  $1,x,\\ldots, x^n$, form the usual basis for the vector space of $n$th degree polynomials $P_N(x)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Example** $P_2(x)$ is the space of all quadratic functions. i.e. $P_2(x) = \\mathrm{span}< 1,x,x^2>$\n",
    "\n",
    "$$\n",
    "    P_2(x) = p_0 + p_1 x + p_2 x^2\n",
    "$$\n",
    "\n",
    "i.e for every vector $\\mathbf{p}\\in\\mathbb{R}^3$, there is a unique quadratic function in $P_2(x)$.  (we say $P_2$ is *isomorphic* to $\\mathbb{R}^3$ and is a three dimensional function space)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**However**,  the monomials are not the only basis for $P_N$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Lagrange 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 $\\mathcal{P}_N(x)$ can be written as\n",
    "\n",
    "$$\\mathcal{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_j}{x_i - x_j} = \\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}}$$\n",
    "\n",
    "are the **Lagrange Polynomials**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Lagrange Polynomials\n",
    "\n",
    "$$\\ell_i(x) = \\prod^{N}_{j=0, j \\neq i} \\frac{x - x_j}{x_i - x_j} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "A Key property of the Lagrange polynomials is that \n",
    "\n",
    "$$\n",
    "    \\ell_i(x_j) = \\delta_{ij} = \\left\\{\\begin{matrix}\n",
    "            0 & i\\neq j \\\\ \n",
    "            1 & i=j\\\\\n",
    "            \\end{matrix}\\right.\n",
    "$$\n",
    "which is why the weights in $P_N(x)$ are simply the $y$ values of the interpolant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Solving for the coefficients of $P_N(x)$\n",
    "\n",
    "In general,   if \n",
    "\n",
    "$$\n",
    "    P_N(x) = \\sum_{n=0}^N w_j\\phi_j(x)\n",
    "$$\n",
    "where $\\phi_j(x)$ is any basis function for $P_N$ (i.e. monomial, Lagrange, and there are many more). Then finding the unique set of weights for the interpolating polynomial through $N+1$ distinct data points $(x_i, y_i)$, just reduces to solving $N+1$ linear equations $y_i = P_N(x_i)$.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For the monomial basis this reduces to the linear system\n",
    "\n",
    "$$\n",
    "    V(\\mathbf{x})\\mathbf{w} = \\mathbf{y}\n",
    "$$\n",
    "\n",
    "What is the matrix for the Lagrange Basis?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Visualizing the Lagrange Polynomials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# ====================================================\n",
    "# Compute the Lagrange basis (\\ell_i(x))\n",
    "def lagrange_basis(x, data):\n",
    "    \"\"\"Compute Lagrange basis at x given data\"\"\"\n",
    "    basis = numpy.ones((data.shape[0], x.shape[0]))\n",
    "    for i in range(data.shape[0]):\n",
    "        for j in range(data.shape[0]):\n",
    "            if i != j:\n",
    "                basis[i, :] *= (x - data[j, 0]) / (data[i, 0] - data[j, 0])\n",
    "    return basis\n",
    "\n",
    "# ====================================================\n",
    "# Calculate full polynomial\n",
    "def poly_interpolant(x, data):\n",
    "    \"\"\"Compute polynomial interpolant of (x,y) using Lagrange basis\"\"\"\n",
    "    P = numpy.zeros(x.shape[0])\n",
    "    basis = lagrange_basis(x, data)\n",
    "    for n in range(data.shape[0]):\n",
    "        P += basis[n, :] * data[n, 1]\n",
    "    return P\n",
    "\n",
    "# ===================================================="
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "x_data = numpy.array([0., 1., 2., 3.])\n",
    "y_data = numpy.ones(x_data.shape)\n",
    "data = numpy.array([x_data, y_data]).T\n",
    "\n",
    "x = numpy.linspace(x_data.min(),x_data.max(),100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# ====================================================\n",
    "# Plot individual basis functions\n",
    "fig = plt.figure(figsize=(8, 6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "basis = lagrange_basis(x, data)\n",
    "for i in range(len(x_data)):\n",
    "    axes.plot(x, basis[i, :], label=\"$\\ell_{%s}(x)$\" % i)\n",
    "\n",
    "axes.set_title(\"Lagrange Basis $\\ell_i(x)$\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$\\ell_i(x)$\")\n",
    "axes.grid()\n",
    "axes.legend(loc='best')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Linear Independence of the Lagrange Polynomials\n",
    "\n",
    "Because the weights of each basis function in the Lagrange basis is just the $y$ value at the interpolation points,  it is straightforward to show that the Lagrange polynomials are linearly independent.  I.e. the statement \n",
    "\n",
    "$$\n",
    "    \\sum_{n=0}^N w_j\\phi_j(x) =0\n",
    "$$\n",
    "is equivalent to interpolating the zero function,  where all the $w_j =0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Example 0 Revisited** In example 0 above the linear polynomial that interpolates the coordinates $(1,3)$ and $(3,7)$ was simply stated as \n",
    "$$\n",
    "    P_1(x) = 2(x-1)+3.\n",
    "$$\n",
    "Another way to look at this example is to first note that when we add two linear polynomials \n",
    "the result is another linear polynomial. The first polynomial to define interpolates $(1,1)$ \n",
    "and $(3,0)$,\n",
    "$$\n",
    "    \\ell_0(x) = \\frac{x-3}{1-3}.\n",
    "$$\n",
    "The second polynomial to define interpolates $(1,0)$ and $(3,1)$,\n",
    "$$\n",
    "    \\ell_1(x) = \\frac{x-1}{3-1}.\n",
    "$$\n",
    "\n",
    "A linear combination of these two functions can be defined that will interpolate the points $(1,3)$ and $(3,7)$,\n",
    "$$\n",
    "    P_1(x) = 3\\cdot\\ell_0(x) + 7\\cdot\\ell_1(x).\n",
    "$$\n",
    "The graphs of these functions are shown below.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# =============================================================\n",
    "# Plot the two example basis functions in the current example\n",
    "x = numpy.linspace(1.0, 3.0, 2)\n",
    "\n",
    "fig = plt.figure(figsize=(8, 6))\n",
    "\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.set_ylim([0,9])\n",
    "axes.plot(x, (x-3)/(-2), color='r', label=\"$\\ell_{%s}(x)$\" % 0)\n",
    "axes.plot(x, (x-1)/(2),  color='b', label=\"$\\ell_{%s}(x)$\" % 1)\n",
    "axes.plot(x, 3*(x-3)/(-2) + 7*(x-1)/(2),color='g',label='interpolant')\n",
    "axes.set_title(\"Interpolant for (1,3) and (3,7)\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "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 $\\mathcal{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",
    "$$\\mathcal{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": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "data = numpy.array([[-1.5, -0.5], [0.0, 0.5]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# data = numpy.array([[-1.5, -0.5], [0.0, 0.5], [-0.5, 1.0]])\n",
    "N = data.shape[0] - 1\n",
    "M = data.shape[0]\n",
    "x = numpy.linspace(-2.0, 2.0, 100)\n",
    "\n",
    "\n",
    "# Plot individual basis functions\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "basis = lagrange_basis(x, data)\n",
    "for i in range(N + 1):\n",
    "    axes.plot(x, basis[i, :], label=\"$\\ell_{%s}(x)$\" % i)\n",
    "axes.grid()\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, poly_interpolant(x, data), 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",
    "axes.grid()\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example 5:  Interpolate four points from $sin(\\pi x)$\n",
    "\n",
    "Use four 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": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "num_points = 21"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# num_points = 5\n",
    "# num_points = 6\n",
    "# num_points = 20\n",
    "data = numpy.empty((num_points, 2))\n",
    "data[:, 0] = numpy.linspace(-1, 1, num_points)\n",
    "data[:, 1] = numpy.sin(2.0 * numpy.pi * data[:, 0])\n",
    "N = data.shape[0] - 1 # Degree of polynomial\n",
    "M = data.shape[0]\n",
    "x = numpy.linspace(-1.0, 1.0, 100)\n",
    "\n",
    "# ====================================================\n",
    "# Plot individual basis functions\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "basis = lagrange_basis(x, data)\n",
    "for i in range(N + 1):\n",
    "    axes.plot(x, basis[i, :], label=\"$\\ell_{%s}(x)$\" % i)\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=1)\n",
    "axes.grid()\n",
    "\n",
    "# Plot full polynomial P_N(x)\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.plot(x, poly_interpolant(x, data), label=\"$P_{%s}(x)$\" % N)\n",
    "axes.plot(x, numpy.sin(2.0 * numpy.pi * x), 'r--', label=\"True $f(x)$\")\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",
    "axes.legend(loc=1)\n",
    "axes.grid()\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "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": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return 1.0 / (1.0 + 25.0 * x**2)\n",
    "\n",
    "x = numpy.linspace(-1, 1, 100)\n",
    "# x = numpy.linspace(-2, 2, 100)\n",
    "\n",
    "num_points = 15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# num_points = 10\n",
    "# num_points = 20\n",
    "data = numpy.empty((num_points, 2))\n",
    "data[:, 0] = numpy.linspace(-1, 1, num_points)\n",
    "data[:, 1] = f(data[:, 0])\n",
    "N = data.shape[0] - 1\n",
    "    \n",
    "# Plot the results\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(x, poly_interpolant(x, data), '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",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example 7: Weierstrass \"Monster\" Function\n",
    "\n",
    "Defined as\n",
    "$$\n",
    "    f(x) = \\sum^\\infty_{n=0} a^n \\cos(b^n \\pi x)\n",
    "$$\n",
    "such that\n",
    "$$\n",
    "    0 < a < 1 \\quad \\text{and} \\quad a b > 1 + \\frac{3\\pi}{2}.\n",
    "$$\n",
    "This function is continuous everywhere but not differentiable anywhere."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def f(x, a=0.9, N=100):\n",
    "    summation = 0.0\n",
    "    b = (1.0 + 3.0 / 2.0 * numpy.pi) / a + 0.01\n",
    "    print(b)\n",
    "    for n in range(N + 1):\n",
    "        summation += a**n * numpy.cos(b**n * numpy.pi * x)\n",
    "    return summation\n",
    "\n",
    "x = numpy.linspace(-1, 1, 1000)\n",
    "# x = numpy.linspace(-2, 2, 100)\n",
    "\n",
    "num_points = 10\n",
    "data = numpy.empty((num_points, 2))\n",
    "data[:, 0] = numpy.linspace(-1, 1, num_points)\n",
    "data[:, 1] = f(data[:, 0])\n",
    "N = data.shape[0] - 1\n",
    "    \n",
    "# Plot the results\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(x, poly_interpolant(x, data), '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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "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) = \\mathcal{P}_N(x) + R_N(x)\n",
    "$$\n",
    "where $\\mathcal{P}_N(x)$ is the interpolating polynomial and\n",
    "$$\n",
    "    R_N(x) = Q(x) \\frac{f^{(N+1)}(c)}{(N+1)!} \\quad \\text{with} \\quad c \\in [-1,1]\n",
    "$$\n",
    "with\n",
    "$$\n",
    "    Q(x) = \\prod^N_{i=0} (x - x_i) = (x-x_0)(x-x_1)\\cdots(x-x_N) .\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "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": "subslide"
    }
   },
   "source": [
    "#### Minimizing $R_N(x)$\n",
    "\n",
    "Minimizing the error $R_N(x)$ in Lagrange's theorem is equivalent to minimizing $|Q(x)|$ for $x \\in [-1, 1]$.\n",
    "\n",
    "Minimizing error $\\Leftrightarrow$ picking roots of $Q(x)$ or picking the points where the interpolant data is located.  How do we this?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "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": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def cheb_poly(x, N):\n",
    "    \"\"\"Compute the *N*th Chebyshev polynomial and evaluate it at *x*\"\"\"\n",
    "    T = numpy.empty((3, x.shape[0]))\n",
    "    T[0, :] = numpy.ones(x.shape)\n",
    "    T[1, :] = x\n",
    "    if N == 0:\n",
    "        return T[0, :]\n",
    "    elif N == 1:\n",
    "        return T[1, :]\n",
    "    else:\n",
    "        for k in range(2, N + 1):\n",
    "            T[2, :] = 2.0 * x * T[1, :] - T[0, :]\n",
    "            T[0, :] = T[1, :]\n",
    "            T[1, :] = T[2, :]\n",
    "    \n",
    "        return T[2, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(-1, 1, 100)\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "for n in range(5):\n",
    "    axes.plot(x, cheb_poly(x, n), label=\"$T_%s$\" % n)\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='best')\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "1. Chebyshev nodes of the 1st kind (roots)\n",
    "$$\n",
    "    x_k = \\cos \\left (\\frac{(2 k - 1) \\pi}{2 N} \\right ) \\quad k = 1, \\ldots, N\n",
    "$$\n",
    "1. Chebyshev nodes of the 2nd kind (extrema)\n",
    "$$\n",
    "    x_k = \\cos \\left( \\frac{k \\pi}{N} \\right) \\quad k = 0, \\ldots, N\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "N = 4  \n",
    "x_extrema = numpy.cos(numpy.arange(N + 1) * numpy.pi / N)\n",
    "x_nodes = numpy.cos((2.0 * numpy.arange(1, N + 1) - 1.0) / (2.0 * N) * numpy.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "# fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "# Plot points\n",
    "axes.plot(x_extrema, numpy.zeros(N+1), 'ro')\n",
    "axes.plot(x_nodes, numpy.zeros(N), 'bo')\n",
    "    \n",
    "# Plot some helpful lines\n",
    "axes.plot((-1.0, -1.0), (-1.1, 1.1), 'k--')\n",
    "axes.plot((1.0, 1.0), (-1.1, 1.1), 'k--')\n",
    "axes.plot((-1.0, 1.0), (0.0, 0.0), 'k--')\n",
    "for i in range(x_extrema.shape[0]):\n",
    "    axes.plot((x_extrema[i], x_extrema[i]), (-1.1, 1.1), 'r--')\n",
    "    axes.plot(x_extrema[i], cheb_poly(x_extrema, N)[i], 'ro')\n",
    "\n",
    "print('Nodes   = {}'.format(numpy.sort(x_nodes)))\n",
    "print('Extrema = {}'.format(numpy.sort(x_extrema)))\n",
    "#print(numpy.cos(x_extrema))\n",
    "    \n",
    "# Plot Chebyshev polynomial\n",
    "x_hat = numpy.linspace(-1, 1, 1000)\n",
    "axes.plot(x_hat, cheb_poly(x_hat, N), 'k')\n",
    "axes.set_xlim((-1.1, 1.1))\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "\n",
    "# Labels\n",
    "axes.set_title(\"Chebyshev Nodes and Extrema, N={}\".format(N), fontsize=\"20\")\n",
    "axes.set_xlabel(\"x\", fontsize=\"15\")\n",
    "axes.set_ylabel(\"$T_{N+1}(x)$\", fontsize=\"15\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# First-kind Nesting (3 x)\n",
    "fig = plt.figure()\n",
    "# fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "N = 5\n",
    "factor = 3\n",
    "x_1 = numpy.cos((2.0 * numpy.arange(1, N + 1) - 1.0) / (2.0 * N) * numpy.pi)\n",
    "x_2 = numpy.cos((2.0 * numpy.arange(1, factor * N + 1) - 1.0) / (2.0 * factor * N) * numpy.pi)\n",
    "\n",
    "axes.plot(x_1, numpy.zeros(N), \"o\", color=\"r\", markerfacecolor=\"lightgray\", markersize=\"15\")\n",
    "axes.plot(x_2, numpy.zeros(N * factor), 'kx', markersize=\"10\")\n",
    "\n",
    "x_hat = numpy.linspace(-1, 1, 1000)\n",
    "axes.plot(x_hat, cheb_poly(x_hat, N), 'k')\n",
    "axes.plot(x_hat, cheb_poly(x_hat, factor * N), 'k')\n",
    "axes.set_xlim((-1.1, 1.1))\n",
    "axes.set_ylim((-1.1, 1.1))\n",
    "axes.set_title(\"Nesting of 1st and 2nd Kind Chebyshev Polynomials\")\n",
    "axes.set_xlabel(\"$x$\")\n",
    "axes.set_ylabel(\"$T_N(x)$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Properties of Chebyshev Polynomials\n",
    "1. Defined by a recurrence relation \n",
    "\n",
    "   $$T_k(x) = 2 x T_{k-1}(x) - T_{k-2}(x)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "2. Leading coefficient of $x^N$ in $T_N(x)$ is $2^{N-1}$ for $N \\geq 1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "3. Extreme values:\n",
    "\n",
    "   $$|T_N(x)| \\leq 1 \\quad \\text{for} \\quad  -1 \\leq x \\leq 1$$\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Properties of Chebyshev Polynomials\n",
    "4. Minimax principle:  The polynomial\n",
    "\n",
    "   $$T(x) = \\frac{T_{N+1}(x)}{2^N}$$\n",
    "   \n",
    "   is a *monic polynomial*, a univariate function with the leading coefficient equal to 1, with the property that\n",
    "   \n",
    "   $$\n",
    "       \\max |T(x)| \\leq \\max |Q(X)| \\quad \\text{for} \\quad  x \\in [-1, 1], \\quad \\text{and}\n",
    "   $$\n",
    "   $$\n",
    "       \\max |T(x)| = \\frac{1}{2^N}\n",
    "   $$\n",
    "   \n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    " Recall that the remainder term in the Lagrange Remainder Theorem was\n",
    "   $$\n",
    "        R_N(x) = Q(x) \\frac{f^{(N+1)}(c)}{(N+1)!} \\quad \\text{with} \\quad c \\in [-1,1]\n",
    "   $$\n",
    "   with\n",
    "   $$\n",
    "        Q(x) = \\prod^N_{i=0} (x - x_i) = (x-x_0)(x-x_1)\\cdots(x-x_N) .\n",
    "   $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Error Analysis Redux\n",
    "\n",
    "Given that the Chebyshev polynomials are a minimum on the interval $[-1, 1]$ we would like $T(x) = Q(x)$.  \n",
    "\n",
    "Since we only know the roots of $Q(x)$ (the points where the interpolant data is located) we require these points to be the roots of the Chebyshev polynomial $T_{N+1}(x)$ therefore enforcing $T(x) = Q(x)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "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 ) \\quad \\text{for} \\quad k=1, \\ldots, N\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": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "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, aspect=\"equal\")\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 range(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": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Summary\n",
    "\n",
    "1. Minimizing the error in Lagrange's theorem is equivalent to minimizing \n",
    "$$\n",
    "    |Q(x)| \\quad \\text{for} \\quad x \\in [-1, 1].\n",
    "$$\n",
    "1. We know Chebyshev polynomials are a minimum on the interval $[-1, 1]$ so we would like to have $T(x) = Q(x)$.\n",
    "1. Since we only know the roots of $Q(x)$ (the points where the interpolant data is located) we require these points to be the roots of the Chebyshev polynomial $T_{N+1}(x)$ therefore enforcing $T(x) = Q(x)$.\n",
    "1. 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 ) \\quad \\text{for} \\quad k=1, \\ldots, N\n",
    "$$\n",
    "These nodal points (sampling the function at these points) can be shown to minimize interpolation error."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Notes\n",
    " - The Chebyshev nodes minimize interpolation error for any polynomial basis (due to uniqueness of the interpolating polynomial, any polynomial that interpolates these points are identical regardless of the basis).\n",
    " - Chebyshev nodes uniquely define the Chebyshev polynomials.\n",
    " - The boundedness properties of Chebyshev polynomials are what lead us to the roots as a minimization but there are other used for these orthogonal polynomials.\n",
    " - There are two kinds of Chebyshev nodes and therefore two definitions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Runge's function again \n",
    "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 = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# ============================================================\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",
    "P_lagrange = poly_interpolant(x, equidistant_data)\n",
    "\n",
    "# ============================================================\n",
    "#  Chebyshev nodes\n",
    "chebyshev_data = numpy.empty((num_points, 2))\n",
    "chebyshev_data[:, 0] = numpy.cos((2.0 * numpy.arange(1, num_points + 1) - 1.0) * numpy.pi / (2.0 * num_points))\n",
    "chebyshev_data[:, 1] = f(chebyshev_data[:, 0])\n",
    "\n",
    "P_cheby1 = poly_interpolant(x, chebyshev_data)\n",
    "    \n",
    "# Fit directly with Chebyshev polynomials\n",
    "coeff = numpy.polynomial.chebyshev.chebfit(chebyshev_data[:, 0], chebyshev_data[:, 1], N)\n",
    "P_cheby2 = numpy.polynomial.chebyshev.chebval(x, coeff)\n",
    "\n",
    "# Check on unique polynomials\n",
    "#print(numpy.allclose(P_cheby1, P_cheby2))\n",
    "\n",
    "# calculate errornorms for different interpolants\n",
    "equidistant_err = numpy.linalg.norm(P_lagrange - f(x))\n",
    "cheb_err = numpy.linalg.norm(P_cheby1 - f(x))\n",
    "\n",
    "# ============================================================\n",
    "# Plot the results\n",
    "fig = plt.figure(figsize=(16,6))\n",
    "fig.subplots_adjust(hspace=.5)\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 at Equispaced Points: err = {}\".format(equidistant_err))\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"y\")\n",
    "axes.legend(loc=8)\n",
    "#print('Equispaced error = {}'.format(numpy.linalg.norm(P_lagrange - f(x))))\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.plot(x, f(x), 'k', label=\"True $f(x)$\")\n",
    "axes.plot(x, P_cheby1, 'b', label=\"$P_%s(x)$\" % N)\n",
    "axes.plot(chebyshev_data[:, 0], chebyshev_data[:, 1], 'ro', label=\"data\")\n",
    "          \n",
    "axes.set_title(\"Interpolation at Chebyshev Points: err = {}\".format(cheb_err))\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"y\")\n",
    "axes.legend(loc=1)\n",
    "#print('Chebyshev error = {}'.format(numpy.linalg.norm(P_cheby1 - f(x))))\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Piece-Wise Polynomial Interpolation\n",
    "\n",
    "Given $N$ points, use lower order polynomial interpolation to fit the function in pieces.  We can choose the order of the polynomials and the continuity.\n",
    "\n",
    " - $C^0$: Interpolant is continuous\n",
    "   - Linear interpolation\n",
    "   - Quadratic interpolation\n",
    " - $C^1$: Interpolation and 1st derivative are continuous\n",
    "   - Cubic Hermite polynomials (PCHiP)\n",
    " - $C^2$: Interpolation, 1st and 2nd derivatives are continuous \n",
    "   - Cubic splines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Piece-Wise Linear\n",
    "\n",
    "Given a segment between point $(x_k, y_k)$a nd $(x_{k+1}, y_{k+1})$ define the segment as\n",
    "\n",
    "$$\\mathcal{P}_k(x) = \\frac{y_{k+1} - y_k}{x_{k+1} - x_k} (x - x_k) + y_k$$\n",
    "\n",
    "The final interpolant $\\mathcal{P}(x)$ is then defined on $[x_k, x_{k+1}]$ using this function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "data = numpy.array([[1.0, 3.0], [2.0, 1.0], [3.5, 4.0], [5.0, 0.0], [6.0, 0.5], [9.0, -2.0], [9.5, -3.0]])\n",
    "x = numpy.linspace(0.0, 10, 100)\n",
    "N = data.shape[0] - 1\n",
    "\n",
    "# Lagrange Basis\n",
    "P_lagrange = poly_interpolant(x, data)\n",
    "\n",
    "# C^0 Piece-wise linear\n",
    "# P_pw_linear = numpy.interp(x, data[:, 0], data[:, 1])\n",
    "P_linear = numpy.zeros(x.shape)\n",
    "for n in range(1, N + 1):\n",
    "    P_linear += ((data[n, 1] - data[n - 1, 1]) / (data[n, 0] - data[n - 1, 0]) * (x - data[n - 1, 0]) \n",
    "                         + data[n - 1, 1]) * (x > data[n - 1, 0]) * (x <= data[n, 0])\n",
    "                 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Add end points for continuity\n",
    "P_linear += numpy.ones(x.shape) * data[0, 1] * (x < data[0, 0])\n",
    "P_linear += numpy.ones(x.shape) * data[-1, 1] * (x >= data[-1, 0])\n",
    "\n",
    "# Plot\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(data[:,0], data[:,1], 'ko')\n",
    "axes.plot(x, P_lagrange, 'b--')\n",
    "axes.plot(x, P_linear, 'r')\n",
    "\n",
    "axes.set_title(\"Interpolated Data - $C^0$ Linear\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$P_1(x)$\")\n",
    "axes.set_xlim([0.0, 10.0])\n",
    "axes.set_ylim([-4.0, 15.0])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Piece-Wise Overlapping Polynomials\n",
    "\n",
    "In sets of three points $(x_{k+1}, y_{k+1})$, $(x_{k}, y_{k})$, and $(x_{k-1}, y_{k-1})$, find quadratic interpolant and define final interpolant $P(x)$ using the quadratic interpolant $\\mathcal{P}_k(x)$ on $[x_{k-1}, x_{k+1}]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "data = numpy.array([[1.0, 3.0], [2.0, 1.0], [3.5, 4.0], [5.0, 0.0], [6.0, 0.5], [9.0, -2.0], [9.5, -3.0]])\n",
    "x = numpy.linspace(0.0, 10, 100)\n",
    "N = data.shape[0] - 1\n",
    "\n",
    "# This isn't overlapping, it's more like C_0 P_2\n",
    "# C^0 Piece-wise quadratic\n",
    "P_quadratic = numpy.zeros(x.shape)\n",
    "for k in range(1, N + 1, 2):\n",
    "    p = numpy.polyfit(data[k - 1:k + 2, 0], data[k - 1:k + 2, 1], 2)\n",
    "    P_quadratic += numpy.polyval(p, x) * (x > data[k - 1, 0]) * (x <= data[k + 1, 0])\n",
    "# Add end points for continuity\n",
    "P_quadratic += numpy.ones(x.shape) * data[0, 1] * (x < data[0, 0])\n",
    "P_quadratic += numpy.ones(x.shape) * data[-1, 1] * (x >= data[-1, 0])\n",
    "    \n",
    "# Plot\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(data[:,0], data[:,1], 'ko')\n",
    "axes.plot(x, P_lagrange, 'b--')\n",
    "axes.plot(x, P_quadratic, 'r')\n",
    "\n",
    "axes.set_title(\"Interpolated Data - $C^0$ Quadratic\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$P_3(x)$\")\n",
    "axes.set_xlim([0.0, 10.0])\n",
    "axes.set_ylim([-4.0, 15.0])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Piece-Wise $C^1$ Cubic Interpolation\n",
    "\n",
    "For the previous two cases we had discontinous 1st derivatives!  We can make this better by constraining the polynomials to be continuous at the boundaries of the piece-wise intervals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Given a segment between points $(x_k, y_k)$ and $(x_{k+1}, y_{k+1})$ we want to fit a cubic function between the two points.\n",
    "\n",
    "$$\\mathcal{P}_k(x) = p_0 + p_1 x + p_2 x^2 + p_3 x^3$$\n",
    "\n",
    "$$\\mathcal{P}_k(x_k) = y_k, \\quad \\mathcal{P}_k(x_{k+1}) = y_{k+1}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Now we have 4 unknowns but only two data points!  Constraining the derivative at each interval end will lead to two new equations and therefore we can solve for the interpolant.\n",
    "\n",
    "$$\\frac{\\text{d}}{\\text{dx}} \\mathcal{P}_k(x_k) = d_k, \\quad \\frac{\\text{d}}{\\text{dx}} \\mathcal{P}_k(x_{k+1}) = d_{k+1}$$\n",
    "\n",
    "where we need to prescribe the $d_k$s.  Since we know the polynomial we can write these 4 equations as\n",
    "\n",
    "$$\\begin{aligned}\n",
    "p_0 + p_1 x_k + p_2 x_k^2 + p_3 x_k^3 &= y_k \\\\\n",
    "p_0 + p_1 x_{k+1} + p_2 x_{k+1}^2 + p_3 x_{k+1}^3 &= y_{k+1} \\\\\n",
    "p_1 +  2p_2 x_k + 3 p_3 x_k^2 &= d_k \\\\\n",
    "p_1 + 2 p_2 x_{k+1} + 3 p_3 x_{k+1}^2 &= d_{k+1}\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Rewriting this as a system we get\n",
    "$$\\begin{bmatrix}\n",
    "    1 & x_k & x_k^2 & x_k^3 \\\\\n",
    "    1 & x_{k+1} & x_{k+1}^2 & x_{k+1}^3 \\\\\n",
    "    0 & 1 & 2 x_k & 3 x_k^2 \\\\\n",
    "    0 & 1 & 2 x_{k+1} & 3 x_{k+1}^2\n",
    "\\end{bmatrix} \\begin{bmatrix}\n",
    "p_0 \\\\ p_1 \\\\ p_2 \\\\ p_3\n",
    "\\end{bmatrix} = \\begin{bmatrix}\n",
    "y_k \\\\ y_{k+1} \\\\ d_k \\\\ d_{k+1}\n",
    "\\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A common simplification to the problem description re-parameterizes the locations of the points such that $s \\in [0, 1]$ and recast the problem with $(0, y_k)$ and $(1, y_{k+1})$.  This simplifies the above system to\n",
    "$$\\begin{bmatrix}\n",
    "    1 & 0 & 0 & 0 \\\\\n",
    "    1 & 1 & 1 & 1 \\\\\n",
    "    0 & 1 & 0 & 0 \\\\\n",
    "    0 & 1 & 2 & 3\n",
    "\\end{bmatrix} \\begin{bmatrix}\n",
    "p_0 \\\\ p_1 \\\\ p_2 \\\\ p_3\n",
    "\\end{bmatrix} = \\begin{bmatrix}\n",
    "y_k \\\\ y_{k+1} \\\\ d_k \\\\ d_{k+1}\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "which can be solved to find\n",
    "\n",
    "$$\\begin{aligned}\n",
    "\\mathcal{P}(s) &= (1-s)^2 (1 + 2s) y_k + s^2 (3 - 2 s) y_{k+1} + s (1 - s)^2 d_k - s^2 (1 - s)d_{k+1}\\\\\n",
    "\\mathcal{P}'(s) &= 6s(s-1) y_k + 6s(1-s) y_{k+1} + (s-1)(3s-1) d_k - s(3s-2) d_{k+1}\\\\\n",
    "\\mathcal{P}''(s) &= 6 (1-2s)(y_{k+1} - y_k) + (6s - 4) d_k + (6s-2) d_{k+1}\n",
    "\\end{aligned}$$\n",
    "\n",
    "Now, how to choose $d_k$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### PCHIP\n",
    "\n",
    "Piecewise Cubic Hermite Interpolation Polynomial\n",
    "\n",
    " - Picks the slope that preserves monotonicity\n",
    " - Also tried to preserve the shape of the data\n",
    " - Note that in general this interpolant is $\\mathcal{P}_k(x) \\in C^1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import scipy.interpolate as interpolate\n",
    "\n",
    "data = numpy.array([[1.0, 3.0], [2.0, 1.0], [3.5, 4.0], [5.0, 0.0], [6.0, 0.5], [9.0, -2.0], [9.5, -3.0]])\n",
    "x = numpy.linspace(0.0, 10, 100)\n",
    "\n",
    "# C^1 Piece-wise PCHIP\n",
    "P_pchip = interpolate.pchip_interpolate(data[:, 0], data[:, 1], x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Plot\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(data[:,0], data[:,1], 'ro')\n",
    "axes.plot(x, P_pchip, 'r')\n",
    "\n",
    "axes.set_title(\"Interpolated Data - $C^1$ Cubic PCHIP\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$P_3(x)$\")\n",
    "axes.set_xlim([0.0, 10.0])\n",
    "axes.set_ylim([-4.0, 15.0])\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Cubic Splines\n",
    "\n",
    "Enfores continuity on second derivatives as well:\n",
    "\n",
    "$$\\mathcal{P}''_{k}(x_{k}) = \\mathcal{P}''_{k-1}(x_k)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "From our generalization before we know\n",
    "\n",
    "$$\\mathcal{P}''(s) = 6 (1-2s)(y_{k+1} - y_k) + (6s - 4) d_k + (6s-2) d_{k+1}$$\n",
    "\n",
    "and our constraint now becomes\n",
    "\n",
    "$$\\mathcal{P}''_{k}(0) = \\mathcal{P}''_{k-1}(1)$$\n",
    "\n",
    "$$\\mathcal{P}''_{k-1}(1) = 6 (1-2 \\cdot 1)(y_{k} - y_{k-1}) + (6\\cdot 1 - 4) d_{k-1} + (6\\cdot 1-2) d_{k}$$\n",
    "\n",
    "$$\\mathcal{P}''_{k}(0) = 6 (1-2 \\cdot 0)(y_{k+1} - y_k) + (6\\cdot 0 - 4) d_k + (6\\cdot 0-2) d_{k+1}$$\n",
    "\n",
    "$$-6(y_{k} - y_{k-1}) + 2 d_{k-1} + 4 d_{k} = 6 (y_{k+1} - y_k) - 4 d_k -2  d_{k+1}$$\n",
    "\n",
    "We now have constraints on choosing the $d_k$ values.  Note that we still need to prescribe them at the boundaries of the full interval."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This forms a linear set of equations for the $d_k$s based on the $y_k$ values and can be reformulated into a tri-diagonal linear system\n",
    "\n",
    "$$\\begin{bmatrix}\n",
    "         & \\ddots  & \\ddots  & \\ddots  & &  &\\\\\n",
    "         & 0 & 2 & 8 & 2 & 0 & & \\\\\n",
    "         & & 0 & 2 & 8 & 2 & 0 & & & \\\\\n",
    "         & &   & 0 & 2 & 8 & 2 & 0 & & \\\\\n",
    "         & &    & & & \\ddots & \\ddots & \\ddots &\n",
    "\\end{bmatrix}\\begin{bmatrix}\n",
    "    \\vdots \\\\ d_{k-1} \\\\ d_{k} \\\\ d_{k+1} \\\\ \\vdots\n",
    "\\end{bmatrix} = \\begin{bmatrix}\n",
    "    \\vdots \\\\ 6 (y_{k} - y_{k-2}) \\\\ 6 (y_{k+1} - y_{k-1}) \\\\ 6 (y_{k+2} - y_{k}) \\\\\\vdots\n",
    "\\end{bmatrix}$$\n",
    "\n",
    "The boundaries are still left unconstrained and we must pick some rule to specify the derivatives there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "import scipy.interpolate as interpolate\n",
    "\n",
    "data = numpy.array([[1.0, 3.0], [2.0, 1.0], [3.5, 4.0], [5.0, 0.0], [6.0, 0.5], [9.0, -2.0], [9.5, -3.0]])\n",
    "x = numpy.linspace(0.0, 10, 100)\n",
    "\n",
    "# C^2 Piece-wise Splines\n",
    "# Note that to get an interpolant we need to set the smoothing \n",
    "# parameters *s* to 0\n",
    "P_spline = interpolate.UnivariateSpline(data[:, 0], data[:, 1], s=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# Plot\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(data[:,0], data[:,1], 'ro')\n",
    "axes.plot(x, P_spline(x), 'r', label = '$C^2$')\n",
    "axes.plot(x, P_pchip, 'b--', label = 'Pchip')\n",
    "axes.set_title(\"Interpolated Data - $C^2$ Cubic Splines\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$P_3(x)$\")\n",
    "axes.set_xlim([0.0, 10.0])\n",
    "axes.set_ylim([-4.0, 15.0])\n",
    "axes.grid()\n",
    "axes.legend(loc='best')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Let's compare all of these methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import scipy.interpolate as interpolate\n",
    "\n",
    "data = numpy.array([[1.0, 3.0], [2.0, 1.0], [3.5, 4.0], [5.0, 0.0], [6.0, 0.5], [9.0, -2.0], [9.5, -3.0]])\n",
    "x = numpy.linspace(0.0, 10, 100)\n",
    "\n",
    "# Lagrange Basis\n",
    "N = data.shape[0] - 1\n",
    "\n",
    "lagrange_basis = numpy.ones((N + 1, x.shape[0]))\n",
    "for i in range(N + 1):\n",
    "    for j in range(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_lagrange = numpy.zeros(x.shape[0])\n",
    "for n in range(N + 1):\n",
    "    P_lagrange += lagrange_basis[n, :] * data[n, 1]\n",
    "\n",
    "# C^0 Piece-wise linear\n",
    "# P_pw_linear = numpy.interp(x, data[:, 0], data[:, 1])\n",
    "P_linear = numpy.zeros(x.shape)\n",
    "for n in range(1, N + 1):\n",
    "    P_linear += ((data[n, 1] - data[n - 1, 1]) / (data[n, 0] - data[n - 1, 0]) * (x - data[n - 1, 0]) \n",
    "                         + data[n - 1, 1]) * (x > data[n - 1, 0]) * (x <= data[n, 0])\n",
    "    \n",
    "# Add end points for continuity\n",
    "P_linear += numpy.ones(x.shape) * data[0, 1] * (x < data[0, 0])\n",
    "P_linear += numpy.ones(x.shape) * data[-1, 1] * (x >= data[-1, 0])\n",
    "# C^0 Piece-wise quadratic\n",
    "P_quadratic = numpy.zeros(x.shape)\n",
    "for k in range(1, N + 1, 2):\n",
    "    p = numpy.polyfit(data[k - 1:k + 2, 0], data[k - 1:k + 2, 1], 2)\n",
    "    P_quadratic += numpy.polyval(p, x) * (x > data[k - 1, 0]) * (x <= data[k + 1, 0])\n",
    "# Add end points for continuity\n",
    "P_quadratic += numpy.ones(x.shape) * data[0, 1] * (x < data[0, 0])\n",
    "P_quadratic += numpy.ones(x.shape) * data[-1, 1] * (x >= data[-1, 0])\n",
    "    \n",
    "# C^1 Piece-wise PCHIP\n",
    "P_pchip = interpolate.pchip_interpolate(data[:, 0], data[:, 1], x)\n",
    "\n",
    "# C^2 Piece-wise Splines\n",
    "P_spline = interpolate.UnivariateSpline(data[:, 0], data[:, 1], s=0)\n",
    "    \n",
    "# Plot\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(data[:,0], data[:,1], 'ko', label=\"Data\")\n",
    "axes.plot(x, P_lagrange, 'y', label=\"Lagrange\")\n",
    "axes.plot(x, P_linear, 'g', label=\"PW Linear\") \n",
    "axes.plot(x, P_quadratic, 'r', label=\"PW Quadratic\")\n",
    "axes.plot(x, P_pchip, 'c', label=\"PW Cubic - PCHIP\")\n",
    "axes.plot(x, P_spline(x), 'b', label=\"PW Cubic - Spline\")\n",
    "axes.grid()\n",
    "\n",
    "axes.set_title(\"Interpolated Data - Method Comparisons\")\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$P(x)$\")\n",
    "axes.legend(loc='best')\n",
    "axes.set_xlim([0.0, 10.0])\n",
    "axes.set_ylim([-4.0, 15.0])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Relationship to Regression\n",
    "\n",
    "What if we have more data and want a lower degree polynomial but do not want to use a piece-wise defined interpolant?  \n",
    "Regression techniques are often used to minimize a form of error between the data points $y_i$ at $x_i$ with an approximating function $f(x_i)$.  Note that this is NOT interpolation anymore!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Least-Squares\n",
    "\n",
    "One way of doing this is to require that we minimize the least-squares error\n",
    "$$\n",
    "    E = \\left( \\sum^m_{i=1} |y_i - f(x_i)|^2 \\right )^{1/2}.\n",
    "$$\n",
    "where as before we have data $y_i$ at locations $x_i$ and an approximating function $f(x_i)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "From the beginning of our discussion we know we can write the interpolant as a system of linear equations which we can then solve for the coefficients of a monomial basis.  If we wanted to fit a line \n",
    "$$\n",
    "    \\mathcal{P}_1(x) = p_0 + p_1 x\n",
    "$$\n",
    "to $N$ data points we would have\n",
    "$$\n",
    "   \\begin{bmatrix}\n",
    "        1 & x_1 \\\\\n",
    "        1 & x_2 \\\\\n",
    "        \\vdots & \\vdots \\\\\n",
    "        1 & x_N\n",
    "   \\end{bmatrix} \\begin{bmatrix}\n",
    "       p_0 \\\\ p_1\n",
    "   \\end{bmatrix} = \\begin{bmatrix}\n",
    "       y_1 \\\\ y_2 \\\\ \\vdots \\\\ y_N\n",
    "   \\end{bmatrix}\n",
    "$$\n",
    "or\n",
    "$$\n",
    "    A p = y\n",
    "$$\n",
    "What's wrong with this system?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "This leads to the likelihood that there is no solution to the system as\n",
    "$$\n",
    "    A \\in \\mathbb{R}^{N \\times 2}, p \\in \\mathbb{R}^{2 \\times 1}, \\text{  and  } y \\in \\mathbb{R}^{N \\times 1}.\n",
    "$$\n",
    "\n",
    "Instead we can solve the related least-squares system\n",
    "$$\n",
    "    A^T A p = A^T y\n",
    "$$\n",
    "whose solution minimizes the least-square error defined before as $E$.\n",
    "\n",
    "Note:  In general, this is not the most stable way to solve least squares problems,  in general, using an orthogonalization technique like $QR$ factorization is better numerically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Linear Least Squares Problem\n",
    "N = 50\n",
    "x = numpy.linspace(-1.0, 1.0, N)\n",
    "y = x + numpy.random.random((N)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "A = numpy.ones((x.shape[0], 2))\n",
    "A[:, 1] = x\n",
    "p = numpy.linalg.solve(numpy.dot(A.transpose(), A), numpy.dot(A.transpose(), y))\n",
    "#p = numpy.linalg.lstsq(A, y, rcond=None)[0]\n",
    "f = lambda x: p[0] + p[1] * x\n",
    "E = numpy.linalg.norm(y - f(x), ord=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(x, y, 'ko')\n",
    "axes.plot(x, f(x), 'r')\n",
    "axes.set_title(\"Least Squares Fit to Data, err={}\".format(E))\n",
    "axes.set_xlabel(\"$x$\")\n",
    "axes.set_ylabel(\"$f(x)$ and $y_i$\")\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Themes and variations\n",
    "\n",
    "You can play all sorts of games, whether they are justified by the data or not,  for example we can fit the same random data with a function like\n",
    "\n",
    "$$\n",
    "    f(x) = p_0 + p_1\\tanh(x)\n",
    "$$\n",
    "\n",
    "which is still a linear problem for the coefficients $p_0$ and $p_1$, however the vandermonde matrix now has columns of $\\mathbf{1}$ and $\\tanh\\mathbf{x}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Linear Least Squares Problem\n",
    "A = numpy.ones((x.shape[0], 2))\n",
    "A[:, 1] = numpy.tanh(x)\n",
    "p = numpy.linalg.solve(numpy.dot(A.transpose(), A), numpy.dot(A.transpose(), y))\n",
    "# p = numpy.linalg.lstsq(A, y)[0]\n",
    "f = lambda x: p[0] + p[1] * numpy.tanh(x)\n",
    "E = numpy.linalg.norm(y - f(x), ord=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(x, y, 'ko')\n",
    "axes.plot(x, f(x), 'r')\n",
    "axes.set_title(\"Least Squares Fit to Data, err = {}\".format(E))\n",
    "axes.set_xlabel(\"$x$\")\n",
    "axes.set_ylabel(\"$f(x)$ and $y_i$\")\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Let ye be warned...\n",
    "\n",
    "![curve-fitting methods](https://imgs.xkcd.com/comics/curve_fitting.png)\n",
    "\n",
    "(Original image can be found at [Curve Fitting](https://xkcd.com/2048/).)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.5"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}