ch08-figures.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sympy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integral"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x, x0 = sympy.symbols(\"x, x_0\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = (x+0.5) ** 3 - 3.5 * (x+0.5) **2 + x  + 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f.subs(x,x0) + f.diff(x).subs(x,x0) * (x - x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f_func = sympy.lambdify(x, f, 'numpy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def arrowify(fig, ax):\n",
    "    xmin, xmax = ax.get_xlim()\n",
    "    ymin, ymax = ax.get_ylim()\n",
    "\n",
    "    # removing the default axis on all sides:\n",
    "    for side in ['bottom','right','top','left']:\n",
    "        ax.spines[side].set_visible(False)\n",
    "\n",
    "    # removing the axis labels and ticks\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.xaxis.set_ticks_position('none')\n",
    "    ax.yaxis.set_ticks_position('none')\n",
    "\n",
    "    # wider figure for demonstration\n",
    "    #fig.set_size_inches(4,2.2)\n",
    "\n",
    "    # get width and height of axes object to compute\n",
    "    # matching arrowhead length and width\n",
    "    dps = fig.dpi_scale_trans.inverted()\n",
    "    bbox = ax.get_window_extent().transformed(dps)\n",
    "    width, height = bbox.width, bbox.height\n",
    "\n",
    "    # manual arrowhead width and length\n",
    "    hw = 1./25.*(ymax-ymin)\n",
    "    hl = 1./25.*(xmax-xmin)\n",
    "    lw = 0.5 # axis line width\n",
    "    ohg = 0.3 # arrow overhang\n",
    "\n",
    "    # compute matching arrowhead length and width\n",
    "    yhw = hw/(ymax-ymin)*(xmax-xmin)* height/width\n",
    "    yhl = hl/(xmax-xmin)*(ymax-ymin)* width/height\n",
    "\n",
    "    # draw x and y axis\n",
    "    ax.arrow(xmin, 0, xmax-xmin, 0., fc='k', ec='k', lw = lw,\n",
    "             head_width=hw, head_length=hl, overhang=ohg,\n",
    "             length_includes_head=True, clip_on=False)\n",
    "\n",
    "    ax.arrow(0, ymin, 0., ymax-ymin, fc='k', ec='k', lw = lw,\n",
    "             head_width=yhw, head_length=yhl, overhang=ohg,\n",
    "             length_includes_head=True, clip_on=False)\n",
    "\n",
    "    \n",
    "    ax.text(xmax*0.95, ymin*0.25, r'$x$', fontsize=22)\n",
    "    ax.text(xmin*0.35, ymax*0.9, r'$f(x)$', fontsize=22)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(8,4))\n",
    "\n",
    "xvec = np.linspace(-1.75, 3.0, 100)\n",
    "ax.plot(xvec, f_func(xvec), 'k', lw=2)\n",
    "\n",
    "xvec = np.linspace(-0.8, 0.85, 100)\n",
    "ax.fill_between(xvec, f_func(xvec), color='lightgreen', alpha=0.9)\n",
    "xvec = np.linspace(0.85, 2.31, 100)\n",
    "ax.fill_between(xvec, f_func(xvec), color='red', alpha=0.5)\n",
    "xvec = np.linspace(2.31, 2.6, 100)\n",
    "ax.fill_between(xvec, f_func(xvec), color='lightgreen', alpha=0.99)\n",
    "\n",
    "ax.text(0.6, 3.5, r\"$\\int_a^b\\!f(x){\\rm d}x$\", fontsize=22)\n",
    "ax.text(-0.88, -0.85, r\"$a$\", fontsize=18)\n",
    "ax.text(2.55, -0.85, r\"$b$\", fontsize=18)\n",
    "ax.axis('tight')\n",
    "\n",
    "arrowify(fig, ax)\n",
    "fig.savefig(\"ch8-illustration-integral.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quadrature rules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from numpy import polynomial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import integrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy import interpolate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = 0\n",
    "b = 1.0\n",
    "\n",
    "def f(x):\n",
    "    return np.exp(-x**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = -1.0\n",
    "b = 1.0\n",
    "\n",
    "def f(x):\n",
    "    return 3 + x + x**2 + x**3 + x**4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = np.linspace(a, b, 100)\n",
    "xx = np.linspace(a-0.2, b+0.2, 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16, 4))\n",
    "\n",
    "npoints = 2\n",
    "npoints = 5\n",
    "\n",
    "X = np.linspace(a, b, npoints)\n",
    "\n",
    "ax1.plot(xx, f(xx), lw=1, color='k')\n",
    "ax1.fill_between(x, f(x), color='lightgreen', alpha=0.9)\n",
    "\n",
    "for n in range(len(X) - 1):\n",
    "    f_mid = f(X[n:n+2].mean())\n",
    "    ax1.plot([X[n], X[n]], [0, f_mid], 'b')\n",
    "    ax1.plot([X[n+1], X[n+1]], [0, f_mid], 'b')\n",
    "    ax1.plot(X[n:n+2], [f_mid] * 2, 'b')\n",
    "    ax1.plot(X[n:n+2].mean(), f_mid, 'xk')\n",
    "    \n",
    "i = (b-a)*f_mid\n",
    "#ax1.text(-1, 7, r'$\\int_{-1}^{\\,1} f(x)dx \\approx %f$' % i, fontsize=18)\n",
    "ax1.plot(X, f(X), 'ro')\n",
    "ax1.set_xlim(xx.min(), xx.max())\n",
    "ax1.set_title('Mid-point rule')\n",
    "ax1.set_xticks([-1, 0, 1])\n",
    "ax1.set_xlabel(r'$x$', fontsize=18)\n",
    "ax1.set_ylabel(r'$f(x)$', fontsize=18)\n",
    "\n",
    "names = [\"Trapezoid rule\", \"Simpson's rule\"]\n",
    "for idx, ax in enumerate([ax2, ax3]):\n",
    "    ax.plot(xx, f(xx), lw=1, color='k')\n",
    "    ax.fill_between(x, f(x), color='lightgreen', alpha=0.9)\n",
    "\n",
    "    i = 0\n",
    "    for n in range(len(X) - 1):\n",
    "        XX = np.linspace(X[n], X[n+1], idx+2)\n",
    "\n",
    "        f_interp = polynomial.Polynomial.fit(XX, f(XX), len(XX)-1)    \n",
    "        ax.plot([X[n], X[n]], [0, f(X[n])], 'b')\n",
    "        ax.plot([X[n+1], X[n+1]], [0, f(X[n+1])], 'b')\n",
    "        XXX = np.linspace(X[n], X[n+1], 50)\n",
    "        ax.plot(XXX, f_interp(XXX), 'b')\n",
    "        \n",
    "        F = f_interp.integ()\n",
    "        i += F(X[n+1]) - F(X[n])\n",
    "    ax.text(-1, 5.5, r'$\\int_a^{\\,b} f(x)dx \\approx %f$' % i, fontsize=18)\n",
    "    ax.plot(X, f(X), 'ro')\n",
    "    ax.set_xlabel(r'$x$', fontsize=18)\n",
    "    ax.set_ylabel(r'$f(x)$', fontsize=18)\n",
    "    ax.set_xlim(xx.min(), xx.max())\n",
    "    ax.set_title(names[idx])\n",
    "    ax.set_xticks([-1, 0, 1])\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.savefig('ch8-quadrature-rules-%d.pdf' % npoints)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# mid-point rule\n",
    "(b-a) * f((b+a)/2.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# trapezoid rule\n",
    "(b-a)/2.0 * (f(a) + f(b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# simpsons rule\n",
    "(b-a)/6.0 * (f(a) + 4 * f((a+b)/2.0) + f(b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# exact result\n",
    "integrate.quad(f, a, b)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "integrate.trapz(f(X), X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "integrate.simps(f(X), X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "integrate.quadrature(f, a, b)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "integrate.romberg(f, a, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "integrate.newton_cotes(2)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.fill_between(x, f(x), alpha=0.25)\n",
    "ax.plot(X, f(X), 'ro')\n",
    "\n",
    "for n in range(len(X) - 1):\n",
    "    f_mid = f(X[n:n+2].mean())\n",
    "    ax.plot([X[n], X[n]], [0, f_mid], 'k')\n",
    "    ax.plot([X[n+1], X[n+1]], [0, f_mid], 'k')\n",
    "    ax.plot(X[n:n+2], [f_mid] * 2, 'k')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.fill_between(x, f(x), alpha=0.25)\n",
    "ax.plot(X, f(X), 'ro')\n",
    "\n",
    "for n in range(len(X) - 1):\n",
    "    f_mid = f(X[n:n+2].mean())\n",
    "    ax.plot([X[n], X[n]], [0, f(X[n])], 'k')\n",
    "    ax.plot([X[n+1], X[n+1]], [0, f(X[n+1])], 'k')\n",
    "    ax.plot(X[n:n+2], [f(X[n]), f(X[n+1])], 'k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.6 (NPM ed2)",
   "language": "python",
   "name": "py3.6-npm-e2"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}