sine-series.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Activity 10: Sine Series Problem \n",
    "\n",
    "* *Computational Physics*: Ch 2.5\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem: Summing Series: sin(x)\n",
    "\n",
    "Evaluate the $\\sin$ function from its series representation\n",
    "$$\n",
    "\\sin x = x - \\frac{x^3}{3!} + \\frac{x^5}{5!} - \\frac{x^7}{7!} + \\dots\n",
    "$$\n",
    "\n",
    "### Explicit sum\n",
    "A naive algorithm is to sum the series up to the $N$th term:\n",
    "$$\n",
    "\\sin x \\approx \\sum_{n=1}^N \\frac{(-1)^{n-1} x^{2n-1}}{(2n - 1)!}\n",
    "$$\n",
    "\n",
    "Problems:\n",
    "\n",
    "- How to decide when to stop summing?\n",
    "- Division of large terms (overflows!)\n",
    "- Powers and factorials are very expensive to compute.\n",
    "\n",
    "### Recursive sum\n",
    "Better approach: Build up series terms $a_n$ using previous term $a_{n-1}$ through a recursion:\n",
    "\n",
    "\\begin{align}\n",
    "a_n &= a_{n-1} \\times q_n\\\\\n",
    "a_n &= \\frac{(-1)^{n-1} x^{2n-1}}{(2n - 1)!} = \\frac{(-1)^{n-2} x^{2n-3}}{(2n - 3)!} \\frac{-x^2}{(2n - 1)(2n - 2)}\\\\\n",
    "a_n & = a_{n-1} \\frac{-x^2}{(2n - 1)(2n - 2)}\n",
    "\\end{align}\n",
    "\n",
    "Accuracy of this approach? Not clear in absolute terms but we can make the assumption that the error is approximately the last term in the sum, $a_N$. Hence we can strive to make the relative error smaller than the machine precision\n",
    "$$\n",
    "\\left| \\frac{a_N}{\\sum_{n=1}^N a_n} \\right| < \\epsilon_m\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task\n",
    "\n",
    "In file `series.py` implement a function `sin_recursive(x, eps=1e-16)` that \n",
    "* computes $\\sin(x)$ to relative accuracy $\\epsilon$ using the **Recursive sum** approach\n",
    "* returns the value of $\\sin(x)$ and the number of iterations $N$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluation and plotting "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import your `series.py` as a module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import series"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: If you changed code in `series` you will need to **reload the module**: execute the following cell whenever you changed `series.py`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from importlib import reload\n",
    "reload(series)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test the implementation against the \"exact\" numpy function `np.sin()`.\n",
    "\n",
    "Report\n",
    "1. `x`\n",
    "2. maximum `n`\n",
    "3. `sin_series(x)`\n",
    "4. relative error `abs(sin_series(x) - sin(x))/abs(sin(x))`\n",
    "\n",
    "Plot against `x` the quantities above and also `sin(x)`.\n",
    "* `x` - `sin_series(x)`\n",
    "* `x` - `sin(x)`\n",
    "* `x` - max `n`\n",
    "* `x` - relative error (semilogy plot!)\n",
    "\n",
    "For a range of numbers look at `Xsmall` and `Xlarge`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xsmall = np.pi*np.arange(-2, 2.05, 0.05)\n",
    "Xlarge = np.pi*np.arange(-50, 50.05, 0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.pi*np.arange(-2, 2.2, 0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Testing code\n",
    "\n",
    "Implementation of the test against the numpy reference implementation `np.sin`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_sin(x):\n",
    "    y, nmax = series.sin_recursive(x)\n",
    "    y0 = np.sin(x)\n",
    "    delta = y - y0\n",
    "    if y0 != 0:\n",
    "        relative_error = delta/y0\n",
    "    else:\n",
    "        relative_error = 0\n",
    "    return x, nmax, y, relative_error "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plotting code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_plot_sine(X, filename=\"sine_error.pdf\"):\n",
    "    results = np.array([test_sin(x) for x in X])\n",
    "    \n",
    "    fig = plt.figure(figsize=(8, 10))\n",
    "    \n",
    "    ax1 = fig.add_subplot(3,1,1)\n",
    "    ax1.plot(results[:, 0], results[:, 2], 'k-', lw=1, label=\"series\")\n",
    "    ax1.plot(results[:, 0], np.sin(results[:, 0]), 'g--', lw=2, label=\"sin x\")\n",
    "    ax1.legend(loc=\"best\")\n",
    "\n",
    "    ax2 = fig.add_subplot(3,1,2)\n",
    "    ax2.plot(results[:, 0], results[:, 1], label=\"max n\")\n",
    "    ax2.legend(loc=\"best\")\n",
    "\n",
    "    ax3 = fig.add_subplot(3,1,3)\n",
    "    ax3.semilogy(results[:, 0], results[:, 3], label=\"rel.error\")\n",
    "    ax3.legend(loc=\"best\")\n",
    "    \n",
    "    fig.suptitle(\"sine series approximation\")\n",
    "    \n",
    "    fig.savefig(filename)\n",
    "    print(\"saved to file {0}\".format(filename))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Analyze recursive implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now test the two ranges of numbers and write to different files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_plot_sine(Xsmall, filename=\"sine_error_Xsmall.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### BONUS: explicit sum \n",
    "\n",
    "Implement the **Explicit sum** as `series.sin_explicit(x, eps=1e-16)` and compare it against `series.sin_recursive()`.\n",
    "\n",
    "(Note: no tests are provided)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}