Series_Expansion_Approximations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Series Expansion Approximations\n",
    "\n",
    "This notebook lets you play around with different ways of approximating analytic functions in a given interval.  The three ways are\n",
    "\n",
    "1. Legendre series\n",
    "2. Taylor series\n",
    "3. Fourier series\n",
    "    \n",
    "You have previously seen Taylor series and Fourier series in your first-year courses, but the Legendre polynomials are new in this course on differential equations.\n",
    "\n",
    "Taylor series are defined by\n",
    "\\begin{equation}\n",
    "    T_f(x) = \\sum_{k=0}^\\infty \\frac{f^{(k)}(a)}{k!} {(x-a)}^k\n",
    "    \\quad\\text{where}\\quad\n",
    "    f^{(k)}(a) = \\left.\\frac{\\mathrm{d}^k f(x)}{\\mathrm{d}x^k}\\right\\lvert_{x=a}\n",
    "\\end{equation}\n",
    "and Fourier series on the interval $[-1, 1)$ are defined by\n",
    "\\begin{equation}\n",
    "    F_f(x) = \\frac12 a_{f,0} + \\sum_{k=1}^\\infty \\biggl(a_{f,k}\\cos(n\\pi x) + b_{f,k}\\sin(n\\pi x)\\biggr),\n",
    "\\end{equation}\n",
    "where the coefficients are\n",
    "\\begin{equation}\n",
    "    a_{f,k} = \\int_{-1}^1 f(x) \\cos(k \\pi x)\\,\\mathrm{d}x \\quad\\text{and}\\quad b_{f,k} = \\int_{-1}^1 f(x) \\sin(k \\pi x)\\,\\mathrm{d}x.\n",
    "\\end{equation}\n",
    "\n",
    "The Legendre series is a similar approximation to the Fourier series, but using the Legendre polynomials $P_n$ instead of the trigonometric basis used in Fourier.  This is a series of the form\n",
    "\\begin{equation}\n",
    "    L_f(x) = \\sum_{k=0}^\\infty \\ell_k P_k(x),\n",
    "\\end{equation}\n",
    "where you have considered how to calculate the coefficients $\\ell_k$ elsewhere in this seminar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialisation code\n",
    "\n",
    "These few cells produce some functions for finding (to a good approximation) the various series that we're interested in.  You don't need to understand the details of how any of this works, but feel free to ask.  The Legendre series and Fourier series calculations are done in the same way they are described above, but the Taylor series ones are not.\n",
    "\n",
    "**You need to run all the cells in this section for anything to work.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "_LXe4fZX8Kpm"
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from matplotlib import pyplot\n",
    "import matplotlib\n",
    "import scipy.special\n",
    "import scipy.integrate\n",
    "import numpy as np\n",
    "import numpy.polynomial\n",
    "import ipywidgets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "YLFFM5OQ8KqG"
   },
   "outputs": [],
   "source": [
    "def legendre_series(f, n):\n",
    "    r\"\"\"\n",
    "    Calculate the terms of the Legendre series expansion of the function\n",
    "    ..math:`f(x)` with the first ..math:`n_terms` terms.  This will be the\n",
    "    terms up to but _excluding_ the coefficient of ..math:`P_n(x)`.\n",
    "\n",
    "    The resultant object can be called like a function to return the value of\n",
    "    the approximation at values of ..math:`x`.\n",
    "    \"\"\"\n",
    "    if n < 1:\n",
    "        raise ValueError(\"'n' must be at least 1.\")\n",
    "    def integrand(x, k):\n",
    "        return scipy.special.eval_legendre(k, x)*f(x)\n",
    "    # Approximate the inner product integral for each of the polynomials,\n",
    "    # including the normalisation factor.  `scipy.integrate.quad` performs\n",
    "    # numerical integration (also called 'quadrature') until a particular\n",
    "    # precision goal is reached.\n",
    "    return np.polynomial.legendre.Legendre(np.array([\n",
    "        scipy.integrate.quad(integrand, -1, 1, args=(k,))[0] * (k + 0.5)\n",
    "        for k in range(n)\n",
    "    ]))\n",
    "\n",
    "\n",
    "def taylor_coefficient(f, k, a=15):\n",
    "    r\"\"\"\n",
    "    Calculate the ..math:`k`th coefficient in the Taylor expansion of\n",
    "    ..math:`f(x)` around the point ..math:`x_0 = 0`.  The first term is\n",
    "    ..math:`k = 0`, as this is the zeroth-order term.\n",
    "\n",
    "    ``a`` is a precision factor, and should probably just be left as-is.\n",
    "    \"\"\"\n",
    "    if k == 0:\n",
    "        return f(0)\n",
    "    # The standard way of defining Taylor series with derivatives and\n",
    "    # factorials doesn't play nicely with numerical methods.  This method is\n",
    "    # based on contour integration (magic).\n",
    "    scale = np.exp(-a/k)\n",
    "    return np.exp(a)/k * sum(\n",
    "        (-1)**n * np.imag(f(scale * np.exp(1j*np.pi*(0.5-n)/k)))\n",
    "        for n in range(1, k+1)\n",
    "    )\n",
    "\n",
    "def taylor_series(f, n, a=15):\n",
    "    r\"\"\"\n",
    "    Calculate the first ..math:`n` terms of the Taylor series expansion of\n",
    "    ..math:`f(x)` around the point ..math:`x_0 = 0` up to but excluding the\n",
    "    term ..math:`x^n`.\n",
    "\n",
    "    The resultant object can be called like a function to return the value of\n",
    "    the approximation at values of ..math:`x`.\n",
    "    \"\"\"\n",
    "    if n < 1:\n",
    "        raise ValueError(\"'n' must be at least 1.\")\n",
    "    return np.polynomial.Polynomial([\n",
    "        taylor_coefficient(f, k, a)\n",
    "        for k in range(n)\n",
    "    ])\n",
    "\n",
    "\n",
    "class fourier_series:\n",
    "    r\"\"\"\n",
    "    Calculate the first ..math:`n` terms of the Fourier series expansion of\n",
    "    ..math:`f(x)` when mapped to the period ..math:`[-1, 1)`.\n",
    "\n",
    "    The terms are \"numbered\" in the order\n",
    "    ..math::\n",
    "        a_0, b_1, a_1, b_2, a_2, \\dotsc\n",
    "    This is by analogy to Taylor series; the first term is the constant, then\n",
    "    the lowest-order odd term, the next-lowest even term, and so on.\n",
    "\n",
    "    The resultant object can be called like a function to return the value of\n",
    "    the approximation at values of ..math:`x`.\n",
    "    \"\"\"\n",
    "    def __init__(self, f, n):\n",
    "        if n < 1:\n",
    "            raise ValueError(\"'n' must be at least 1.\")\n",
    "        self._n_a = (n + 1) // 2\n",
    "        self._n_b = n - self._n_a\n",
    "        self.a = np.empty((self._n_a,), dtype=np.float64)\n",
    "        # To keep the labelling clear I store the `b[0] = 0` too.\n",
    "        self.b = np.empty((self._n_b + 1,), dtype=np.float64)\n",
    "        self.a[0] = 0.5 * scipy.integrate.quad(f, -1, 1)[0]\n",
    "        self.b[0] = 0\n",
    "        def cosint(x, k): return f(x) * np.cos(k*np.pi*x)\n",
    "        def sinint(x, k): return f(x) * np.sin(k*np.pi*x)\n",
    "        for k in range(1, self._n_a):\n",
    "            self.a[k] = scipy.integrate.quad(cosint, -1, 1, args=(k,))[0]\n",
    "        for k in range(1, self._n_b):\n",
    "            self.b[k] = scipy.integrate.quad(sinint, -1, 1, args=(k,))[0]\n",
    "\n",
    "    def __call__(self, xs):\n",
    "        out = np.zeros_like(xs)\n",
    "        for k in range(self._n_a):\n",
    "            out += self.a[k] * np.cos(k*np.pi * xs)\n",
    "        for k in range(1, self._n_b):\n",
    "            out += self.b[k] * np.sin(k*np.pi * xs)\n",
    "        return out\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "478tRHOy8Kqq"
   },
   "outputs": [],
   "source": [
    "def series_plot(f, terms):\n",
    "    r\"\"\"\n",
    "    Plot a function ``f`` and its Legendre-, Taylor- and Fourier-series\n",
    "    approximations of the given order on the interval ..math:`[-1, 1)`.\n",
    "    \"\"\"\n",
    "    xs = np.linspace(-1, 1, 301)\n",
    "    cm = matplotlib.cm.viridis\n",
    "    _, axes = pyplot.subplots()\n",
    "    fs = np.array([f(x) for x in xs])\n",
    "    axes.plot(xs, fs, label=\"Exact\", color='black', linewidth=2, dashes=(5, 5))\n",
    "    axes.plot(xs, legendre_series(f, terms)(xs), label=\"Legendre\", color=cm(0))\n",
    "    axes.plot(xs, taylor_series(f, terms)(xs), label=\"Taylor\", color=cm(0.4))\n",
    "    axes.plot(xs, fourier_series(f, terms)(xs), label=\"Fourier\", color=cm(0.8))\n",
    "    axes.set_xlim((xs[0] - 0.03*(xs[-1]-xs[0]), xs[-1] + 0.03*(xs[-1]-xs[0])))\n",
    "    # Base the y axis on the actual function rather than all the data, so the\n",
    "    # Taylor series doesn't blow out the entire scale when it's wrong, and it\n",
    "    # stays the same no matter the number of terms used in the plot.\n",
    "    mn, mx = np.min(fs), np.max(fs)\n",
    "    axes.set_ylim((mn - 0.1*(mx-mn), mx + 0.1*(mx-mn)))\n",
    "    axes.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Investigating series behaviour at low orders\n",
    "\n",
    "I have defined some mathematical functions in Python which are hopefully quite clear from the code.  Below them, I have plotted the function, and the series expansions of the function with a controllable `terms` parameter.  Feel free to copy-paste to see what happens when you use your own functions.\n",
    "\n",
    "Think about these questions, and then discuss them with your group.\n",
    "\n",
    "1. Which series would you call the \"most accurate\"?  Why?  Does it depend on the function?\n",
    "2. What sorts of functions are the different expansions best at approximating?  Which are they bad at?\n",
    "3. The Legendre series often seems to \"give up\" in the middle of some shapes at low orders (_e.g._ sinusoids).  Why is this, and are there any things the series is still useful for?\n",
    "4. The Taylor series almost invariably has the largest pointwise error.  Why is this?  What is the Taylor series useful for?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Approximating a polynomial\n",
    "\n",
    "$f(x) = c_0 + c_1x + c_2x^2 + \\dotsb + c_{12}x^{12}$ for arbitrary constants $\\{c_n\\}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7d1ad04b283a4ce09b4ca048dd8287ef",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=8, description='terms', max=15, min=1), Output()), _dom_classes=('widget…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def high_order_polynomial(x):\n",
    "    return np.polynomial.Polynomial([\n",
    "        -0.0372875, 0.674885, 1.34898, -12.652, -7.15369, 57.7268,\n",
    "        8.73373, -104.258, 10.0257, 79.9955, -21.4594, -21.5587, 8.38861,\n",
    "    ])(x)\n",
    "\n",
    "ipywidgets.interact(\n",
    "    lambda terms: series_plot(high_order_polynomial, terms),\n",
    "    terms=(1, 15, 1),\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Approximating a logistic function\n",
    "\n",
    "$f(x) = {\\bigl(1 + e^{-5x}\\bigr)}^{-1}$, which makes a sigmoid (\"s\") shape, asymptotically tending to $0$ and $1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "id": "BT26OANA8Kq_",
    "outputId": "8b44a71b-975f-4560-dd80-f4423b5122b7"
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "026949bcd43d4e28a56f14ad20f1cb15",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=6, description='terms', max=12, min=1), Output()), _dom_classes=('widget…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def logistic(x):\n",
    "    return 1 / (1 + np.exp(-5*x))\n",
    "\n",
    "ipywidgets.interact(\n",
    "    lambda terms: series_plot(logistic, terms),\n",
    "    terms=(1, 12, 1),\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Approximating a Lorentzian function\n",
    "\n",
    "$f(x) = {\\bigl(\\varGamma\\pi + {(x/\\varGamma)}^2\\bigr)}^{-1}$.  This sort of shape is common when talking about linewidths of atomic transitions and lasers, and is a valid probability distribution.  Unlike Gaussians, they can be quite tricky to work with, because no integer statistical moments of them (other than the zeroth order) exist—the series only work here because the domain is limited to $[-1, 1]$ rather than being infinite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "id": "yd5G9xp79BgX"
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d568a1f2fa19474b8312882b9662d534",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=6, description='terms', max=12, min=1), Output()), _dom_classes=('widget…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def lorentzian(x, c=0.2):\n",
    "    return 1 / (c*np.pi + (x/c)**2)\n",
    "\n",
    "ipywidgets.interact(\n",
    "    lambda terms: series_plot(lorentzian, terms),\n",
    "    terms=(1, 12, 1),\n",
    ");"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "Series Expansion Approximations.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}