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",
    "Interpolation is a fundamental numerical problem that is central to many other numerical algorithms.  Simply put, given a finite number of points where values are known,  find an *interpolant*, a (usually) continuous function that returns values everywhere and is guaranteed to pass through the known data.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Objectives: \n",
    "\n",
    "* Define an interpolant\n",
    "* Understand **Polynomial** Interpolation and Approximation\n",
    "    * Define the interpolating polynomial $P_N(x)$ that exactly interpolates $N+1$ points\n",
    "    * Calculation (in the monomial basis)\n",
    "    * Uniqueness\n",
    "    * Other bases (Lagrange, Newton?)\n",
    "    * Error Analysis\n",
    "    * Chebyshev Polynomials and optimal interpolation\n",
    "* Understand Other interpolants\n",
    "    * Piecewise Polynomial interpolation\n",
    "    * Overlapping Polynomial interpolation\n",
    "    * Cubic Splines and other smooth interpolants\n",
    "    * Higher dimensional Interpolation schemes and scipy.interpolate\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Interpolation (vs  Fitting)\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$).\n",
    "\n",
    "A visual example for 3 random points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.interpolate import pchip_interpolate, make_interp_spline\n",
    "\n",
    "N= 3\n",
    "x = [ 1, 2.3, 4.8 ]\n",
    "y = [ 1., 0.5, 5.0]\n",
    "p2 = numpy.polyfit(x,y,2)\n",
    "p1 = numpy.polyfit(x,y,1)\n",
    "\n",
    "b_spline = make_interp_spline(x, y, bc_type=\"natural\") \n",
    "\n",
    "xa = numpy.linspace(0.,5)\n",
    "xs = numpy.linspace(x[0],x[-1])\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(x,y,'ro',markersize=10,label='data')\n",
    "axes.plot(xa,numpy.polyval(p2, xa), label='quadratic interpolant')\n",
    "axes.plot(x,y,'g', label='Piecewise linear interpolant')\n",
    "axes.plot(xs,pchip_interpolate(x,y,xs),'k', label='Piecewise cubic interpolant')\n",
    "axes.plot(xs,b_spline(xs),'c', label='B-spline')\n",
    "\n",
    "\n",
    "axes.plot(xa,numpy.polyval(p1, xa), label='best fit line')\n",
    "axes.legend(loc='best')\n",
    "axes.set_xlim(min(x) - 0.5,max(x) + 0.5)\n",
    "axes.set_xlabel('x',fontsize=16)\n",
    "axes.set_ylabel('f(x)',fontsize=16)\n",
    "\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Comment**:\n",
    "In general a polynomial of degree $N$ can be used to interpolate $N+1$ data points. There are many choices 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_0, y_1, \\ldots, y_N $ at *distinct* points $x_0, x_1, \\ldots, x_N$.\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 two points $(x_0, y_0)$ and $(x_1, y_1)$, There is a unique line \n",
    "\n",
    "$$P_1(x) = p_0 + p_1 x$$\n",
    "\n",
    "that connects them. We simply need to use the data to find $p_0$ and $p_1$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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",
    "$$\n",
    "    \\begin{matrix}\n",
    "       P_1(x_0) = & p_0 + p_1 x_0  =& y_0  \\\\\n",
    "       P_1(x_1) =& p_0 + p_1 x_1 = &   y_1 \\\\\n",
    "    \\end{matrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We can also (and should) write this problem as a small $2\\times2$ linear algebra problem $A\\mathbf{x}=\\mathbf{b}$\n",
    "\n",
    "$$\n",
    "    \\begin{bmatrix}\n",
    "       1 & x_0 \\\\\n",
    "       1 & x_1 \\\\\n",
    "    \\end{bmatrix}\n",
    "    \\begin{bmatrix}\n",
    "       p_0 \\\\\n",
    "       p_1 \\\\\n",
    "    \\end{bmatrix}=\\begin{bmatrix}\n",
    "       y_0 \\\\\n",
    "       y_1 \\\\\n",
    "    \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Question:**  What are the unknowns,  and where does the data sit in $A$ and $\\mathbf{b}$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "With a bit of algebra you should be able to show that the solution of this problem is\n",
    "\n",
    "$$\n",
    "    p_1 = \\frac{y_1 - y_0}{x_1 - x_0}, \\quad\\quad\n",
    "    p_0 = y_0 - \\frac{y_1 - y_0}{x_1 - x_0} x_0 \n",
    "$$\n",
    "    \n",
    "and \n",
    "\n",
    "$$P_1(x)  = \\frac{y_1 - y_0}{x_1 - x_0} (x - x_0) + y_0$$\n",
    "\n",
    "which is just the equation of the straight line with slope $p_1$, that connects the two points.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Example 2:  3 Points\n",
    "\n",
    "Given three points  $(x_0, y_0)$, $(x_1, y_1)$, and $(x_2, y_2)$  leads to a quadratic polynomial:\n",
    "$$P_2(x) = p_0 + p_1 x + p_2x^2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Again,  setting $P_2(x_i) = y_i$  yields 3 equations in 3 unknowns\n",
    "\n",
    "$$ p_0 + p_1 x_0 + p_2 x_0^2  = y_0$$\n",
    "\n",
    "$$ p_0 + p_1 x_1 + p_2 x_1^2  = y_1$$\n",
    "\n",
    "$$ p_0 + p_1 x_2 + p_2 x_2^2  = y_2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Which reduces to the $3\\times3$ linear system\n",
    "$ \n",
    "    A(\\mathbf{x})\\mathbf{p} = \\mathbf{y}\n",
    "$\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",
    "\n",
    "for $a_i\\in \\mathbb{C}$.  Therefore, a _non-trivial_ $n$th order polynomial can only be zero at at most $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",
    "$$\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: First define the following vectors\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 a *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 \n",
    "\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$$\n",
    "\n",
    "as $V(\\mathbf{x}) \\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$. i.e $V_j = \\mathbf{x}^j$.\n",
    "\n",
    "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 also 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": [
    "### Visualizing the Lagrange Polynomials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "# ====================================================\n",
    "# Compute the Lagrange basis (\\ell_i(x))\n",
    "def lagrange_basis(x, data):\n",
    "    \"\"\"Compute Lagrange basis at x given N data points\n",
    "    \n",
    "    params:\n",
    "    -------\n",
    "        x: ndarray\n",
    "            1-d Array of floats\n",
    "        data: ndarray of shape (N,2)\n",
    "            2-d Array of data where each row is [ x_i, y_i ]\n",
    "    returns:\n",
    "    --------\n",
    "        basis: ndarray of shape (N, x.shape)\n",
    "            : 2-D array of lagrange basis functions evaluated at x\n",
    "            \n",
    "    \"\"\"\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# ====================================================\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": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x_data = numpy.array([0., 1., 2.])\n",
    "#x_data = numpy.array([0., 1., 2., 3.])\n",
    "y_data = numpy.ones(x_data.shape)\n",
    "\n",
    "data = numpy.array([x_data, y_data]).T\n",
    "x = numpy.linspace(x_data.min(),x_data.max(),100)\n",
    "basis = lagrange_basis(x, data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# ====================================================\n",
    "# Plot individual basis functions\n",
    "\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",
    "axes.plot(x_data,numpy.zeros(x_data.shape),'ko')\n",
    "axes.plot(x_data,y_data,'k+',markersize=10)\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": [
    "### Solving for the coefficients of $P_N(x)$\n",
    "\n",
    "In general,   if \n",
    "\n",
    "$$\n",
    "    P_N(x) = \\sum_{j=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 $P_N(x_i) = y_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",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "What is the matrix $V$, for the Lagrange Basis?  What are the weights $\\mathbf{w}$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "A little pen-and-paper work here can be instructive$\\ldots$"
   ]
  },
  {
   "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\\ell_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 :  $N = 1$ 1st order general Lagrange Polynomial as linear interpolant\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",
    "and the 1st order interpolating Polynomial is simply\n",
    "\n",
    "$$\\mathcal{P}_1(x) = y_0\\ell_0(x)+ y_1\\ell_1(x)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The behavior of $\\mathcal{P}_1(x)$ becomes clearer if we note that on the interval $x\\in[x_0,x_1]$ that \n",
    "\n",
    "$$\n",
    "\\ell_1(x) = \\frac{x- x_0}{x_1 - x_0} = s\n",
    "$$ \n",
    "is simply the fractional distance across the interval.\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We should also note that \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    \\ell_0(x) &= \\frac{x - x_1}{x_0 - x_1}\\\\\n",
    "    &= 1-\\ell_1(x) \\\\\n",
    "    &= 1-s\\\\\n",
    "    \\end{align}\n",
    "$$ \n",
    "\n",
    "(show this)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Then the interpolating polynomial is simply\n",
    "$$\n",
    "\\begin{align}\n",
    "\\mathcal{P}_1(s) &= y_0(1-s) + y_1 s \\\\&= y_0 + s(y_1 - y_0)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "for $s\\in[0,1]$ which is just the linear line segment that connects points $y_0$ and $y_1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "As a specific example we will plot $\\ell_0(x)$, $\\ell_1(x)$ and $\\mathcal{P}_1$ for the interval \n",
    "$x=[1,3]$  and $y_0=2$, $y_1=3$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x0, y0 = (1., 2)\n",
    "x1, y1 = (3., 3.)\n",
    "ell0 = lambda x: (x - x1)/(x0 - x1)\n",
    "ell1 = lambda x: (x - x0)/(x1 - x0)\n",
    "P1 = lambda x:  y0*ell0(x) + y1*ell1(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "x = numpy.linspace(x0,x1)\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(x,ell0(x),'r',label='$\\ell_0$')\n",
    "axes.plot(x,ell1(x),'b',label='$\\ell_1$')\n",
    "axes.plot(x,P1(x),'g',label='$P_1(x)$')\n",
    "axes.plot((x0,x1),(y0,y1),'go')\n",
    "axes.set_xlabel('x')\n",
    "axes.grid()\n",
    "axes.legend(loc='best',fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Example: Interpolate six points from $\\sin(2\\pi x)$\n",
    "\n",
    "Use six evenly spaced points to approximate $\\sin$ on the interval $x \\in [-1, 1]$.  What is the behavior as $N \\rightarrow \\infty$?  Also plot the absolute 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 = 6"
   ]
  },
  {
   "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",
    "px = poly_interpolant(x, data)\n",
    "f = numpy.sin(2.0 * numpy.pi * x)\n",
    "err = numpy.abs(f - px)\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)$\",fontsize=16)\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$\\ell_i(x)$\")\n",
    "axes.legend(loc='best',fontsize=14)\n",
    "axes.grid()\n",
    "\n",
    "# Plot full polynomial P_N(x)\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.plot(x, px , label=\"$P_{%s}(x)$\" % N)\n",
    "axes.plot(x, numpy.sin(2.0 * numpy.pi * x), 'r--', label=\"True $f(x)$\")\n",
    "axes.plot(x, err, 'k', linestyle='dotted',label='abs error')\n",
    "for point in data:\n",
    "    axes.plot(point[0], point[1], 'ko')\n",
    "axes.set_title(\"$P_N(x)$\",fontsize=16)\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"$P_N(x)$\")\n",
    "axes.legend(loc='best',fontsize=14)\n",
    "axes.grid()\n",
    "    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "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",
    "num_points = 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "\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_{{{name}}}(x)$\".format(name=num_points-1))\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\",fontsize=18)\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"y\")\n",
    "axes.legend(loc='best',fontsize=14)\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "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": "skip"
    }
   },
   "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",
    "$$\n",
    "    f(x) = \\mathcal{P}_N(x) + R_N(x)\n",
    "$$\n",
    "\n",
    "where $\\mathcal{P}_N(x)$ is the interpolating polynomial and\n",
    "\n",
    "$$\n",
    "    R_N(x) = Q(x) \\frac{f^{(N+1)}(c)}{(N+1)!} \\quad \\text{with} \\quad c \\in [-1,1]\n",
    "$$\n",
    "where\n",
    "$$\n",
    "    Q(x) = \\prod^N_{i=0} (x - x_i) = (x-x_0)(x-x_1)\\cdots(x-x_N) .\n",
    "$$\n",
    "\n",
    "is a *monic* polynomial of order $N+1$ with exactly $N+1$ roots at the nodes $x_i$"
   ]
  },
  {
   "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$$\n",
    "\n",
    "In general,  the Chebyshev polynomials are generated by a recurrence relation \n",
    "\n",
    "$$T_k(x) = 2 x T_{k-1}(x) - T_{k-2}(x)$$"
   ]
  },
  {
   "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": [
    "#### Chebyshev nodes\n",
    "\n",
    "Chebyshev polynomials have many special properties and locations including the location of their roots and extrema known as Chebyshev nodes\n",
    "\n",
    "\n",
    "* Chebyshev nodes of the 1st kind (roots)\n",
    "\n",
    "$$\n",
    "    x_k = \\cos \\left (\\frac{(2 k - 1) \\pi}{2 N} \\right ) \\quad k = 1, \\ldots, N\n",
    "$$\n",
    "\n",
    "* Chebyshev nodes of the 2nd kind (extrema)\n",
    "\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": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "N = 6  \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",
    "   $$ T_0(x) = 1,\\quad T_1(x) = x$$\n",
    "   $$T_k(x) = 2 x T_{k-1}(x) - T_{k-2}(x),\\quad k\\geq2$$"
   ]
  },
  {
   "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 $Q(x) = T(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 $Q(x)=T(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 uses for these orthogonal polynomials.\n",
    " - There are two kinds of Chebyshev nodes and therefore two definitions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example:  Chebyshev Interpolation of Runge's function\n",
    "\n",
    "Comparison between interpolation at Chebyshev Nodes vs equally spaced points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "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 = 6"
   ]
  },
  {
   "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_{{ {} }}(x)$'.format(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(\"Equispaced Points: err = {:5.5g}\".format(equidistant_err),fontsize=18)\n",
    "axes.set_xlabel(\"x\",fontsize=16)\n",
    "axes.set_ylabel(\"y\",fontsize=16)\n",
    "axes.grid()\n",
    "axes.legend(loc=8,fontsize=16)\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_{{ {} }}(x)$'.format(N))\n",
    "axes.plot(chebyshev_data[:, 0], chebyshev_data[:, 1], 'ro', label=\"data\")\n",
    "          \n",
    "axes.set_title(\"Chebyshev Points: err = {:5.5g}\".format(cheb_err),fontsize=18)\n",
    "axes.set_xlabel(\"x\",fontsize=16)\n",
    "axes.set_ylabel(\"y\",fontsize=16)\n",
    "axes.legend(loc=1,fontsize=16)\n",
    "axes.grid()\n",
    "#print('Chebyshev error = {}'.format(numpy.linalg.norm(P_cheby1 - f(x))))\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Affine Transformation of intervals\n",
    "\n",
    "While the chebyshev polynomials (and their roots) are defined on the interval \n",
    "$$\n",
    "    s \\in[-1, 1]\n",
    "$$ \n",
    "\n",
    "it turns out we can use them to approximate functions on any arbitrary interval $x\\in[a,b]$ as long as we can map the interval $[-1,1]$ onto $[a,b]$.   This is easy to do given the  *affine transformation*\n",
    "\n",
    "$$\n",
    "    x(s) = \\alpha + \\beta s\n",
    "$$\n",
    "\n",
    "defined by a scaling $\\beta$ and a shift $\\alpha$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "To solve for $\\alpha$ and $\\beta$, we simply transform the two endpoints such that\n",
    "$$\n",
    "\\begin{align}\n",
    "    x(-1) &= \\alpha - \\beta  = a \\\\\n",
    "    x(1) &= \\alpha + \\beta  = b \\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or\n",
    "\n",
    "$$\n",
    "    \\begin{bmatrix} 1 & -1 \\\\ 1 & 1\\end{bmatrix}\\begin{bmatrix} \\alpha \\\\ \\beta \\end{bmatrix}=\\begin{bmatrix} a \\\\ b \\end{bmatrix}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "or solving\n",
    "\n",
    "$$\n",
    "    \\begin{align}\n",
    "        \\alpha &= \\frac{a+b}{2} \\\\\n",
    "        \\beta &= \\frac{b -a}{2}\\\\\n",
    "        \\end{align}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Affine Transformation of intervals\n",
    "\n",
    "Thus if we know the Chebyshev nodes for $s\\in[-1, 1]$\n",
    "$$\n",
    "    s_k = \\cos\\left( \\frac{(2k - 1) \\pi}{2 N} \\right ) \\quad \\text{for} \\quad k=1, \\ldots, N\n",
    "$$ \n",
    "\n",
    "then the correct interpolation points on $[a,b]$ are simply $x_k = x(s_k)$\n",
    "\n",
    "For example:  if $s_k=0$ (in the center of the interval), then \n",
    "\n",
    "$$x_k = x(0) = \\frac{a+b}{2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Note\n",
    "\n",
    "Affine Transformation of the interval is useful for may problems and applications."
   ]
  },
  {
   "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)$and $(x_{k+1}, y_{k+1})$ we can do linear interpolation for all $x\\in[x_k, x_{k+1}]$ as \n",
    "\n",
    "$$\n",
    "    \\mathcal{P}_k(x) = y_k(1 - s) + y_{k+1}s\n",
    "$$\n",
    "\n",
    "where \n",
    "$$\n",
    "    s = \\frac{x - x_k}{x_{k+1} - x_k}\n",
    "$$\n",
    "\n",
    "is the fractional distance of $x$ within the segment (i.e. $s=\\ell_1(x)$ for $x\\in[x_k, x_{k+1}]$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def piece_wise_linear_interp(x, data):\n",
    "    \"\"\" return the piecewise linear interpolant of x given an array of data points\"\"\"\n",
    "    x_data = data[:,0]\n",
    "    y_data = data[:,1]\n",
    "    \n",
    "    y = numpy.zeros(x.shape)\n",
    "    \n",
    "    # loop over segments\n",
    "    n_segments = len(data) - 1\n",
    "    for k in range(n_segments):\n",
    "        # define endpoints of the segment k\n",
    "        xk  = x_data[k]\n",
    "        xkp = x_data[k+1]\n",
    "        \n",
    "        # find the index of all x data within segment k\n",
    "        indx = (x >= xk )*( x <= xkp)\n",
    "        \n",
    "        # calculate the fractional distance within the element (l_1(x))\n",
    "        s = (x[indx] - xk)/(xkp - xk)\n",
    "        \n",
    "        # linear interpolation within the segment \n",
    "        y[indx] = y_data[k]*(1. - s) + y_data[k+1]*s\n",
    "        \n",
    "    return y\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "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",
    "\n",
    "# Interpolating polynmial through the data\n",
    "P_lagrange = poly_interpolant(x, data)\n",
    "\n",
    "# Piecewise linear interpolatnt between the data\n",
    "P_linear = piece_wise_linear_interp(x, data)\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "scrolled": false,
    "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 = len(data)-1\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')\n",
    "axes.plot(x, P_lagrange, 'b--',label='$P_{{{name}}}$'.format(name=N))\n",
    "axes.plot(x, P_linear, 'r',label='Piecewise Linear')\n",
    "\n",
    "axes.set_title(\"Interpolated Data - $C^0$ Linear\",fontsize=18)\n",
    "axes.set_xlabel(\"x\",fontsize=16)\n",
    "axes.set_ylabel(\"$P_1(x)$\",fontsize=16)\n",
    "axes.set_xlim([0.0, 10.0])\n",
    "axes.set_ylim([-4.0, 15.0])\n",
    "axes.legend(loc='best',fontsize=14)\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Piece-Wise Quadratic $C^0$ Polynomials\n",
    "\n",
    "Use every three points $(x_{k+1}, y_{k+1})$, $(x_{k}, y_{k})$, and $(x_{k-1}, y_{k-1})$, to 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": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "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(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(data[:,0], data[:,1], 'ko')\n",
    "axes.plot(x, P_lagrange, 'b--',label='$P_{{{name}}}$'.format(name=N))\n",
    "axes.plot(x, P_quadratic, 'r',label='Piecewise quadratic')\n",
    "\n",
    "axes.set_title(\"Interpolated Data - $C^0$ Quadratic\",fontsize=18)\n",
    "axes.set_xlabel(\"x\",fontsize=16)\n",
    "axes.set_ylabel(\"$P_2(x)$\",fontsize=16)\n",
    "axes.set_xlim([0.0, 10.0])\n",
    "axes.set_ylim([-4.0, 15.0])\n",
    "axes.legend(loc='best',fontsize=14)\n",
    "axes.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Transformation of intervals\n",
    "\n",
    "The previous algorithms are quite direct but can look a bit messy because every interval (or set of points) appears different...An important trick that is used in finite-element analysis is to transform each interval (or element) \n",
    "\n",
    "$$\n",
    "    x \\in[x_k, x_{k+1}]\n",
    "$$ \n",
    "\n",
    "to the unit interval\n",
    "\n",
    "$$\n",
    "    s \\in [0,1]\n",
    "$$\n",
    "\n",
    "by the affine transformation\n",
    "\n",
    "$$ x(s) = x_k + s(x_{k+1} - x_k)$$\n",
    "\n",
    "and do all the interpolation in the transformed frame."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Example: Linear Interpolation\n",
    "\n",
    "for Linear interpolation over an arbitrary interval $[x_k, x_{k+1}]$ we can use a Lagrange interpolant\n",
    "\n",
    "$$\n",
    "    P_k(x) = y_k\\ell_0(x) + y_{k+1}\\ell_1(x)\n",
    "$$\n",
    "\n",
    "where \n",
    "$$ \n",
    "    \\ell_0(x) = \\frac{x - x_{k+1}}{x_k - x_{k+1}},~\\quad~  \\ell_1(x) = \\frac{x - x_{k}}{x_{k+1} - x_{k}}\n",
    "$$\n",
    "\n",
    "and do all the interpolation in the transformed frame."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Substituting \n",
    "$$ x(s) = x_k + s(x_{k+1} - x_k)$$\n",
    "\n",
    "into the definitions of $\\ell_0(x)$ and $\\ell_1(x)$ show that within this element\n",
    "\n",
    "$$\n",
    "    \\ell_0(s) = 1-s,~\\quad~\\ell_1(s) = s\n",
    "$$\n",
    "\n",
    "and \n",
    "$$\n",
    "    P_k(s) = y_k(1-s)  + y_{k+1}s\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "xk = [ .5, 1., 3., 5.]\n",
    "yk = [ .5, 2., 2.5, 1. ]\n",
    "ell0 = lambda x: (x - x1)/(x0 - x1)\n",
    "ell1 = lambda x: (x - x0)/(x1 - x0)\n",
    "P1 = lambda x:  y0*ell0(x) + y1*ell1(x)\n",
    "\n",
    "x = numpy.linspace(x0,x1)\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1,1,1)\n",
    "axes.plot(x,ell0(x),'r',label='$\\ell_0$')\n",
    "axes.plot(x,ell1(x),'b',label='$\\ell_1$')\n",
    "axes.plot(x,P1(x),'g',label='$P_1(x)$')\n",
    "axes.plot((x0,x1),(y0,y1),'go')\n",
    "axes.set_xlabel('x')\n",
    "axes.set_xlim(0.,5.)\n",
    "axes.grid()\n",
    "axes.legend(loc='best',fontsize=16)\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",
    "where $d^*_k=P'(0)$ is the slope at point $k$ in the transformed frame (likewise $d^*_{k+1}=P'(1)$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### 2 Questions\n",
    "\n",
    "* what is the relationship between $d^*_k = P'(s)$ in the transformed frame $d_k=P'(x)$ in the original frame?\n",
    "* how to choose $d_k$ in general?\n"
   ]
  },
  {
   "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": [
    "from scipy.interpolate import pchip_interpolate\n",
    " \n",
    "data = numpy.array([[1.0, 3.0], [2.0, 1.0], [3.5, 4.0], [5.0, 4.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 = 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",
    "# Interpolating polynmial through the data\n",
    "P_lagrange = poly_interpolant(x, data)\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "axes.plot(data[:,0], data[:,1], 'ko')\n",
    "axes.plot(x, P_lagrange, 'b--',label='$P_{{{name}}}$'.format(name=N))\n",
    "axes.plot(x, P_pchip, 'r',label='PChip') \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.legend()\n",
    "axes.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "#### Cubic Splines\n",
    "\n",
    "Enforces 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",
    "or rearrange knowns and unknowns to get.\n",
    "\n",
    "$$\n",
    "    2 d_{k-1} + 8 d_k + 2 d_{k+1} =  6 (y_{k+1} - y_{k-1}) \n",
    "$$ \n",
    "\n",
    "We now have constraints on choosing the $d_k$ values for all interior values of $k$.  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": [
    "from scipy.interpolate import UnivariateSpline\n",
    "\n",
    "data = numpy.array([[1.0, 3.0], [2.0, 1.0], [3.5, 4.0], [5.0, 4.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 = UnivariateSpline(data[:, 0], data[:, 1], s=0) \n",
    "\n",
    "# C^1 Piece-wise PCHIP for reference\n",
    "P_pchip = pchip_interpolate(data[:, 0], data[:, 1], x)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Plot and compare to pchip\n",
    "\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.legend(loc='best')\n",
    "axes.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide_input": false,
    "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": "fragment"
    }
   },
   "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",
    "print('Normal Equations: p = {}'.format(p))\n",
    "p = numpy.linalg.lstsq(A, y, rcond=None)[0]\n",
    "print('Numpy Lstsq     : p = {}'.format(p))\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": true,
    "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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "#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] * numpy.tanh(x) \n",
    "E = numpy.linalg.norm(y - f(x), ord=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": true,
    "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/).)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Higher Dimensional interpolation\n",
    "\n",
    "Much of the machinery for interpolation in 1-dimension can be extended in straightforward ways to higher dimensions.  Again, many of these methods are provided by `scipy.interpolate` but their mechanisms are straightforward to derive.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Here we will consider two relative simple extensions\n",
    "\n",
    "* bi-polynomial (bi-linear) interpolation on a regular mesh\n",
    "* linear interpolation on Triangles"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.7"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}