lab01b_trap_integration.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# ASTR 21100 \n",
    "\n",
    "# *\"Computational Techniques in Astrophysics\"*\n",
    "\n",
    "## Lab session 1 (Thursday, Jan 16, 2020)\n",
    "\n",
    "### Instructor: Andrey Kravtsov\n",
    "\n",
    "### office: ERC 415; email: kravtsov@uchicago.edu\n",
    "### office hours: Tue, 10:30am-noon (unless noted otherwise)\n",
    "\n",
    "### Teaching Assistants: \n",
    "### Dimitrios Tanoglidis (dtanoglidis@uchicaago.edu)\n",
    "### Georgios Zakharegkas (gzakharegkas@uchicago.edu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Background: from trapezoidal to Simpson integration scheme\n",
    "\n",
    "One of the simplest methods of numerical integration is *\"trapezoidal integration.\"* It can be used to get reasonably accurate numerical integral estimate and easiest to remember. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Specifically, to evaluate integral of $f(x)$ (area under $f(x)$ for $x\\in[a,b]$):\n",
    "\n",
    "Split the integration interval into $N$ equal-size sub-intervals of size $h=(b-a)/N=x_{i+1}-x_i$. For convenience, let's denote $f_i = f(x_i)$, $f_{i+1}=f(x_{i+1})$, etc. There will be $N+1$ points at which the integrated function needs to be evaluated, including interval ends $a$ and $b$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Approximate area under $f(x)$ in each interval, $A_i$, can then be approximated by the area of the trapezoid formed by vertices $(x_i, 0)$, $(x_{i+1},0)$, $(x_i, f_i)$, $(x_{i+1}, f_{i+1})$, $T_i$:\n",
    "\n",
    "$$T_i = \\frac{1}{2}(f_i + f_{i+1})\\, h.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The total area under $f(x)$ in the interval $[a,b]$ is then:\n",
    "\n",
    "$$A = \\int\\limits_a^b f(x)\\, dx = \\sum\\limits_{i=0}^{N-1} A_i \\approx \\sum\\limits_{i=0}^{N-1} T_i= \\frac{h}{2}\\sum\\limits_{i=0}^{N-1}(f_i + f_{i+1}),$$\n",
    "\n",
    "which can be recast as\n",
    "\n",
    "$$A \\approx h\\left[\\frac{1}{2}(f_0 + f_{N}) + \\sum\\limits_{i=1}^{N-1}f_i,\\right]$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In class on Wednesday, January 15, I showed how trapezoidal scheme can be used to get a scheme of \"higher order\" - that is a scheme in which fractional error scales as higher power of step size, $h$. \n",
    "\n",
    "Namely, if I denote trapezoidal scheme for a step size $h$ as\n",
    "\n",
    "$$T(h) = \\frac{h}{2}\\sum\\limits_{i=0}^{N-1}(f_i + f_{i+1}) = h\\left[\\frac{1}{2}(f_0 + f_{N}) + \\sum\\limits_{i=1}^{N-1}f_i,\\right],$$\n",
    "\n",
    "I showed in class (see <a href=\"https://github.com/a-kravtsov/a211w20/blob/master/02_integration_approximation.ipynb\">02_integration_approximation.ipynb</a> notebook in the github repository) that the following construct\n",
    "\n",
    "$$R_1(h) = \\frac{1}{3}\\left[4T(h/2) - T(h)\\right]$$\n",
    "\n",
    "has fractional error that scales as $\\epsilon_{R_1}\\propto h^4$ instead of $\\epsilon_T\\propto h^2$ of the trapezoidal scheme. \n",
    "\n",
    "Usually, it is said that the scheme $R_1(h)$ has a higher order of accuracy  than the trapezoidal scheme. This means that it will give a more accurate result for a given step size $h$. $R_1(h)$ is actually the integration scheme known as the *Simpson rule*. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A more efficient calculation of $R_1(h)$\n",
    "\n",
    "Computation of $R_1(h)$ involves computation of $T(h)$ and $T(h/2)$ for the same interval and function. \n",
    "\n",
    "However, computation of $T(h/2)$ replicates computation of function values at all of the nodes that were used to compute $T(h)$. In other words, computation of $T(h)$ is *embedded* in computation of $T(h/2)$. It is thus inefficient/redundant to compute $T(h)$ and $T(h/2)$ separately via two calls, as I have done in the <a href=\"https://github.com/a-kravtsov/a211w20/blob/master/02_integration_approximation.ipynb\">02_integration_approximation.ipynb</a> notebook. \n",
    "\n",
    "If we denote $h_j = (b-a)/2^j$ and trapezoidal scheme with step $h_j$ as $T_j$, then it is easy to show that trapezoidal scheme for $j=1, 2, \\ldots$ can be written as\n",
    "\n",
    "$$T_j = \\frac{1}{2} T_{j-1} + h_j\\sum\\limits_{i=1}^{2^{j-1}}f(a+(2i-1)h_j).$$\n",
    "\n",
    "This expression allows to reuse the previous computations used to estimate $T_{j-1}$ in computing $T_j$. Taking this into account, the computation of $R_1(h)$ described above, can be written as \n",
    "\n",
    "$$R_1(h) = \\frac{1}{3}\\left[4T(h/2)-T(h)\\right]=\\frac{1}{3}\\left[T(h)+2h\\sum\\limits_{i=1}^{2^{j-1}}f\\left(a+(2i-1)\\frac{h}{2}\\right)\\right],$$\n",
    "\n",
    "where $j=\\log_2[2(b-a)/h]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Tasks for the lab\n",
    "\n",
    "* Write routine for trapezoidal integration that takes input of the function to be integrated, end points $a$ and $b$ ($b>a$), defining the integration range, and step size $h\\leq b-a$, and returns the estimate of the integral obtained using $R_1(h)$ using the efficient formula above. \n",
    "\n",
    "First, implement the routine in python the way with which you are most familiar. When you get your routine working and after you test it, try to write another version that uses numpy vector operations and avoiding any explicit loops in computing $R_1(h)$. You can use the routine <tt>trapzd</tt> in the <a href=\"https://github.com/a-kravtsov/a211w20/blob/master/02_integration_approximation.ipynb\">02_integration_approximation.ipynb</a> as an example.  It may be helpful to use numpy functions <a href=\"https://docs.scipy.org/doc/numpy/reference/generated/numpy.arange.html\">numpy.arange</a> and <a href=\"https://docs.scipy.org/doc/numpy/reference/generated/numpy.sum.html\">numpy.sum</a>. The former is useful to set up a vector of $x_i$ values at which to evaluate $f$ and the latter can be used for vectorized summation of the result.  \n",
    "\n",
    "\n",
    "* Test your routine by comparing its result to the result of direct calculation of the integral for a function, for which the \"closed form\" expression for the integral exists - that is there is an analytic expression for the integral, such as $\\int_0^1 e^x dx$, $\\int_0^\\pi\\sin x dx$, etc.\n",
    "\n",
    "* *The fractional error* of your integral estimate is computed as\n",
    "$$\\epsilon = \\vert 1-R_1(h)/A_{\\rm direct}\\vert,$$\n",
    "where $R_1(h)$ is the estimate of the integral using trapezoidal method with step $h$ and $A_{\\rm direct}$ is the result of direct calculation of the integral. \n",
    "\n",
    "    Compute the fractional error of the test integral for a grid of values of $h$ and record it in an array or list. Make sure your $h$-grid spans sufficiently wide range of values, so that fractional error reaches values at least as low as $\\approx 10^{-16}$. \n",
    "\n",
    "\n",
    "* Plot the $\\log_{10}$ of the fractional error $\\epsilon$ as a function of $\\log_{10}$ of step size $h$ using routine <tt>plot_line</tt> below. Examine the plot and try answer the following questions:\n",
    "\n",
    "    1) does the fractional error decrease monotonically with decreasing $h$ for all values of $h$? If not, what do you think prevents it from decreasing? \n",
    "    \n",
    "    2) In the range of $h$ values where fractional error is monotonically decreasing with decreasing $h$, can you figure out the functional form of the dependence of $\\epsilon$ on $h$ from your plot?\n",
    "    \n",
    "    3) Suppose instead of the test integral for which the correct answer can be calculated using \"closed form\" expression of the integral, we are dealing with an integral for which not such closed form expression exists. Can you think of a way to evaluate the accuracy of your estimate of the integral in this case? \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Matplotlib\n",
    "\n",
    "Matplotlib is a package for (mostly) 2D plots built upon the Numpy and Scipy libraries. It was conceived by <a href=\"https://en.wikipedia.org/wiki/John_D._Hunter\">John Hunter</a> at U.Chicago in 2002, developed by him and others over the subsequent decade into a full-fledged library. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# import matplotlib's PyPlot library denoting it as plt for brevity\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def plot_pretty(dpi=175,fontsize=9):\n",
    "    # import pyplot and set some parameters to make plots prettier\n",
    "    plt.rc(\"savefig\", dpi=dpi)\n",
    "    plt.rc(\"figure\", dpi=dpi)\n",
    "    plt.rc('font', size=fontsize)\n",
    "    plt.rc('xtick', direction='in') \n",
    "    plt.rc('ytick', direction='in')\n",
    "    plt.rc('xtick.major', pad=5) \n",
    "    plt.rc('xtick.minor', pad=5)\n",
    "    plt.rc('ytick.major', pad=5) \n",
    "    plt.rc('ytick.minor', pad=5)\n",
    "    plt.rc('lines', dotted_pattern = [0.5, 1.1])\n",
    "\n",
    "    return\n",
    "\n",
    "plot_pretty(dpi=250, fontsize=9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### LaTeX\n",
    "\n",
    "Although not strictly required, I highly recommend installing <a href=\"https://www.latex-project.org/get/\">a LaTeX distribution</a> on your laptop. I use LaTeX commands to format plot labels in matplotlib. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "#if you don't have LaTeX installed on your laptop and this statement \n",
    "# generates error, comment it out\n",
    "plt.rc('text', usetex=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "plot_pretty(dpi=150, fontsize=12)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}