chapter4-3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c52c1c95",
   "metadata": {},
   "source": [
    "# 4.3 Hermite interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dcd3b3e",
   "metadata": {},
   "source": [
    "In polynomial interpolation, our starting point has been the $x$- and $y$-coordinates of some data we want to interpolate. Suppose, in addition, we know the derivative of the underlying function at these $x$-coordinates. Our new data set has the following form.\n",
    "\n",
    "**Data:**\n",
    "- $x_{0},x_{1},\\ldots,x_{n}$\n",
    "- $y_{0},y_{1},\\ldots,y_{n};\\, y_{i}=f(x_{i})$\n",
    "- $y_{0}',y_{1}',\\ldots,y_{n}';\\, y'_{i}=f'(x_{i})$\n",
    "\n",
    "We seek a polynomial that fits the $y$ and $y'$ values, that is, we seek a polynomial $H(x)$ such that $H(x_{i})=y_{i}$ and $H'(x_{i})=y'_{i}$,\n",
    "$i=0,1,\\ldots,n.$ This makes $2n+2$ equations, and if we let\n",
    "\\begin{equation*}\n",
    "H(x)=a_{0}+a_{1}x+\\cdots+a_{2n+1}x^{2n+1},\n",
    "\\end{equation*}\n",
    "then there are $2n+2$ unknowns, $a_0,\\ldots,a_{2n+1}$, to solve for. The following theorem shows that there is a unique solution to this system of equations; a proof can be found in Burden et al. [2016]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c64873d",
   "metadata": {},
   "source": [
    "\\begin{theorem}\\label{theo:chap3_thm61}\n",
    "If $f \\in C^1[a,b]$ and $x_{0},\\ldots,x_{n}\\in[a,b]$\n",
    "are distinct, then there is a unique polynomial $H_{2n+1}(x)$, of degree at most $2n+1$, agreeing with\n",
    "$f$ and $f'$ at $x_{0},\\ldots,x_{n}$. The polynomial can be written as:\n",
    "\\begin{equation*}\n",
    "H_{2n+1}(x)=\\sum_{i=0}^{n}y_{i}h_{i}(x)+\\sum_{i=0}^{n}y_{i}'\\tilde{h_{i}}(x)\n",
    "\\end{equation*}\n",
    "where\n",
    "\\begin{align*}\n",
    "h_{i}(x) & =\\left(1-2(x-x_{i})l'_{i}(x_{i})\\right)(l_{i}(x))^{2}\\\\\n",
    "\\widetilde{h}_{i}(x) & =(x-x_{i})(l_{i}(x))^{2}.\n",
    "\\end{align*}\n",
    "Here $l_{i}(x)$ is the $i$th Lagrange basis function for the nodes\n",
    "$x_{0},\\ldots,x_{n}$, and $l'_{i}(x)$ is its derivative. $H_{2n+1}(x)$ is called the Hermite interpolating\n",
    "polynomial.\n",
    "\\end{theorem}\n",
    "\n",
    "The only difference between Hermite interpolation and polynomial interpolation is that in the former, we have the derivative information, which can go a long way in capturing the shape of the underlying function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5b06f27",
   "metadata": {},
   "source": [
    "\\begin{example}\\label{example:chap3_exa62}\n",
    "We want to interpolate the following data:\n",
    "\\begin{align*}\n",
    "x \\text{-coordinates}&: -1.5, 1.6, 4.7\\\\\n",
    "y \\text{-coordinates}&: 0.071,-0.029,-0.012.\n",
    "\\end{align*}\n",
    "The underlying function the data comes from is $\\cos x$, but we pretend we do not know this. The figure below plots the underlying function, the data, and the polynomial interpolant for the data. Clearly, the polynomial interpolant does not come close to giving a good approximation to the underlying function $\\cos x$.\n",
    "\\end{example}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4feb8e29",
   "metadata": {},
   "source": [
    "![](./images/hermite_cos1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d10be8e0",
   "metadata": {},
   "source": [
    "Now let's assume we know the derivative of the underlying function at these nodes:\n",
    "\\begin{align*}\n",
    "x \\text{-coordinates}&: -1.5, 1.6, 4.7\\\\\n",
    "y \\text{-coordinates}&: 0.071,-0.029,-0.012\\\\\n",
    "y' \\text{-values}&: 1, -1, 1.\n",
    "\\end{align*}\n",
    "\n",
    "We then construct the Hermite interpolating polynomial, incorporating the derivative information. The figure below plots the Hermite interpolating polynomial, together with the polynomial interpolant, and the underlying function.\n",
    "\n",
    "It is visually difficult to separate the Hermite interpolating polynomial from the underlying function $\\cos x$ in the following figure. Going from polynomial interpolation to Hermite interpolation results in rather dramatic improvement in approximating the underlying function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb50d320",
   "metadata": {},
   "source": [
    "![](./images/hermite_cos2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bacc99b5",
   "metadata": {},
   "source": [
    "## Computing the Hermite polynomial"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fdee514",
   "metadata": {},
   "source": [
    "We do not use Theorem \\ref{theo:chap3_thm61} to compute the Hermite polynomial: there\n",
    "is a more efficient method using divided differences for this computation.\n",
    "\n",
    "We start with the data:\n",
    "\\begin{align*}\n",
    "&x_{0},x_{1},\\ldots,x_{n}\\\\\n",
    "&y_{0},y_{1},\\ldots,y_{n}; y_{i}=f(x_{i})\\\\\n",
    "&y_{0}',y_{1}',\\ldots,y_{n}';y'_{i}=f'(x_{i})\\\\\n",
    "\\end{align*}\n",
    "and define a sequence $z_{0},z_{1},\\ldots,z_{2n+1}$ by\n",
    "\\begin{align*}\n",
    "z_{0} & =x_{0},z_{2}=x_{1},z_{4}=x_{2},\\ldots,z_{2n}=x_{n}\\\\\n",
    "z_{1} & =x_{0},z_{3}=x_{1},z_{5}=x_{2},\\ldots,z_{2n+1}=x_{n}\n",
    "\\end{align*}\n",
    "i.e., $z_{2i}=z_{2i+1}=x_{i},$ for $i=0,1,\\ldots,n.$\n",
    "\n",
    "Then the Hermite polynomial can be written as:\n",
    "\\begin{align*}\n",
    "H_{2n+1}(x) =& f[z_{0}]+f[z_{0},z_{1}](x-z_{0})+f[z_{0},z_{1},z_{2}](x-z_{0})(x-z_{1})+\\\\\n",
    " & \\cdots+f[z_{0},z_{1},\\ldots,z_{2n+1}](x-z_{0})\\cdots(x-z_{2n})\\\\\n",
    "  =&f[z_{0}]+\\sum_{i=1}^{2n+1}f[z_{0},\\ldots,z_{i}](x-z_{0})(x-z_{1})\\cdots(x-z_{i-1}).\n",
    "\\end{align*}\n",
    "There is a little problem with some of the first divided differences above: they are undefined! Observe that\n",
    "\\begin{equation*}\n",
    "f[z_{0},z_{1}]=f[x_{0},x_{0}]=\\frac{f(x_{0})-f(x_{0})}{x_{0}-x_{0}}\n",
    "\\end{equation*}\n",
    "or, in general,\n",
    "\\begin{equation*}\n",
    "f[z_{2i},z_{2i+1}]=f[x_{i},x_{i}]=\\frac{f(x_{i})-f(x_{i})}{x_{i}-x_{i}}\n",
    "\\end{equation*}\n",
    "for $i=0,\\ldots,n.$\n",
    "\n",
    "From Theorem 1 Chapter4-2, we know $f[x_{0},\\ldots,x_{n}]=\\frac{f^{(n)}(\\xi)}{n!}$\n",
    "for some $\\xi$ between the min and max of $x_{0},\\ldots,x_{n}.$ From\n",
    "a classical result by Hermite \\& Gennochi (see Atkinson [1989], page 144), divided differences are\n",
    "continuous functions of their variables $x_{0},\\ldots,x_{n}.$ This implies\n",
    "we can take the limit of the above result as $x_{i}\\rightarrow x_{0}$\n",
    "for all $i,$ which results in\n",
    "\\begin{equation*}\n",
    "f[x_{0},\\ldots,x_{0}]=\\frac{f^{(n)}(x_{0})}{n!}.\n",
    "\\end{equation*}\n",
    "Therefore in the Hermite polynomial coefficient calculations, we will\n",
    "put\n",
    "\\begin{equation*}\n",
    "f[z_{2i},z_{2i+1}]=f[x_{i},x_{i}]=f'(x_{i})=y'_{i}\n",
    "\\end{equation*}\n",
    "for $i=0,1,\\ldots,n.$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6183c132",
   "metadata": {},
   "source": [
    "\\begin{example}\\label{example:chap3_exa63}\n",
    "Let's compute the Hermite polynomial of Example \\ref{example:chap3_exa62}. The data is:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{array}{cccc}\n",
    "i & x_{i} & y_{i} & y_{i}'\\\\\n",
    "0 & -1.5 & 0.071 & 1\\\\\n",
    "1 & 1.6 & -0.029 & -1\\\\\n",
    "2 & 4.7 & -0.012 & 1\\\\\n",
    "\\end{array}\n",
    "\\end{equation*}\n",
    "Here $n=2,$ and $2n+1=5,$ so the Hermite polynomial is\n",
    "\\begin{equation*}\n",
    "H_{5}(x)=f[z_{0}]+\\sum_{i=1}^{5}f[z_{0},\\ldots,z_{i}](x-z_{0})\\cdots(x-z_{i-1}).\n",
    "\\end{equation*}\n",
    "The divided differences are:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{array}{ccccccc}\n",
    "{\\small{}z} & {\\small{}f(z)} & {\\small{}\\text{1st diff}} & {\\small{}\\text{2nd diff}} & {\\small{}\\text{3rd diff}} & {\\small{}\\text{4th diff}} & {\\small{}\\text{5th diff}}\\\\\n",
    "\\hline\n",
    "{\\small{}z_{0}=-1.5} & {\\small{}\\color{blue}{0.071}} &  &  &  &  & \\\\\n",
    " &  & {\\small{}f'(z_{0})=\\color{blue}{1}} &  &  &  & \\\\\n",
    "{\\small{}z_{1}=-1.5} & {\\small{}0.071} &  & {\\small{}\\frac{-0.032-1}{1.6+1.5}=\\color{blue}{-0.33}} &  &  & \\\\\n",
    " &  & {\\small{}f[z_{1},z_{2}]=-0.032} &  & {\\small{}\\color{blue}{0.0065}} &  & \\\\\n",
    "{\\small{}z_{2}=1.6} & {\\small{}-0.029} &  & {\\small{}\\frac{-1+0.032}{1.6+1.5}=-0.31} &  & {\\small{}\\color{blue}{0.015}} & \\\\\n",
    " &  & {\\small{}f'(z_{2})=-1} &  & {\\small{}0.10} &  & {\\small{}\\color{blue}{-0.005}}\\\\\n",
    "{\\small{}z_{3}=1.6} & {\\small{}-0.029} &  & {\\small{}\\frac{0.0055+1}{4.7-1.6}=0.32} &  & {\\small{}-0.016} & \\\\\n",
    " &  & {\\small{}f[z_{3},z_{4}]=0.0055} &  & {\\small{}0} &  & \\\\\n",
    "{\\small{}z_{4}=4.7} & {\\small{}-0.012} &  & {\\small{}\\frac{1-0.0055}{4.7-1.6}=0.32} &  &  & \\\\\n",
    " &  & {\\small{}f'(z_{4})=1} &  &  &  & \\\\\n",
    "{\\small{}z_{5}=4.7} & {\\small{}-0.012} &  &  &  &  & \\\\\n",
    "\\hline\n",
    "\\end{array}\n",
    "\\end{equation*}\n",
    "\n",
    "Therefore, the Hermite polynomial is:\n",
    "\\begin{align*}\n",
    "H_{5}(x) =&\\color{blue}{0.071}+\\color{blue}{1}(x+1.5)\\color{blue}{-0.33}(x+1.5)^{2}+\\color{blue}{0.0065}(x+1.5)^{2}(x-1.6)+\\\\\n",
    " & +\\color{blue}{0.015}(x+1.5)^{2}(x-1.6)^{2}\\color{blue}{-0.005}(x+1.5)^{2}(x-1.6)^{2}(x-4.7).\n",
    "\\end{align*}\n",
    "\\end{example}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2477524b",
   "metadata": {},
   "source": [
    "## Python code for computing Hermite interpolating polynomial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39f959c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1eec6415",
   "metadata": {},
   "source": [
    " The following function **hdiff** computes the divided\n",
    "differences needed for Hermite interpolation. It is based on the\n",
    "function **diff** for computing divided differences for Newton\n",
    "interpolation. The inputs to **hdiff** are the $x$-coordinates,\n",
    "the $y$-coordinates, and the derivatives **yprime}**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd922ab3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def hdiff(x, y, yprime):\n",
    "    m = x.size # here m is the number of data points. Note n=m-1\n",
    "    # and 2n+1=2m-1\n",
    "    l = 2*m\n",
    "    z = np.zeros(l)\n",
    "    a = np.zeros(l)\n",
    "    for i in range(m):\n",
    "        z[2*i] = x[i]\n",
    "        z[2*i+1] = x[i]\n",
    "    for i in range(m):\n",
    "        a[2*i] = y[i]\n",
    "        a[2*i+1] = y[i]\n",
    "    for i in np.flip(np.arange(1, m)): # computes the first divided\n",
    "        # differences using derivatives\n",
    "        a[2*i+1] = yprime[i]\n",
    "        a[2*i] = (a[2*i]-a[2*i-1]) / (z[2*i]-z[2*i-1])\n",
    "    a[1] = yprime[0]\n",
    "    for j in range(2, l): # computes the rest of the divided differences\n",
    "        for i in np.flip(np.arange(j, l)):\n",
    "            a[i]=(a[i]-a[i-1]) / (z[i]-z[i-j])\n",
    "    return a"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d65eff21",
   "metadata": {},
   "source": [
    "Let's compute the divided differences of \\ref{example:chap3_exa62}."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "737db491",
   "metadata": {},
   "outputs": [],
   "source": [
    "hdiff(np.array([-1.5, 1.6, 4.7]),\n",
    "      np.array([0.071,-0.029,-0.012]), \n",
    "      np.array([1,-1,1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e74585f7",
   "metadata": {},
   "source": [
    "Note that in the hand-calculations of \\ref{example:chap3_exa63}, where two-digit rounding\n",
    "was used, we obtained 0.0065 as the first third divided difference. In the\n",
    "Python output above, this divided difference is 0.0067.\n",
    "\n",
    "The following function computes the Hermite interpolating polynomial,\n",
    "using the divided differences obtained from **hdiff**, and then\n",
    "evaluates the polynomial at $w$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76180e57",
   "metadata": {},
   "outputs": [],
   "source": [
    "def hermite(x, y, yprime, w):\n",
    "    m = x.size # here m is the number of data points. not the\n",
    "    # degree of the polynomial\n",
    "    a = hdiff(x, y, yprime)\n",
    "    z = np.zeros(2*m)\n",
    "    for i in range(m):\n",
    "        z[2*i] = x[i]\n",
    "        z[2*i+1] = x[i]\n",
    "    sum = a[0]\n",
    "    pr = 1.0\n",
    "    for j in range(2*m-1):\n",
    "        pr *= w-z[j]\n",
    "        sum += a[j+1]*pr\n",
    "    return sum"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a72d697",
   "metadata": {},
   "source": [
    " Let's recreate the Hermite interpolating polynomial plot of \\ref{example:chap3_exa62}."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca3829bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "xaxis = np.linspace(-np.pi/2, 3*np.pi/2, 120)\n",
    "x = np.array([-1.5, 1.6, 4.7])\n",
    "y = np.array([0.071,-0.029,-0.012])\n",
    "yprime = np.array([1, -1, 1])\n",
    "funct = np.cos(xaxis)\n",
    "interp = hermite(x, y, yprime, xaxis)\n",
    "plt.plot(xaxis, interp, label='Hermite interpolation')\n",
    "plt.plot(xaxis, funct, label=\"cos(x)\")\n",
    "plt.plot(x, y, 'o', label='data')\n",
    "plt.legend(loc='upper right');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87402cd3",
   "metadata": {},
   "source": [
    "**Exercise 4.3-1**\n",
    "\n",
    "The following table gives the values of $y=f(x)$ and $y'=f'(x)$ where $f(x)=e^x+\\sin 10x$. Compute the Hermite interpolating polynomial and the polynomial interpolant for the data in the table. Plot the two interpolating polynomials together with $f(x)=e^x+\\sin 10x$ on $(0,3)$.\n",
    "\\begin{equation*}\n",
    "\\begin{array}{ccccccc}\n",
    "x & 0 & 0.4 & 1 & 2 & 2.6 & 3\\\\\n",
    "\\hline\n",
    "y & 1 & 0.735 & 2.17 & 8.30 & 14.2 & 19.1\\\\\n",
    "\\hline\n",
    "y' & 11 & -5.04 & -5.67 & 11.5 & 19.9 &21.6\\\\\n",
    "\\end{array}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd2bdeea",
   "metadata": {},
   "source": [
    "**References**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89bd1d53",
   "metadata": {},
   "source": [
    "K.E. Atkinson. An Introduction to Numerical Analysis. John Wiley & Sons, 2nd edition, 1989.\n",
    "\n",
    "R.L. Burden, D. Faires, and A.M. Burden. Numerical Analysis. Cengate, 10th edition, 2016."
   ]
  }
 ],
 "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.9.7"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}