lab01_trap_integration_plots.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# ASTR 21100 \n",
    "\n",
    "# *\"Computational Techniques in Astrophysics\"*\n",
    "\n",
    "## Lab session 1\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": [
    "### Part I. Trapezoidal integration\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."
   ]
  },
  {
   "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-1}) + \\sum\\limits_{i=1}^{N-2}f_i,\\right]$$"
   ]
  },
  {
   "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 the trapezoidal method. Try to carry out computation of the integral using function <a href=\"https://docs.scipy.org/doc/numpy/reference/generated/numpy.sum.html\">numpy.sum</a> from the  numpy library.  \n",
    "\n",
    "    If you are unsure how to write such routine, please let us know. \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-A_{\\rm trap}(h)/A_{\\rm direct}\\vert,$$\n",
    "where $A_{\\rm trap}(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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# use jupyter \"magic\" command to tell it to embed plot into the notebook \n",
    "%matplotlib inline\n",
    "\n",
    "def plot_line(x, y, figsize=6, xlabel=' ', ylabel=' ', col= 'darkslateblue', legend=None, figsave = None):\n",
    "    plt.figure(figsize=(figsize,figsize))\n",
    "    plt.xlabel(xlabel); plt.ylabel(ylabel)\n",
    "    \n",
    "    if legend:\n",
    "        plt.plot(x, y, lw = 1., c=col, label = legend)\n",
    "        plt.legend(frameon=False, loc='lower left')\n",
    "    else:\n",
    "        plt.plot(x, y, lw = 1., c=col)\n",
    "\n",
    "    if figsave:\n",
    "        plt.savefig(figsave, bbox_inches='tight')\n",
    "\n",
    "    plt.show()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "### Making \"publication quality\" plots\n",
    "\n",
    "We will discuss the basic guidelines for plots after you make your plots. \n",
    "\n",
    "Matplotlib allows for very detailed customization of plots. It may be tricky to find customization you want, but options in the routine <tt>plot_pretty</tt> above allow to control the main features of the plot. "
   ]
  }
 ],
 "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
}