poly_interp2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "98a1809d",
   "metadata": {},
   "source": [
    "# 3. Interpolation II"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0440e900",
   "metadata": {},
   "source": [
    "Earlier, we discussed two options for the choice of interpolation basis (monomial/Vandermonde and Lagrange). Now let's take a look at the other parameter of interpolation: the choice of collocation nodes. \n",
    "\n",
    "We established that the error in numerical interpolation comes from two terms,\n",
    "\n",
    "$$ ||f - \\tilde{p}_n|| \\leq ||f - p_n|| + ||p_n - \\tilde{p}_n||,$$ \n",
    "\n",
    "where the former is what we'll focus on now.\n",
    "To minimize errors from the second term (numerical evaluation of the interpolating polynomial), we'll use the Lagrange basis and define\n",
    "\n",
    "$$ p(x) = p_n(x) := (L_n f)(x),$$\n",
    "\n",
    "where $L_n$ is the Lagrange interpolating operator. Let's start with a _bad_ example: equispaced nodes on the interval $[-1, 1]$ (which we can always rescale), $x_j = -1 + \\tfrac{2j}{n}$ for $j = 0, 1, \\ldots, n$. In the example below, I attempt to interpolate $f(x) = \\frac{1}{1 + 25x^2}$ using an increasing number of interpolation nodes. \n",
    "\n",
    "```{warning}\n",
    "Click the rocket (launch) button on the top right of this page to interact with the plots below.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d80bf8f3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1adb613652b749f9b6477132a2dddd06",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=42, description='n', max=50, min=1), FloatSlider(value=0.2, description=…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.draw_equiint(n, x0)>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.interpolate as sciint\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact, IntSlider, FloatSlider\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore') # Don't try this at home :)\n",
    "\n",
    "\n",
    "# The function we wish to interpolate on [-1, 1]\n",
    "f = lambda x, x0: 1/(1 + (x/x0)**2)\n",
    "\n",
    "# We choose the collocation nodes to be equidistant \n",
    "# nodes = lambda n: np.linspace(-1, 1, n)\n",
    "# Try commenting the above line and uncommenting the below one \n",
    "# for Chebyshev nodes\n",
    "nodes = lambda n: np.cos(np.linspace(0, np.pi, n))\n",
    "\n",
    "# Target nodes\n",
    "x_star = np.linspace(-1, 1, 1000)\n",
    "\n",
    "# Barycentric interpolation from given nodes\n",
    "baryint = lambda nodes, n, x0, f: sciint.barycentric_interpolate(nodes(n), f(nodes(n), x0), x_star)\n",
    "\n",
    "def draw_equiint(n, x0):\n",
    "    fig, ax = plt.subplots(1, 1, figsize = (6, 6))\n",
    "    ax.plot(x_star, f(x_star, x0), '--', label = '$f(x)$', color = 'black')\n",
    "    ax.set_ylim(1/(1+(1/x0)**2)-0.1, 1.1)\n",
    "    ax.set_xlim(-1.1, 1.1)\n",
    "    ax.set_xlabel(\"$x$\")\n",
    "    ax.plot(x_star, baryint(nodes, n, x0, f), label = \"$f_{\"+\"{}\".format(n)+\"}(x)$\")\n",
    "    ax.plot(nodes(n), f(nodes(n), x0), '.')\n",
    "    ax.legend(loc = 'lower center')\n",
    "\n",
    "\n",
    "interact(draw_equiint, n = IntSlider(min = 1, max = 50, value = 42),\\\n",
    "                     x0 = FloatSlider(min = 0.2 , max = 3, step = 0.2),\\\n",
    "                     continuous_update = False)            "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b0dc75c",
   "metadata": {},
   "source": [
    "```{admonition} Question\n",
    "What do you observe as $n$ increases with equispaced nodes? What happens if you switch to Chebyshev nodes?\n",
    "```\n",
    "\n",
    "With increasing $n$, we see that\n",
    "\n",
    "$$ |f - L_nf| \\begin{cases}\n",
    "&\\to 0 \\text{ as } n \\to \\infty \\text{ in the central region,}\\\\\n",
    "&\\text{ oscillates between interpolation nodes near $\\pm 1$ as } n \\to \\infty. \n",
    "\\end{cases} $$\n",
    "\n",
    "This is called the **Runge phenomenon**. \n",
    "If we \"bunch\" interpolation nodes near the edges of the interval instead, e.g. use Chebyshev nodes\n",
    "$$ x_j = - \\cos(j\\pi/n) \\quad j = 0, 1, \\ldots, n,  $$\n",
    "then we observe\n",
    "\n",
    "$$ \\max_{x \\in [-1, 1]}| f - L_nf | \\equiv || f - L_nf||_{L_\\infty} \\to 0 \\quad \\text{as } n \\to \\infty,$$\n",
    "\n",
    "i.e. we see **uniform convergence**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b41c8b2",
   "metadata": {},
   "source": [
    "If you now started to change the $x_0$ parameter in the above $f(x) = \\frac{1}{1 + \\left(\\frac{x}{x_0}\\right)^2}$ from its default value $x_0 = 0.2$ to something like $x_0 = 2.0$, then even at large $n$, interpolation with equispaced nodes appears to yield an accurate answer and no Runge phenomenon is seen.\n",
    "\n",
    "From these experiments it's clear that:\n",
    "- the choice of collocation nodes has an effect on the convergence rate of interpolation,\n",
    "- and so do properties of the function $f(x)$, and it turns out that for some choices of $f$, we don't have to worry about non-convergence. \n",
    "\n",
    "```{admonition} Question\n",
    "Can you guess what property of $f(x)$ determines the convergence rate of interpolation from the above? (Hint: what did $x_0$ control?)\n",
    "```\n",
    "To see this, let's derive a (loose) error bound."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a792178",
   "metadata": {},
   "source": [
    "## 3.1 Error of interpolation for analytic $f$\n",
    "\n",
    "```{admonition} Theorem\n",
    ":class: hint\n",
    "If we assume nothing about the collocation nodes, but that $f \\in C^{(n+1)}[a, b]$ ($n + 1$ times continuously differentiable on the interpolation interval), then \n",
    "$$ || f - L_nf ||_{L_{\\infty}} \\leq \\frac{||f^{(n+1)}||_{L_{\\infty}}}{(n+1)!} H^{n+1}, $$\n",
    "where $H$ is the length of the interpolation interval $[a, b]$, and $f^{(n+1)}$ is the $(n+1)$th derivative of $f$.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee3b3f9c",
   "metadata": {},
   "source": [
    "**Proof**: First consider the **remainder** function, defined as\n",
    "\n",
    "$$ R_nf := f - L_nf. $$\n",
    "\n",
    "We'll first prove that \n",
    "\n",
    "$$ (R_n f)(x) = \\frac{f^{(n+1)}(\\xi) }{(n+1)!} \\prod_{j = 0}^{n} (x - x_j), $$\n",
    "\n",
    "for $x \\in [a, b]$ and some $\\xi \\in [a, b]$.\n",
    "The expression is clearly true if $x = x_j$ (since both sides are zero), so let's assume $x$ is not one of the nodes. Then let\n",
    "\n",
    "$$ q_{n+1}(x) := \\prod_{j = 0}^n (x - x_j),$$\n",
    "\n",
    "and consider the function \n",
    "\n",
    "$$ g(y) := f(y) - (L_nf )(y) - q_{n+1}(y) \\frac{f(x) - (L_n f)(x)}{q_{n+1}(x)}. $$\n",
    "\n",
    "If $f$ was $C^{n+1}$, then so if $g$, since $L_nf$ and $q_{n+1}$ are polynomials (and hence $C^{\\infty}$). $g$ therefore has at least $n+2$ roots on the inteprolation interval, $g'$ has at least $n+1$, and by induction, $g^{(n+1)}$ has at least one. Let that root be at $\\xi$. Differentiating both sides of the expression with respect to $y$ $(n+1)$ times and evaluating at $\\xi$ gives\n",
    "\n",
    "$$ 0 = f^{(n+1)}(\\xi) - 0 + (n+1)!\\frac{(R_nf)(x)}{q_{n+1}(x)}, $$\n",
    "\n",
    "where I used the definition of the remainder and that $L_nf$ is a degree-$n$ polynomial. Rearranging gives the desider expression for $R_n f$. Now using \n",
    "\n",
    "$$ \\prod_{j = 0}^n (x - x_j) \\leq (b - a)^{n+1} = H^{n+1} $$\n",
    "\n",
    "and taking $L_{\\infty}$ norms, we complete the proof.\n",
    "\n",
    "Given an interpolation interval of length $H$, it's now clear that the factor that determines the convergence rate is $||f^{(n+1)}/(n+1)!||$...\n",
    "\n",
    "```{admonition} Question\n",
    "What is this?\n",
    "```\n",
    "... the magnitude of the $(n+1)$th Taylor series coefficient of $f$. Just how large can it get?\n",
    "\n",
    "We are now going to make an even stricter assumption of $f$: instead of asking for it to be $C^{n+1}$, we'll assume it's **analytic** in a region of the complex plane.\n",
    "\n",
    "```{admonition} Definition\n",
    ":class: hint\n",
    "$f$ is analytic in an open domain if at every point in the domain, it is locally given by a convergent power series. \n",
    "```\n",
    "\n",
    "For example, if $f$ is analytic at $0$, then there exists $\\rho > 0$ such that \n",
    "\n",
    "$$ f(z) = \\sum_{n = 0}^\\infty a_n z^n, \\quad a_n = \\frac{f^{(n)}(0)}{n!}, $$\n",
    "\n",
    "is an an absolutely convergent series for all $z$ in the disc $|z| < \\rho$, and is absolutely divergent outside of it ($|z| > \\rho$), with $f$ being analytic in $|z| < \\rho$.\n",
    "For absolute convergence $\\forall \\varepsilon > 0$, in $|z| = \\rho - \\varepsilon$ the $a_n$ must decay as $n \\to \\infty$, so are bounded by some constant $c$:\n",
    "\n",
    "$$|a_n|(\\rho - \\varepsilon)^n \\leq c, $$\n",
    "\n",
    "i.e. $|a_n| \\leq \\frac{c}{(\\rho - \\varepsilon)^n}$.\n",
    "\n",
    "Let's assume $f$ is analytic in a \"stadium\" defined as $\\{ z \\in \\mathbb{C}, \\text{dist}(x, [a, b]) \\leq R\\}$, a set of points that are within a distance $R$ of the interpolation interval. Then the Taylor coefficients satisfy \n",
    "\n",
    "$$|a_n| \\leq \\frac{c}{R^n} $$\n",
    "\n",
    "(with the expansion center being in $[a, b]$). Then \n",
    "\n",
    "$$ ||f - L_n f||_{L_{\\infty}}  \\leq \\frac{c}{R^n}H^{n+1} = \\tilde{c} \\left( \\frac{H}{R} \\right)^n.$$\n",
    "\n",
    "So for functions $f$ that are analytic in some region, it's the width of this region relative to the interpolation interval that determines the convergence rate. If it is larger than the interval width, we'll get _geometric_/_exponential_ convergence _regardless_ of the choice of nodes!\n",
    "\n",
    "```{admonition} Question\n",
    "What determines the region of analyticity of $f$ in the code experiment?\n",
    "```\n",
    "\n",
    "That region in turn is determined by the nearest singularity (pole) in $f$ (unless there is something more pathological going on), which were at $\\pm i x_0$, a distance $x_0$ away from $[a, b]$.  With this in mind, revisit the experiment and play around with the $x_0$ setting.\n",
    "\n",
    "There are, however, some good and bad news (standard results) about the convergence of interpolation and choice of nodes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "972d555d",
   "metadata": {},
   "source": [
    "## 3.2 The good and the bad\n",
    "\n",
    "I'll state these without proof.\n",
    "\n",
    "```{admonition} Theorem (Faber)\n",
    ":class: hint\n",
    "If we constructed a sequence of interpolating operators $L_n$, each with nodes $\\{ x_j^{(n)}\\}_{j = 0}^n$ for $n = 0, 1, \\ldots,$ then there always $\\exists f \\in C[a, b]$ such that $L_n f$ does not converge to $f$ uniformly on $[a, b]$.\n",
    "```\n",
    "On the contrary,\n",
    "\n",
    "```{admonition} Theorem (Marcinkiewicz)\n",
    ":class: hint\n",
    "For each function $f \\in C[a, b]$, $\\exists$ a sequence of interpolating nodes $\\{ x_j^{(n)}\\}_{j = 0}^n$ for $n = 0, 1, \\ldots,$, such that for the resulting sequence of inteprolating operators $L_n$, $L_n f$ converges to $f$ uniformly on $[a, b]$.\n",
    "```\n",
    "\n",
    "With this in mind, why is it worth clustering the interpolating nodes at the edges of the interval?\n",
    "\n",
    "## 3.2 Why cluster nodes at the edges?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "493287e6",
   "metadata": {},
   "source": [
    "To explain this, we invoke some physics (potential theory) that will be useful later on when considering boundary intergal equations. \n",
    "\n",
    "Consider the **potential** function\n",
    "\n",
    "$$\\phi_{(n+1)}(z) := \\frac{1}{n+1} \\ln \\left| \\prod_{j = 0}^n (z - x_j) \\right|  = - \\frac{1}{n+1} \\sum_{j = 0}^n \\ln \\frac{1}{|z - x_j|},$$\n",
    "\n",
    "where the term in the absolute value sign on the left is none other than $q_{n+1}(z)$.\n",
    "\n",
    "```{admonition} Question\n",
    "What is this interpreted as in physics? (Hint: think electrostatics in 2D.)\n",
    "```\n",
    "\n",
    "Physically, the right-hand-side denotes the potential generated by $n+1$ charges located at each of the $x_j$, in two dimensions, $\\mathbb{R}^2 = \\mathbb{C}$. Each of the charges has strength $1/(n+1)$, so that they add up to $1$. Rearranging gives\n",
    "\n",
    "$$q_{n+1}(z) = e^{(n+1)\\phi_{(n+1)}(z)}. $$\n",
    "\n",
    "As $n \\to \\infty$, let's say these discrete charges tend to a fixed density function $\\rho(x) > 0$ on the interval $[-1, 1]$. We keep them normalized so that the total charge is still $1$,\n",
    "\n",
    "$$ \\int_{-1}^1 \\rho(x)\\mathrm{d}x = 1. $$\n",
    "\n",
    "What happens to the potential they generate? It becomes \n",
    "\n",
    "$$ \\phi_{(n+1)}(z) \\to \\phi(z) = - \\int_{-1}^1 \\rho(x) \\ln \\frac{1}{|z - x|} \\mathrm{d} x. $$\n",
    "\n",
    "Now identifying the charges as the interpolation nodes, we can inspect what happens at $n \\to \\infty$. For equidistant charges, the charge density in the limit $n \\to \\infty$ is uniform, i.e.\n",
    "\n",
    "$$ \\rho(x) = \\frac{1}{2}, $$\n",
    "\n",
    "(to ensure proper normalization), and (check this yourself)\n",
    "\n",
    "\\begin{align}\n",
    "\\phi(0) &= -1, \\\\\n",
    "\\phi(\\pm 1) & -1 + \\ln 2.\n",
    "\\end{align}\n",
    "\n",
    "This means that $|q_{n+1}|$ is $\\approx e^{(n+1) \\ln 2} = 2^{n+1}$ times larger at the edges of the interval than at the center. The same applies to the basis polynomials $l_k$, which is what causes the issues associated with equispaced nodes.\n",
    "\n",
    "Ideally, then, we seek a charge distribution $\\rho(x)$, such that the generated potential $\\phi(z)$ is uniform (constant) over the interpolation interval. The right distribution turns out to be\n",
    "\n",
    "$$ \\rho(x) = \\frac{1}{\\pi\\sqrt{1 - x^2}}, $$\n",
    "\n",
    "which is the limiting distribution of the nodes $x_j = -\\cos(j \\pi/n)$ for $j = 0, 1, \\ldots, n$, i.e. the Chebyshev nodes. For these set of nodes, the neareast singularity in $f$ may be _arbitrarily close_ to the interval, and we still get geometric convergence.\n",
    "\n",
    "Check it in the following experiment, showing the convergence of (barycentric) interpolation with Chebyshev nodes as the pole of $f$ moves closer to the interval:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "cb3b6669",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "93a98cc14c2248f2ba14f79c24e6a98b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.65, description='x0', max=2.0, min=0.01, step=0.01), Output()), _dom…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<function __main__.draw_convergence(x0)>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compute n+1 Chebyshev polynomial roots\n",
    "chebx = lambda n: np.cos(np.linspace(0, np.pi, n+1))\n",
    "\n",
    "# Equidistant nodes\n",
    "equinodes = lambda n: np.linspace(-1, 1, n)\n",
    "\n",
    "def draw_convergence(x0):\n",
    "    ns = np.arange(1, 200, dtype = int)\n",
    "    equierr = np.array([np.max(np.abs(f(x_star, x0) - baryint(equinodes, n, x0, f))) for n in ns])\n",
    "    chebyerr = np.array([np.max(np.abs(f(x_star, x0) - baryint(chebx, n, x0, f))) for n in ns])\n",
    "    fig, ax = plt.subplots(1, 1, figsize = (6, 6))\n",
    "    ax.semilogy(ns, equierr, label = \"Equidistant\")\n",
    "    ax.semilogy(ns, chebyerr,'--', label = \"Chebyshev\")\n",
    "    ax.set_ylim(1e-16, 1e0)\n",
    "    ax.set_xlim(1,100)\n",
    "    ax.set_ylabel(\"Error\")\n",
    "    ax.set_xlabel(\"$n$\")\n",
    "    ax.legend(loc = \"lower left\")\n",
    "    \n",
    "\n",
    "interact(draw_convergence, x0 = FloatSlider(min = 0.01 , max = 2, step = 0.01, value = 0.65),\\\n",
    "         continuous_update = False) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76e9cef1",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}