phys512/integration_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cd079071",
   "metadata": {},
   "source": [
    "# Integration"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0431bce",
   "metadata": {},
   "source": [
    "## Newton-Cotes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "18ed0344",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3 0.34075853066724404 0.0519405510314801 0.0022798774922103693\n",
      "5 0.1834653418221377 0.012884199027224597 0.00013458497419382986\n",
      "9 0.0949599423108507 0.003214828113830337 8.295523967971619e-06\n",
      "17 0.04828406569741284 0.0008033195149277361 5.166847063531321e-07\n",
      "33 0.024342886926189244 0.00020080567998093102 3.2265001115305836e-08\n",
      "65 0.012221646395186303 5.0199907898895724e-05 2.0161288194486815e-09\n",
      "129 0.006123373269068644 1.2549882473678053e-05 1.2600120946615334e-10\n",
      "257 0.003064824111058906 3.137464712255067e-06 7.87503395827116e-12\n",
      "513 0.0015331964220766103 7.843658088590999e-07 4.922728891187944e-13\n",
      "1025 0.0007667943025135848 1.960914292054028e-07 3.086420008457935e-14\n",
      "2049 0.00038344617411567583 4.902285577479404e-08 1.9984014443252818e-15\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGiCAYAAAAfnjf+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACGK0lEQVR4nOzdd3QUZd/G8e/upnfSewMSegIBQq/B0EITBCvV3hELFsAKShERHrsCVkSlSIdApPdekpCQQEjvyW6STdl9/4jmeXiJChqYlN/nnD3HnZ3MXhuRvZx75r5VRqPRiBBCCCFEA6dWOoAQQgghRF2QUiOEEEKIRkFKjRBCCCEaBSk1QgghhGgUpNQIIYQQolGQUiOEEEKIRkFKjRBCCCEaBSk1QgghhGgUTJQOcLsYDAbS0tKwtbVFpVIpHUcIIYQQN8BoNFJcXIynpydq9V+fi2kypSYtLQ0fHx+lYwghhBDiH0hJScHb2/sv92kypcbW1hao/qXY2dkpnEYIIYQQN6KoqAgfH5+a7/G/0qBKzYYNG3juuecwGAy8+OKLTJs27YZ/9o8hJzs7Oyk1QgghRANzI5eONJhSU1lZyfTp09m1axf29vaEhYUxevRonJyclI4mhBBCiHqgwdz9dPjwYdq2bYuXlxc2NjYMGTKEbdu2KR1LCCGEEPXEbSs1u3fvJioqCk9PT1QqFWvXrr1un2XLluHv74+FhQXh4eEcPny45rW0tDS8vLxqnnt5eZGamno7ogshhBCiAbhtpUan0xESEsKyZctqfX3VqlVMnz6d2bNnc/z4cUJCQoiMjCQrK+sfvZ9er6eoqOiahxBCCCEar9tWaoYMGcJbb73F6NGja3190aJFPPjgg0yePJk2bdrw8ccfY2VlxZdffgmAp6fnNWdmUlNT8fT0/NP3mzt3Lvb29jUPuZ1bCCGEaNzqxTU15eXlHDt2jIiIiJptarWaiIgIDhw4AEDXrl05e/YsqampaLVaNm/eTGRk5J8ec+bMmRQWFtY8UlJSbvnnEEIIIYRy6sXdTzk5OVRVVeHm5nbNdjc3N2JjYwEwMTFh4cKF9O/fH4PBwAsvvPCXdz6Zm5tjbm5+S3MLIYQQov6oF6XmRo0YMYIRI0YoHUMIIYQQ9VC9GH5ydnZGo9GQmZl5zfbMzEzc3d0VSiWEEEKIhqRelBozMzPCwsKIjo6u2WYwGIiOjqZ79+4KJhNCCCFEQ3Hbhp+0Wi0JCQk1z5OSkjh58iSOjo74+voyffp0Jk6cSOfOnenatSuLFy9Gp9MxefLk2xVRCCGEEA3YbSs1R48epX///jXPp0+fDsDEiRNZvnw548ePJzs7m1mzZpGRkUFoaChbtmy57uJhIYQQQojaqIxGo1HpELdDUVER9vb2FBYW3pIFLau0WjQ2NnV+XCGEEKIpu5nv73pxTU1DZ6yo4GLPXiQMGEhlTs5/t1dVKZhKCCGEaFqk1PxLVQYjD81bi1GvpzS/kH05VWQVlWE0Gsmav4CEAQPJX71a6ZhCCCFEo9eg5qmpjy7n6tius2LfsDfx0OWRuPwoAI7WZrwTcwC/tDSOXCnC/WoBLV1t0WRlcOWBB7Ds2BHPBfNRqVQKfwIhhBCicZBS8y+52lnw2QOdicso4kJGMaQXkZSjI09XznMh99HcP5UrKdYULN2HWgWjiuJ4MC2NXI0F585n0trDDi8HSzLfeJ2qvHycHnwQy/btlP5YQgghRIMjpeZfsjE3YVAbNwa1+e9dWmUVVVzM1BKbUURsRjG2GUVcSC8mT1fOZgs/Eno+gomhiuNfHwPA2kzDlxu3Yqcr4EjnO/C29STY3RaT82fJW74c6149aXbXXUp9RCGEEKJBkFJzC1iYamjvbU97b/uabUajkWytnriMYmLTi7mQUURZejEJWVp0+kreCR1Pi4KrbDitR3++ehHPqSl7GHtsG3FZWkqCe9HK3ZYAZ2ty3p2HibMLDuPGYtKsmVIfUwghhKhXpNTcJiqVCldbC1xtLejd0qVme0WVgeQcHRcyOhGbXkTPjGLiMopJLSgl2q4FBW2HcdXclUPfnwDAjgq+X/sNaoxs9OxIYLAfrdztsDpzjNKz57Du0QPLdm2V+phCCCGEYqTUKMxUo6almy0t3WwZEeJZs72wtKL6rE5GEWbpxVRkFBGXUUylrpSv2g7FU5vDkr2ZsLd6vawXzv5C/4T97DmZTOWUR2jpZkuLZubo/rMU85YtsY8ajspE/nULIYRovORbrp6ytzSla4AjXQMca7YZDEau5pdyIaMXcRnFDM0oIja9mKRcHUft/FB76dhW5sTxn04D4FeUzsc7v6TMzJJNpi1p6WZHCzcbPPdvR5OVie0dg7Bo1UqpjyiEEELUKSk1DYharcLXyQpfJysi2/539fLS8iriM3sSl1FMSEYRFhnFXMzSUq4zZV1gL1RGI5/tTa7Zf+7elYTmJPJVkh79QAMtXW0IUmlxXLUcuw7tcXzgfgU+nRBCCPHvSKlpBCzNNIT4OBDi43DN9oKSci5mRXExU8vkrGIuZmq5mFVMjHcnMq2c+LXckZS9SQB0TzvLrMO/cmH/Cb4tD6KFqw1Bbra03bASBypwu/9eLIKCFPh0QgghxI2RUtOIOViZ0cXfkS7+jtdsLyztS0JWMe6ZWi5maYnPLEaLNytaD6bY1JK9CTnsTahe7mHF1o2YlRYwtcgTY7s8Wrra0LEgmebbVuPQpzc+D09T4qMJIYQQ15FS0wTZW5oS5udImN//lp1wisvu4mKWlr6/n9GJzyhmbd4oHNIvc8LEiZLEXPYn5qK9+BvB5w6zKb2Mz7L9aOFqQ0tXW4Z9+y6Wdja4v/A8bsGBin0+IYQQTZOUGlHD1sKUTr7N6OT7P3PfTA1Hq69kZJaWi5nV8+pku1Sx0t6WeKM1OdpycrR5HL+Yyfgzx9BgZPBHfVA5JdLC1YYhlw/R/lg0qjuG4PPQFFxszFGpVBiNRlkiQgghRJ2SUiP+lo25CaE+DoT+cc3O0NbACErKK0nI0nIxU0tCWj6/WD1P1eUkCixsMerKOZyUR9eTp7C+nMAPv51jRXY0DlamBDtZ8NIXMyh396L0jQU0D/TA08ESVVkpKnNzVBqNkh9XCCFEAyWlRvxjVmYmdPB2oIO3A4R5Q1R7AF4qryIxu3oI62qIHevOhHPRaIMKKCipIDUzFYuSYgxXkpi0+gKoYjE3UTM9dgM9Yvdwcfh9mNx9H81drAloZok6JRkzf3/UZmaKfl4hhBD1m5QaUecszTS087KnnZc9dPSGu3sD1WtiJWZruZhawJ4unuReSSfY1o6kHB36SgOW2WmYVFawJVnLtt9nUPbQ5fDl9nmUm5qz7q1vaO5uS3MXG/xyU7C30GDRvDlqKyslP64QQoh6QkqNuG0sTDW09bSnrac9dPEDYCZQZTByNb+EhPSOHIlNxK1EQ2cdJGRrccwtpsTEnAxLR748cLnmWHMOfEF45gXW9buX3IFRNHe1oaVFFb5nDuIe0gbbLp0V+pRCCCGUIqVGKE6jVuHnZI2fkzW087zmtTxdPxIy70Z3OZNpOkjM1pKYrUNnZkmeuS37Kmw5c+wqAJ2y4nh7/2ccsHVl/t1v0NzFhuau1oSd2oWrhRrvqCE4BPgq8RGFEELcBlJqRL3maG1G10AnugY6XbO97Nk+JOXoeDirmMTsEhKztRhOZnDcvRVpls24mFU9Bw/n4NMd36PRZjP5ZAnpLTrQ3NWazvpMOh/bjmWnTnjffw9uduZyN5YQQjRwUmrqwvonybXsRrNeI1Bb2iqdpkmwMNXQ2sOO1h52/914d0cMhimkFZbSNVtHYpaWhGwtCRk9yEhNJtnOnfyiMjKKynBIOszgU9EcSkxj7FU3rM00NHe14eGdn2FjYUrFpIfwC2mNn5M1ZiZq5T6oEEKIGyal5t/KOIP+6E/8lD0Ei7VbGNv/GNY97wGPDkona5LUahXezazwbmZF3yCX6o2j5wEwsqSCxBwtiVlask+q2GuvIl5jj0atQldexdmUfLxjj2NmqGSyUx8ydmajUasYn3OKqNNbyOjSF/3
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def do_int(n):\n",
    "    # integration limits\n",
    "    x1 = 0.0\n",
    "    x2 = np.pi/2\n",
    "    x = np.linspace(x1, x2, n)\n",
    "    dx = (x2-x1)/(n-1)\n",
    "    \n",
    "    # Choose the function to integrate\n",
    "    if True:\n",
    "        f = np.cos(x)\n",
    "        I0 = 1.0\n",
    "    else: \n",
    "        # Polynomial\n",
    "        a = 4\n",
    "        I0 = (np.pi/2)**(a+1)/(a+1)  # The analytic value of the integral\n",
    "        f = x**a  # The function to integrate\n",
    "    \n",
    "    I1 = np.sum(f[:-1])*dx    # rectangle rule\n",
    "    I2 = np.sum(f[1:-1]*dx) + 0.5*dx*(f[0]+f[-1])    # trapzoidal rule\n",
    "    I3 = np.sum(4*f[1:-1:2]*dx) + np.sum(2*f[2:-1:2]*dx) + f[0]*dx + f[-1]*dx  #simpson's rule\n",
    "    I3 = I3/3.0\n",
    "    return abs((I1-I0)/I0), abs((I2-I0)/I0), abs((I3-I0)/I0)\n",
    "\n",
    "# try different numbers of points, generated as 2^n + 1 so the number is odd\n",
    "n_vals = [2**(i+1) + 1 for i in range(11)]\n",
    "err1=np.array([])\n",
    "err2=np.array([])\n",
    "err3=np.array([])\n",
    "for n in n_vals:\n",
    "    e1, e2, e3 = do_int(n)\n",
    "    print(n, e1, e2, e3)\n",
    "    err1 = np.append(err1, e1)\n",
    "    err2 = np.append(err2, e2)\n",
    "    err3 = np.append(err3, e3)\n",
    "    \n",
    "plt.plot(n_vals, err1, label='Rectangular')\n",
    "plt.plot(n_vals, err2, label='Trapezoidal')\n",
    "plt.plot(n_vals, err3, label='Simpson')\n",
    "\n",
    "n = np.array(n_vals)\n",
    "dx = (np.pi/2)/(n-1)\n",
    "# Plot the Euler Maclaurin error formulas\n",
    "plt.plot(n_vals, dx/2, \":\")\n",
    "plt.plot(n_vals, dx**2/12, \":\")\n",
    "plt.plot(n_vals, dx**4/180, \":\")\n",
    "\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0bb8f63",
   "metadata": {},
   "source": [
    "Notes:\n",
    "- The scalings with the number of points are $1/N$ for rectangular, $1/N^2$ for trapezoidal, and $1/N^4$ for Simpson's rule. (You might have been expecting $1/N^3$ for Simpson's rule but because of the way it was constructed using double intervals, the third order term is antisymmetric over the double interval and cancels.)\n",
    "- If you try integrating polynomials, you should find that the error goes to zero (machine precision) for cubic polynomials and below (Simpson), linear or below (trapezoid) or constant (rectangular).\n",
    "- Some special cases can give surprising results. For example, for $\\int \\sin^2(x) dx$, all the methods can give exact results. To see why, you can rewrite $\\sin^2 x=(1-\\cos(2x))/2$. If you have an odd number of sample points, the cos term averages to zero and the remaining term is a constant, which all three methods will fit exactly. \n",
    "- It is possible to derive analytic expressions for the errors in the different approximations (these are known as [Euler Maclaurin](https://en.wikipedia.org/wiki/Euler–Maclaurin_formula) formulae). For the trapezoidal rule, the leading term is \n",
    "$$\\epsilon \\approx {1\\over 12} (\\Delta x)^2 \\left[f^\\prime(a) - f^\\prime(b)\\right]$$\n",
    "and for Simpson's rule,\n",
    "$$\\epsilon \\approx {1\\over 180} (\\Delta x)^4 \\left[f^{\\prime\\prime\\prime}(a) - f^{\\prime\\prime\\prime}(b)\\right].$$\n",
    "These are plotted as dotted lines in the Figure, you can see that they agree remarkably well with the numerical results. (Note that when plotting we use the fact that the derivative terms are equal to unity for $\\cos(x)$). "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5381b56",
   "metadata": {},
   "source": [
    "## Gaussian quadrature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7eb969c8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "quad gives I =   1.493648265624854 1.6582826951881447e-14\n",
      "1 1 2.0 0.33900332898208696\n",
      "2 3 1.4330626211475785 -0.0405621898217985\n",
      "3 5 1.4986795956600294 0.0033684838331537233\n",
      "4 7 1.4933346224495394 -0.00020998462792936129\n",
      "5 9 1.493663920702629 1.048110062810725e-05\n",
      "6 11 1.4936476141506054 -4.361630937756038e-07\n",
      "7 13 1.493648288869414 1.5562271646873954e-08\n",
      "8 15 1.4936482648990144 -4.859508274756656e-10\n",
      "9 17 1.4936482656450039 1.3490379379714472e-11\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGjCAYAAADUwuRbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABFWUlEQVR4nO3dZ3gU5cLG8f/sppEKIRAIJPReEkoIIYCiUUSkiAoqagAFRQQ0loN6jujxVTxybEioFlBUEJUgKKKiSEkwEAjSewmEJNQ0SN19P3jMeXlBaZtMsrl/17Ufdnayc496ufc188zzGHa73Y6IiIhIJWcxO4CIiIiII6jUiIiIiFNQqRERERGnoFIjIiIiTkGlRkRERJyCSo2IiIg4BZUaERERcQoqNSIiIuIUVGpERETEKajUiIiIiFNQqRERERGnUKlKzdKlS2nRogXNmjXjvffeMzuOiIiIVCBGZVnQsri4mNatW/Pzzz/j5+dHp06dSEhIoGbNmmZHExERkQqg0lypSUpKok2bNtSrVw9vb2/69OnD999/b3YsERERqSBcyutAq1atYvLkySQnJ3Ps2DEWLVrEwIEDz9snLi6OyZMnk56eTmhoKO+++y5dunQBIC0tjXr16pXuW69ePY4ePXrZx7fZbKSlpeHj44NhGA45JxERESlbdrudnJwcgoKCsFj++lpMuZWavLw8QkNDGTFiBIMGDbrg8wULFhAbG8uMGTOIiIjg7bffpnfv3uzatYvatWtf8fEKCgooKCgofX/06FFat259TecgIiIi5khNTaV+/fp/uU+5lZo+ffrQp0+fP/38zTffZOTIkQwfPhyAGTNm8M033/DBBx8wYcIEgoKCzrsyc/To0dKrOBczadIkXnrppQu2p6am4uvrew1nIiIiIuUlOzub4OBgfHx8LrmvKQOFDcM47/ZTYWEhnp6efPHFF+fdkoqJieHMmTMsXryY4uJiWrVqxcqVKy9roPD/v1Lzxz+UrKwslRoREZFKIjs7Gz8/v8v6/S63KzV/5cSJE5SUlBAYGHje9sDAQHbu3AmAi4sLb7zxBr169cJms/HMM8/85ZNP7u7uuLu7l2luERERqTgqRKm5XP3796d///5mxxAREZEKqEI80h0QEIDVaiUjI+O87RkZGdSpU8ekVCIiIlKZVIhS4+bmRqdOnVixYkXpNpvNxooVK4iMjDQxmYiIiFQW5Xb7KTc3l71795a+P3DgACkpKfj7+xMSEkJsbCwxMTF07tyZLl268Pbbb5OXl1f6NJSIiIjIXym3UrNhwwZ69epV+j42Nhb4/QmnOXPmMGTIEI4fP84LL7xAeno6YWFhfPfddxcMHhYRERG5mEqz9tO1upJHwkRERKRiuJLf7woxpkZERETkWqnUiIiIiFNQqRERERGnoFIjIiIiTkGlRkRERJyCSs01stvtjPlkIws3pFJiqxIPkomIiFRIKjXX6McdmXyz5RhPf/Ebvd9exXdb06kiT8mLiIhUKCo116hHswCeu7Ul1T1d2ZuZyyPzkhk4LYG1e0+YHU1ERKRK0eR7jvr+/CJmr9rP+2sOcLawBICopjV5pndLQoOrO/x4IiIiVcGV/H6r1DjY8ZwC4n7eyye/HqKo5Pd/tLe0qcNTvZvTtLZPmR1XRETEGanUXER5L5OQeuosb/+4h0WbjmCzg8WAOzrW5/GbmlOverUyP76IiIgzUKm5CLPWftqdkcO/l+/i++0ZALhZLQztGsKYXk0J8HYvtxwiIiKVkUrNRZi9oOWmw6d5/btdJO4/CYCXm5UHezRmZI9G+Hi4lnseERGRykCl5iLMLjXw+5w2a/ae4PXvdrHlaBYANTxdGdOrKfd1bYCHq9WUXCIiIhWVSs1FVIRS8we73c53W9P59/e72Hc8D4C6fh6Mv7EZd3aqj4tVT9qLiIiASs1FVaRS84fiEhtfbTzK2z/uJi0rH4DGAV48eXML+rStg8VimJxQRETEXCo1F1ERS80f8otK+OTXw8T9vJdTeYUAtK3ny9O9W9KzWQCGoXIjIiJVk0rNRVTkUvOHnPwi3l9zgPdWHyC3oBiAro39eeaWlnQMqWFyOhERkfKnUnMRlaHU/OFkbgHTVu7j43WHKCy2ARDdKpCne7egRR1N4CciIlWHSs1FVKZS84e0M+d458c9LExOxWYHw4Dbw+rxxE3NCfb3NDueiIhImVOpuYjKWGr+sDczlzd/2MW3W9IBcLUa3NslhMduaEYtH03gJyIizkul5iIqc6n5w5YjWby+fCer9/y+Ang1VysjujdkVM8m+FXTBH4iIuJ8VGouwhlKzR8S9v0+gV9K6hkA/Kq5Mvr6JsRENqSamybwExER56FScxHOVGrg9wn8vt+ewb+X72JPZi4AtX3cGR/djMGdg3HVBH4iIuIEVGouwtlKzR9KbHbiNx3lrR93c+T0OQAa1vTkiZua0699kCbwExGRSk2l5iKctdT8oaC4hM9+PczUn/dyIvf3Cfxa1fXlmd4tuL5FLU3gJyIilZJKzUU4e6n5Q15BMR+uPcDMX/aT858J/MIb1uCZW1oS3tDf5HQiIiJXRqXmIqpKqfnD6bxCZvyyjzkJByn4zwR+vVrU4uneLWkd5PznLyIizkGl5iKqWqn5Q3pWPlN+2sOC9amU2H7/V90/NIjYm5rTMMDL5HQiIiJ/TaXmIqpqqfnDgRN5vPnDbpZsTgPAxWIwJDyYcTc2I9DXw+R0IiIiF6dScxFVvdT8YevRLP79/S5W7joOgIerhZhuDRl9XROqe7qZnE5EROR8KjUXoVJzvqQDp3j9u51sOHQaAB8PFx65rgnDoxri6eZicjoREZHfqdRchErNhex2Oz/vyuT173axMz0HgABvd8be0JR7uoTg5qIJ/ERExFwqNRehUvPnbDY7S35L443vd3P41FkAgv2r8UR0cwaE1cOqCfxERMQkKjUXoVJzaYXFNhZsSGXKij0czykAoEWgD0/e3JybWgdqAj8RESl3KjUXoVJz+c4VljAn4SDTV+4lO//3Cfw6hlTnnwPa0raen8npRESkKlGpuQiVmiuXdbaImav28cHaA+QX2bAYENOtIU/e3AJvdw0mFhGRsqdScxEqNVcvIzuf//lmR+kcN3V8PZjYrzW3tK2jW1IiIlKmruT3W4+3yCUF+nrw7j0d+GhEFxrU9CQ9O5/Rn2zkwbkbSP3PwGIRERGzqdTIZevZvBbLH+/J2Bua4mo1+GlnJje99QvTV+6jqMRmdjwREaniVGrkini4Wnny5hYsG9+TiEb+5BfZ+Nd3O+k7ZTUbDp4yO56IiFRhKjVyVZrW9mb+qK68cVco/l5u7M7I5c4ZiUz48jfOnC00O56IiFRBKjVy1QzD4I5O9VkRex13hwcDMH99Kje88QtfJh+hioxBFxGRCkKlRq5ZDS83XrujPQsfiaR5oDen8gp5cuFm7pm9jr2ZuWbHExGRKkKlRhwmvKE/S8f24G+3tMTD1cK6/afo884q3vx+F/lFJWbHExERJ6dSIw7l5mJh9PVN+OGJ6+jVohZFJXam/LSX3m+vYtXu42bHExERJ6ZSI2Ui2N+TD4aFM31oRwJ93Tl08iwPfJDE2M82kZmTb3Y8ERFxQio1UmYMw6BPu7r8GHsdw6MaYjFgyeY0bnzjFz5OPEiJTQOJRUTEcbRMgpSbLUeyeD5+C78dyQIgNLg6rwzUIpkiIvLntEyCVEjt6vux6NEoXurfBm93FzannqH/1DW8vHQ7uQXFZscTEZFKrtKUmtTUVK6//npat25N+/btWbhwodmR5CpYLQYx3Rqy4snr6Nu+LjY7vL/mADe9+QvfbU3X3DYiInLVKs3tp2PHjpGRkUFYWBjp6el06tSJ3bt34+XldVl/r9tPFdPKXZn8Y/FWUk+dAyC6VW1e7N+G+jU8TU4mIiIVgVPefqpbty5hYWEA1KlTh4CAAE6d0lpDld31LWrz/ePXMaZXE1ytBj/uyOSmN1cx8xctkikiIlfGYaVm1apV9OvXj6CgIAzDID4+/oJ94uLiaNiwIR4eHkRERJCUlHRVx0pOTqa
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Use Gaussian quadrature to evalute the integral of polynomials from x=-1 to +1\n",
    "\n",
    "import numpy as np\n",
    "import scipy.integrate\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Some choices for functions to integrate:\n",
    "# The first one is a random polynomial (the size parameter sets the degree)\n",
    "#func = np.polynomial.Polynomial(np.random.randint(-10,high=10,size=20)) \n",
    "#func = np.cos\n",
    "#func = np.exp\n",
    "func = lambda x: np.exp(-x**2)\n",
    "\n",
    "# Integrate this with quad to get the true value\n",
    "I0, err0 = scipy.integrate.quad(func, -1.0, 1.0)\n",
    "print('quad gives I =  ', I0, err0)\n",
    "\n",
    "# then use Gaussian quadrature with different numbers of points\n",
    "n_vec = np.array([])\n",
    "err_vec = np.array([])\n",
    "\n",
    "for npoints in range(1,10):\n",
    "    I = 0.0\n",
    "    x, w = np.polynomial.legendre.leggauss(npoints)\n",
    "    for i in range(npoints):\n",
    "        I += w[i] * func(x[i])\n",
    "    err = (I-I0)/I0\n",
    "    print(npoints, 2*npoints-1, I, err)\n",
    "    n_vec = np.append(n_vec, npoints)\n",
    "    err_vec = np.append(err_vec, abs(err))\n",
    "    \n",
    "plt.plot(n_vec, err_vec)\n",
    "plt.yscale('log')\n",
    "plt.xscale('linear')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22c1aadf",
   "metadata": {},
   "source": [
    "- You can use the code above to check that the answer is exact for polynomials of degree $2n-1$ and smaller\n",
    "- The scaling of the error with number of points depends on the function you are integrating, but decreases approximately exponentially (rather than power law with the Newton-Cotes methods)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1650be55",
   "metadata": {},
   "source": [
    "Let's try the Gauss-Hermite quadrature which has weighting function $W(x)=e^{-x^2}$ for integrals from $-\\infty$ to $+\\infty$. In the next example I try $x^4 e^{-x^2}$ which becomes exact once we go to 3 terms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "38475475",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "quad gives I =   1.3293403881791366 1.5859180523983767e-08\n",
      "1 1 0.0 -1.0\n",
      "2 3 0.4431134627263788 -0.6666666666666666\n",
      "3 5 1.3293403881791368 1.6703367090890604e-16\n",
      "4 7 1.3293403881791377 8.351683545445302e-16\n",
      "5 9 1.3293403881791368 1.6703367090890604e-16\n",
      "6 11 1.3293403881791355 -8.351683545445302e-16\n",
      "7 13 1.3293403881791372 5.011010127267181e-16\n",
      "8 15 1.3293403881791364 -1.6703367090890604e-16\n",
      "9 17 1.3293403881791372 5.011010127267181e-16\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGdCAYAAADqsoKGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8K0lEQVR4nO3dfXDU5333+8/uSloJaSWQhFYISciAJGywJYcHWdhuoFbCkV1uPxw7nI5zqkJDp2eM41TTTCG5C83JJKR3OwynY2qadGziJnZo7/s2SesEu1bsEAcaDETEqW2BQMYY0OoBpNWu0Era/Z0/pF0QEqCH3f3tw/s1s+PZB+3vu4DZD9d1fa/LYhiGIQAAgDhnNbsAAACAcCDUAACAhECoAQAACYFQAwAAEgKhBgAAJARCDQAASAiEGgAAkBAINQAAICGkmF1AtAQCAV28eFEOh0MWi8XscgAAwCQYhqG+vj4VFRXJar31WEzShJqLFy+qpKTE7DIAAMA0nD9/XsXFxbd8TdKEGofDIWnkFyU7O9vkagAAwGS43W6VlJSEvsdvJWlCTXDKKTs7m1ADAECcmczSERYKAwCAhECoAQAACYFQAwAAEgKhBgAAJARCDQAASAiEGgAAkBAINQAAICEQagAAQEIg1AAAgIRAqAEAAAmBUAMAABICoQYAACSEpDnQMlI+vdKvF945o5yM1LG3WWPvZ9lTJnUYFwAAmB5CzQx9euWqfvjrT277OpvVouz0lFDIyb4xBN1wy74uHDkIRAAA3BahZobm5aTryw+Vy311SL03uQ0OB+QPGLrSP6Qr/UNTvobVonEh6HahKPgahz1FViuBCACQ+Ag1M7QgL1ONn6u45WsGhvxjg07/+OBzs1DkGw4oYEg9/UPqmWYgcqTfYiToFjdHOoEIABA/CDVRkJ5qU3qqTc7s9Cn/7MCQ/5ajQLcKRQNDI4EoeH+qLBbJYU/R7FlpunOeQ3//h/fKnmKb8vsAABANhJoYFwxEBdMNRAMjgadngtGh2wUiw5DcA8NyDwzrk8v9Otp2WQ+Wz43ApwQAYOYINQksFIgcUw9EvmF/KPB8498+0C9Pd+mUy0OoAQDELPapwYTsKSNhaHGBQ/eWzJYknXb1mVsUAAC3QKjBbVUUOiRJpwg1AIAYRqjBbVU4R0LNaZdHhmGYXA0AABMj1OC2yvIylWK1qM83rEu9A2aXAwDAhAg1uK20FKsWzs2UJLUwBQUAiFGEGkxKeWgKilADAIhNhBpMSuVoqGlp95hcCQAAEyPUYFIqnFmSpNMdjNQAAGIToQaTcn0HVCBABxQAIPYQajApC/IylZZi1dUhvz69ctXscgAAGCeuQs3jjz+uOXPm6MknnzS7lKRjs1q0aO7IFBSb8AEAYlFchZrnnntOL7/8stllJK3K0XU1tHUDAGJRXIWaNWvWyOFwmF1G0qKtGwAQy8IWag4dOqT169erqKhIFotFBw4cGPeaPXv2qKysTOnp6aqpqdHRo0fDdXlEQait20VbNwAg9oQt1Hi9XlVVVWnPnj0TPr9//341NjZqx44dOnHihKqqqrRu3Tp1dHSEXlNdXa1ly5aNu128eDFcZWIGgh1QZzo98tMBBQCIMSnheqP6+nrV19ff9Pldu3Zp8+bN2rhxoyRp7969ev311/Xiiy9q69atkqTm5uZwlSOfzyefzxe673a7w/beyap4ToYyUm26OuTXuW6vFo4uHAYAIBZEZU3N4OCgjh8/rrq6umsXtlpVV1enI0eOROSaO3fuVE5OTuhWUlISkeskE6vVonInHVAAgNgUlVDT1dUlv98vp9M55nGn06n29vZJv09dXZ2eeuop/fSnP1VxcfEtA9G2bdvU29sbup0/f37a9eOa8oKRKahTrKsBAMSYsE0/RcNbb7016dfa7XbZ7fYIVpOcKgtp6wYAxKaojNTk5+fLZrPJ5XKNedzlcqmwsDAaJSBMaOsGAMSqqISatLQ0LV++XE1NTaHHAoGAmpqaVFtbG40SECbBtu6znV4NDgdMrgYAgGvCNv3k8XjU2toaut/W1qbm5mbl5uaqtLRUjY2Namho0IoVK7Rq1Srt3r1bXq831A2F+DAvJ10Oe4r6fMP6uNsbavMGAMBsYQs1x44d09q1a0P3GxsbJUkNDQ3at2+fNmzYoM7OTm3fvl3t7e2qrq7WwYMHxy0eRmyzWEY6oE580qNTrj5CDQAgZoQt1KxZs0aGcesN2bZs2aItW7aE65IwSYXTMRJq2vuke8yuBgCAEXF19hNiQ3CxMG3dAIBYQqjBlFWGQg0dUACA2EGowZRVjO4q/HG3VwNDfpOrAQBgBKEGUzbXYdfsWakKGCOHWwIAEAsINZgyi8WiioLgJnyEGgBAbCDUYFoqCjnYEgAQWwg1mJYKFgsDAGIMoQbTwmndAIBYQ6jBtAQ7oD653K/+wWGTqwEAgFCDacrLsis/K02S1NrBaA0AwHyEGkxbcF1NSzvragAA5iPUYNqCoeY0IzUAgBhAqMG00QEFAIglhBpMW3Cx8CmmnwAAMYBQg2kLntZ9sXdAfQNDJlcDAEh2hBpMW05Gqgqz0yWxXw0AwHyEGsxI+egU1GnW1QAATEaowYxUBtu6CTUAAJMRajAjobZupp8AACYj1GBGKgoZqQEAxAZCDWakvGBkTU1nn089/YMmVwMASGaEGsxIpj1F82dnSKIDCgBgLkINZqySKSgAQAwg1GDGaOsGAMQCQg1mrJLTugEAMYBQgxm7/mBLwzBMrgYAkKwINZixxQVZslikK/1D6vLQAQUAMAehBjOWnmrTgtxZklhXAwAwT9yEmp6eHq1YsULV1dVatmyZvve975ldEq5Tft0UFAAAZkgxu4DJcjgcOnTokGbNmiWv16tly5bpiSeeUF5entmlQSOLhf/jA5da2KsGAGCSuBmpsdlsmjVrZIrD5/PJMAwWpcYQ2roBAGYLW6g5dOiQ1q9fr6KiIlksFh04cGDca/bs2aOysjKlp6erpqZGR48endI1enp6VFVVpeLiYn31q19Vfn5+mKrHTF2/AR9hEwBghrCFGq/Xq6qqKu3Zs2fC5/fv36/Gxkbt2LFDJ06cUFVVldatW6eOjo7Qa4LrZW68Xbx4UZI0e/ZsnTx5Um1tbXrllVfkcrnCVT5m6I78TNmsFvUNDMvl9pldDgAgCYVtTU19fb3q6+tv+vyuXbu0efNmbdy4UZK0d+9evf7663rxxRe1detWSVJzc/OkruV0OlVVVaVf/vKXevLJJyd8jc/nk8937cvV7XZP8pNgOuwpNt2Rn6nWDo9aXH0qzEk3uyQAQJKJypqawcFBHT9+XHV1ddcubLWqrq5OR44cmdR7uFwu9fWNrNfo7e3VoUOHVFlZedPX79y5Uzk5OaFbSUnJzD4EbquCdTUAABNFJdR0dXXJ7/fL6XSOedzpdKq9vX1S73Hu3Dk9+OCDqqqq0oMPPqhnn31Wd999901fv23bNvX29oZu58+fn9FnwO2VF9DWDQAwT9y0dK9atWrS01OSZLfbZbfbI1cQxrm2WJi2bgBA9EVlpCY/P182m23cwl6Xy6XCwsJolIAoCE4/tbr6FAjQAQUAiK6ohJq0tDQtX75cTU1NoccCgYCamppUW1sbjRIQBQvyMpVms8o76NeFnqtmlwMASDJhm37yeDxqbW0N3W9ra1Nzc7Nyc3NVWlqqxsZGNTQ0aMWKFVq1apV2794tr9cb6oZC/Eu1WbVwbqY+au/T6Y4+lYyeBwUAQDSELdQcO3ZMa9euDd1vbGyUJDU0NGjfvn3asGGDOjs7tX37drW3t6u6uloHDx4ct3gY8a3C6dBH7X1qaffo95fwewsAiJ6whZo1a9bcdifZLVu2aMuWLeG6JGIQbd0AALPEzdlPiA+h07o7CDUAgOgi1CCsKkdDzWmXR346oAAAUUSoQViV5M6SPcUq33BA5y/3m10OACCJEGoQVjarReWj62paWFcDAIgiQg3CrqIgOAVFqAEARA+hBmFXwXEJAAATEGoQdrR1AwDMQKhB2AVP6z7b6dWQP2ByNQCAZEGoQdj
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Gauss-Hermite quadrature\n",
    "\n",
    "import numpy as np\n",
    "import scipy.integrate\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Function to integrate\n",
    "func = lambda x: x**4 * np.exp(-x**2)\n",
    "\n",
    "# Integrate this with quad to get true value\n",
    "I0, err0 = scipy.integrate.quad(func, -np.inf , np.inf)\n",
    "print('quad gives I =  ', I0, err0)\n",
    "\n",
    "# Now try the Gauss-Hermite quadrature with different numbers of points and plot the error\n",
    "n_vec = np.array([])\n",
    "err_vec = np.array([])\n",
    "\n",
    "for npoints in range(1,10):\n",
    "    I = 0.0\n",
    "    x, w = np.polynomial.hermite.hermgauss(npoints)\n",
    "    for i in range(npoints):\n",
    "        # note that on the next line we include only the part of the function multiplying \n",
    "        # the exp(-x^2), because the exp(-x^2) is already taken into account in the weights\n",
    "        I += w[i] * func(x[i]) /  np.exp(-x[i]**2)  \n",
    "    err = (I-I0)/I0\n",
    "    print(npoints, 2*npoints-1, I, err)\n",
    "    n_vec = np.append(n_vec, npoints)\n",
    "    err_vec = np.append(err_vec, abs(err))\n",
    "    \n",
    "plt.plot(n_vec, err_vec)\n",
    "plt.yscale('log')\n",
    "plt.xscale('linear')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7749401",
   "metadata": {},
   "source": [
    "## Maxwell-Boltzmann distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d023c71",
   "metadata": {},
   "source": [
    "First summarize some results for the Maxwell-Boltzmann distribution from the Wikipedia page: \n",
    "\n",
    "In 1D, the distribution of velocities is\n",
    "\n",
    "$$f(v)  = \\left({m\\over 2\\pi k_BT}\\right)^{3/2} e^{-mv^2/2k_BT} ,$$\n",
    "\n",
    "which is normalized so that\n",
    "\n",
    "$$\\int_0^\\infty d^3v  f(v) =1.$$\n",
    "\n",
    "The mean speed is\n",
    "\n",
    "$$\\langle v\\rangle = \\int_0^\\infty d^3v\\ v f(v).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "563e6fa9",
   "metadata": {},
   "source": [
    "To get this into a simpler form, change integration variable to $x=(m/2k_BT)^{1/2} v$ and use the spherical volume element $d^3 v= 4\\pi v^2 dv$. This gives\n",
    "\n",
    "$$\\langle v\\rangle = \\left({k_BT\\over \\pi m}\\right)^{1/2}\\ 4\\sqrt{2}\\ I\\hspace{1cm} I = \\int_0^\\infty dx\\ x^3\\ e^{-x^2}.$$\n",
    "\n",
    "We can evaluate the integral $I$ numerically (aiming for an error of 0.1%). (Spoiler: the analytic answer is 1/2).\n",
    "\n",
    "For Gaussian quadrature, you might think of using the Gauss-Hermite polynomials since we have an $e^{-x^2}$ in the integral, but the limits are from zero to infinity, not -infinity to infinity. One way to do it is to write the integral as\n",
    "\n",
    "$${1\\over 2}\\int_{-\\infty}^\\infty dx\\ \\left|x\\right|^3\\ e^{-x^2}.$$\n",
    "\n",
    "Note that this is not the same as having a cubic polynomial, so the usual rule about how many terms we need does not apply.\n",
    "\n",
    "The other way to do it is to rewrite\n",
    "\n",
    "$$\\int_0^\\infty dx\\ x^3\\ e^{-x^2}= \\int_0^\\infty dy\\ {y\\over 2} e^{-y}$$\n",
    "\n",
    "and use Laguerre polynomials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "id": "94aa2a9a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "quad gives 0.5 with error 3.36729e-05; number of evaluations=75\n",
      "Simpson gives 0.499296 with error 0.00140896; number of points=21\n",
      "Gauss-Hermite gives 0.500992 with error 0.0019831; number of terms=15\n",
      "Gauss-Laguerre gives 0.5 with error 0; number of terms=1\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.integrate\n",
    "\n",
    "func = lambda x: x**3 * np.exp(-x**2)\n",
    "I0 = 0.5  # the analytic result\n",
    "\n",
    "# First use quad\n",
    "# We'll set the relative error we are looking for to 1e-3 and also ask for \n",
    "# full output which gives us information such as how many function evaluations were done\n",
    "I1, err1, info = scipy.integrate.quad(func,0.0,np.inf, full_output=True, epsrel=1e-3)\n",
    "print(\"quad gives %lg with error %lg; number of evaluations=%d\" % (I1, err1, info['neval']))\n",
    "\n",
    "# Now Simpson's rule\n",
    "# Sample the function first\n",
    "npoints = 21\n",
    "xp = np.linspace(0.0,7.0,npoints)  # integrate to x=7\n",
    "fp = func(xp)\n",
    "I2 = scipy.integrate.simpson(fp, xp)\n",
    "print(\"Simpson gives %lg with error %lg; number of points=%d\" %(I2, abs(I2-I0)/I0, npoints))\n",
    "\n",
    "# Now use Gaussian quadrature\n",
    "npoints = 15\n",
    "x, w = np.polynomial.hermite.hermgauss(npoints)\n",
    "I3 = 0.5 * np.sum(w * abs(x)**3)\n",
    "print(\"Gauss-Hermite gives %lg with error %lg; number of terms=%d\" %(I3, abs(I3-I0)/I0, npoints))\n",
    "\n",
    "# Gauss-Laguerre\n",
    "npoints = 1\n",
    "x, w = np.polynomial.laguerre.laggauss(npoints)\n",
    "I4 = 0.5 * np.sum(w * x)\n",
    "print(\"Gauss-Laguerre gives %lg with error %lg; number of terms=%d\" %(I4, abs(I4-I0)/I0, npoints))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f9b0ea2",
   "metadata": {},
   "source": [
    "In the Simpson's method, we integrate to $x=7$. You can experiment with this and see how much it changes the answer. At $x=7$, the integrand becomes comparable to the machine precision, so it seems a reasonable place to stop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da407d88",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}