{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simpson rule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def simpson(a,b,n,f,df3):\n",
    "    h = (b-a)/n\n",
    "    x = np.linspace(a,b,n+1)\n",
    "    y = f(x)\n",
    "    res = 4.0*np.sum(y[1:n:2]) + 2.0*np.sum(y[2:n-1:2]) + y[0] + y[n]\n",
    "    est = -(h**4/180.0)*(df3(b) - df3(a))\n",
    "    return (h/3.0)*res, est"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example\n",
    "$$\n",
    "f(x) = \\exp(x)\\cos(x), \\qquad x \\in [0,\\pi]\n",
    "$$\n",
    "The exact integral is $-\\frac{1}{2}(1+\\exp(\\pi))$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     2    -4.77506763e-01        0.00000    -1.63300203e+00\n",
      "     4    -8.54022966e-02        5.59126    -1.02062627e-01\n",
      "     8    -6.13735917e-03       13.91515    -6.37891419e-03\n",
      "    16    -3.94993112e-04       15.53789    -3.98682137e-04\n",
      "    32    -2.48603363e-05       15.88849    -2.49176335e-05\n",
      "    64    -1.55645818e-06       15.97238    -1.55735210e-06\n",
      "   128    -9.73205481e-08       15.99311    -9.73345060e-08\n",
      "   256    -6.08319262e-09       15.99827    -6.08340663e-09\n",
      "   512    -3.80214971e-10       15.99935    -3.80212914e-10\n",
      "  1024    -2.37658782e-11       15.99836    -2.37633071e-11\n"
     ]
    }
   ],
   "source": [
    "# Function f\n",
    "f = lambda x: np.exp(x)*np.cos(x)\n",
    "# Third derivative of f\n",
    "df3 = lambda x: -2.0*np.exp(x)*(np.cos(x) + np.sin(x))\n",
    "\n",
    "qe = -0.5*(1.0 + np.exp(np.pi)) # Exact integral\n",
    "\n",
    "n,N = 2,10\n",
    "e = np.zeros(N)\n",
    "for i in range(N):\n",
    "    integral,est = simpson(0.0,np.pi,n,f,df3)\n",
    "    e[i] = qe - integral\n",
    "    if i > 0:\n",
    "        print('%6d %18.8e %14.5f %18.8e'%(n,e[i],e[i-1]/e[i],est))\n",
    "    else:\n",
    "        print('%6d %18.8e %14.5f %18.8e'%(n,e[i],0,est))\n",
    "    n = 2*n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The convergence rate is $O(h^4)$ and the error estimate becomes very good as $n$ increases."
   ]
  }
 ],
 "metadata": {
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}