{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a398cc49",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Week III - Numerical integration: deterministic methods"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b25ff85",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 1D integrals with trapezoidal and Simpson's rule\n",
    "\n",
    "A common task in numerical analysis is to perform the integration of functions over a certain interval. \n",
    "\n",
    "In the lecture of yesterday, you have seen a number of **deterministic methods** to integrate one-dimensional (1D) functions, in particular your first saw these two:\n",
    "- the **trapezoidal rule**: \n",
    "$$\\int_{x_0}^{x_n} f(x) dx = \\frac{1}{n}\\left[\\frac{f(x_0)}{2} + \\sum_{i=1}^{n-1}f(x_i) + \\frac{f(x_n)}{2}\\right] +\\mathcal{O}(n^{-2}) $$\n",
    "- the **Simpson's rule**\n",
    "$$\\int_{x_0}^{x_n} f(x) dx = \\frac{1}{n}\\left[\\frac{f(x_0)}{3} + \\frac{4f(x_1)}{3}  + \\frac{2f(x_2)}{3} + \\frac{4f(x_3)}{3} + \\cdots + \\frac{2f(x_{n-2})}{3} + \\frac{4f(x_{n-1})}{3} + \\frac{f(x_n)}{3}\\right] + \\mathcal{O}(n^{-5}) $$\n",
    "\n",
    "where the $\\{x_i\\}$ uniformly sample the integration domain from $x_0$ to $x_n$ with a discretization interval equal to $\\Delta x = (x_n-x_0)/n$ .\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b682b59",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "Here below you find a basic implementation for the numerical integration of the following integral\n",
    "$$\\int_0^1 e^x dx = e−1 =1 .718282... $$\n",
    "by using either the trapezoidal or the Simpson's rule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31e48d1a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Loading Numpy module\n",
    "import numpy as np\n",
    "\n",
    "# Here we code integrator\n",
    "def integrate_1D_determ(x_range,n_points,method='trapez'):\n",
    "    '''\n",
    "    This functions integrates exp(x) between x_range[0]\n",
    "    and x_range[1] with a number of intervals equal to n_points\n",
    "    either using the trapezoidal rule or the Simpson' rule\n",
    "    '''    \n",
    "    x = np.linspace(x_range[0],x_range[1], n_points+1,endpoint=True)\n",
    "    dx =  (x_range[1]-x_range[0])/n_points \n",
    "    f = np.exp(x)\n",
    "    #\n",
    "    def trapez():\n",
    "        '''\n",
    "        Trapezoidal rule\n",
    "        '''\n",
    "        #\n",
    "        # Notice how slicing and intrinsic functions allow \n",
    "        # to avoid the explicit \"for\" loops \n",
    "        # with is replaced by a single line of code\n",
    "        #\n",
    "        return dx * (0.5*(f[0]+f[-1]) + np.sum(f[1:-1]))\n",
    "    #\n",
    "    def simpson():\n",
    "        '''\n",
    "        Simpson's rule\n",
    "        '''\n",
    "        #\n",
    "        # Notice how slicing and intrinsic functions allow \n",
    "        # to avoid the explicit \"for\" loops \n",
    "        # with could replaced by a single line of code\n",
    "        # (here a few lines for the sake of clarity)\n",
    "        #\n",
    "        ext = f[0]+f[-1]\n",
    "        ff = f[1:-1]\n",
    "        a = 4.0 * np.sum(ff[::2])\n",
    "        b = 2.0 * np.sum(ff[1::2])\n",
    "        return ( ext + a + b )*(dx/3.0)\n",
    "    \n",
    "    if method=='trapez':\n",
    "        return trapez()\n",
    "    elif method=='simpson':\n",
    "        return simpson()\n",
    "    else:\n",
    "        print('ERROR: your method \"{}\" is not available, please choose either \"trapez\" or \"simpson\"'.format(method))\n",
    "        return None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e620a8e6",
   "metadata": {},
   "source": [
    "-------------\n",
    "\n",
    "### Numerical experiments \n",
    "\n",
    "Now let's play with the code!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa938e53",
   "metadata": {},
   "outputs": [],
   "source": [
    "integrate_1D_determ([0,1],100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f854514",
   "metadata": {},
   "outputs": [],
   "source": [
    "integrate_1D_determ([0,1],100,method='simpson')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67682efc",
   "metadata": {},
   "outputs": [],
   "source": [
    "integrate_1D_determ([0,1],100,method='other')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ef45e84",
   "metadata": {},
   "source": [
    "### Error convergence: trapezoidal VS Simpson's rule\n",
    "\n",
    "Now we investigate the convergence of the absolute error for the numerical integration and compare the two algorithms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d450fd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "n_list = [2**i for i in range(10)]\n",
    "ints_trapez = [integrate_1D_determ([0,1],n_points,method = 'trapez') for n_points in n_list]\n",
    "ints_simpson= [integrate_1D_determ([0,1],n_points,method = 'simpson') for n_points in n_list]\n",
    "ref= np.e - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ad16cb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "[integrate_1D_determ([0,1],n_points,method = 'simpson') for n_points in n_list]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aca73c0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "n_list = [2**i for i in range(12)]\n",
    "ints_trapez = [integrate_1D_determ([0,1],n_points,method = 'trapez') for n_points in n_list]\n",
    "ints_simpson = [integrate_1D_determ([0,1],n_points,method = 'simpson') for n_points in n_list]\n",
    "ref= np.e - 1\n",
    "\n",
    "plt.xlabel('n')\n",
    "plt.ylabel('absolute error')\n",
    "plt.loglog(n_list,[np.abs(ref-res) for res in ints_trapez],'-o', c='blue', label='trapezoidal rule')\n",
    "plt.loglog(n_list,[n**(-2) for n in n_list],'--', c='grey', label='$n^{-2}$')\n",
    "plt.loglog(n_list,[np.abs(ref-res) for res in ints_simpson],'-o', c='red', label=\"Simpson's rule\")\n",
    "plt.loglog(n_list,[n**(-4) for n in n_list],'--', c='k', label='$n^{-4}$')\n",
    "\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08d57f95",
   "metadata": {},
   "source": [
    "### Numerical integration in SciPy\n",
    "\n",
    "The SciPy package offers a powerful interface to numerical integration libraries, including the Fortran library QUADPACK."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7af21073",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.integrate as integrate\n",
    "\n",
    "# Here we precompute the function on a discretized grid\n",
    "n_points = 100\n",
    "x = np.linspace(0,1, n_points+1,endpoint=True)\n",
    "y = np.exp(x)\n",
    "integrate.simpson(y,x)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b6dd9f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we give the function object as input as well as the integration domain\n",
    "integrate.quad(np.exp,0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8febe4e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "integrate.quad?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c280085",
   "metadata": {},
   "source": [
    "Let us compare the implementations with respect to the convergence of absolute errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41471b1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "n_list = [2**i for i in range(12)]\n",
    "ints_trapez = [integrate_1D_determ([0,1],n_points,method = 'trapez') for n_points in n_list]\n",
    "ints_simpson = [integrate_1D_determ([0,1],n_points,method = 'simpson') for n_points in n_list]\n",
    "ints_simpson_scipy = []\n",
    "for n in n_list:\n",
    "    x = np.linspace(0,1, n+1,endpoint=True)\n",
    "    y = np.exp(x)\n",
    "    ints_simpson_scipy.append(integrate.simpson(y,x))\n",
    "ref= np.e - 1\n",
    "\n",
    "plt.xlabel('n')\n",
    "plt.ylabel('absolute error')\n",
    "plt.loglog(n_list,[np.abs(ref-res) for res in ints_trapez],'-o', c='blue', label='trapezoidal rule')\n",
    "plt.loglog(n_list,[n**(-2) for n in n_list],'--', c='grey', label='$n^{-2}$')\n",
    "plt.loglog(n_list,[np.abs(ref-res) for res in ints_simpson],'-o', c='red', label=\"Simpson's rule\")\n",
    "plt.loglog(n_list,[n**(-4) for n in n_list],'--', c='k', label='$n^{-4}$')\n",
    "plt.loglog(n_list,[np.abs(ref-res) for res in ints_simpson_scipy],'-o', c='green', label=\"Simpson's rule (scipy)\")\n",
    "\n",
    "plt.axhline(y=integrate.quad(np.exp,0,1)[1],label='scipy \"quad\" estimated error')\n",
    "plt.axhline(y = np.abs(integrate.quad(np.exp,0,1)[0]-ref),c = 'orange',label='scipy \"quad\" true error')\n",
    "\n",
    "plt.legend()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}