L05.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Lecture 5 - Integration: Romberg, Adaptive quadrature, Gauss quadrature\n",
    "\n",
    "1. Implement these three methods to estimate: \n",
    "$$ \\int_0^{0.8} 0.2+25𝑥−200𝑥^2+675𝑥^3−900𝑥^4+400𝑥^5 dx $$\n",
    "\n",
    "The exact value is 1.640533. Use our `trapezoid` function from last time.\n",
    "\n",
    "Test against built-in `scipy.integrate` functions.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def f(x):\n",
    "    return 0.2 + 25*x - 200*x**2 + 675*x**3 - 900*x**4 + 400*x**5\n",
    "\n",
    "def trapezoid(f,a,b,n):\n",
    "    x = np.linspace(a,b,n+1) # n segments, n+1 points between a and b\n",
    "    I = (b-a)/(2*n) * (f(x[0]) + 2*np.sum(f(x[1:-1])) + f(x[-1]))\n",
    "    return I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def romberg(f, a, b, tol=1e-4, maxiter=10):\n",
    "    # estimate integral of f(x) from a to b\n",
    "    # stop when tolerance reached or maxiter exceeded\n",
    "    # thanks to Simeon for improving this! (Winter 2025)\n",
    "    \n",
    "    R = np.zeros((maxiter, maxiter))  # Romberg matrix\n",
    "    \n",
    "    # Add a new row j every iteration, and update all columns in that row\n",
    "    for j in range(maxiter):\n",
    "        n = 2**j # n = 1, 2, 4, 8 etc. segments\n",
    "        R[j,0] = trapezoid(f,a,b,n) # first column always from trapezoid rule  \n",
    "\n",
    "        # update all other columns in row j\n",
    "        for k in range(1,j+1):\n",
    "            R[j,k] = (4**k * R[j,k-1] - R[j-1,k-1]) / (4**k - 1)\n",
    "        \n",
    "        # how much has the bottom-right value changed from the last column?\n",
    "        if j > 0 and abs(R[j,k] - R[j,k-1]) < tol:\n",
    "            break\n",
    "    \n",
    "    I = R[j,k] # return bottom-right value\n",
    "    num_evals = n+1\n",
    "    return I, num_evals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Romberg: I=1.6405333, num_evals = 9\n"
     ]
    }
   ],
   "source": [
    "a = 0\n",
    "b = 0.8\n",
    "I, num_evals = romberg(f, a, b, tol=1e-4)\n",
    "print('Romberg: I=%0.7f, num_evals = %d' % (I, num_evals))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare to built-in function. If we want to see the full matrix, run with `show=True`. The full documentation is [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.romberg.html#scipy.integrate.romberg)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.6405333333333363\n",
      "Romberg integration of <function vectorize1.<locals>.vfunc at 0x1417e7e20> from [0, 0.8]\n",
      "\n",
      " Steps  StepSize   Results\n",
      "     1  0.800000  0.172800 \n",
      "     2  0.400000  1.068800  1.367467 \n",
      "     4  0.200000  1.484800  1.623467  1.640533 \n",
      "     8  0.100000  1.600800  1.639467  1.640533  1.640533 \n",
      "\n",
      "The final result is 1.6405333333333363 after 9 function evaluations.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/3q/q4pw14qn3r95vp8fslqnfl340000gn/T/ipykernel_34838/391658350.py:2: DeprecationWarning: `scipy.integrate.romberg` is deprecated as of SciPy 1.12.0and will be removed in SciPy 1.15.0. Please use`scipy.integrate.quad` instead.\n",
      "  result = integrate.romberg(f, a, b, tol=1e-4)\n",
      "/var/folders/3q/q4pw14qn3r95vp8fslqnfl340000gn/T/ipykernel_34838/391658350.py:4: DeprecationWarning: `scipy.integrate.romberg` is deprecated as of SciPy 1.12.0and will be removed in SciPy 1.15.0. Please use`scipy.integrate.quad` instead.\n",
      "  integrate.romberg(f, a, b, tol=1e-4, show=True) # print full matrix\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "np.float64(1.6405333333333363)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy import integrate\n",
    "result = integrate.romberg(f, a, b, tol=1e-4)\n",
    "print(result) # print answer only\n",
    "integrate.romberg(f, a, b, tol=1e-4, show=True) # print full matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Adaptive quadrature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def quadapt(f, a, b, tol=1e-4):\n",
    "    # returns integral estimate and number of evaluations\n",
    "    I1 = trapezoid(f, a, b, 1)\n",
    "    I2 = trapezoid(f, a, b, 2)\n",
    "    \n",
    "    if abs(I2-I1) < tol:\n",
    "        return (4*I2 - I1)/3, 2+1 # romberg correction\n",
    "    else:\n",
    "        midpoint = (a+b)/2\n",
    "        Ia,na = quadapt(f, a, midpoint, tol)\n",
    "        Ib,nb = quadapt(f, midpoint, b, tol)\n",
    "        return Ia + Ib, na+nb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adaptive: I=1.6405331, num_evals = 129\n"
     ]
    }
   ],
   "source": [
    "I,num_evals = quadapt(f, a, b, tol=1e-4)\n",
    "print('Adaptive: I=%0.7f, num_evals = %d' % (I, num_evals))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare to built in `quad`. Documentation [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html#scipy.integrate.quad)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "result, err, infodict = integrate.quad(f, a, b, epsabs=1e-4, full_output=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.6405333333333307\n"
     ]
    }
   ],
   "source": [
    "print(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.8213578793317207e-14\n"
     ]
    }
   ],
   "source": [
    "print(err)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "21\n"
     ]
    }
   ],
   "source": [
    "print(infodict['neval'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The built-in function is much more efficient and accurate than our simplified version."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Gauss Quadrature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss2(f,a,b):\n",
    "    c = np.array([1,1])\n",
    "    x = np.array([-1/np.sqrt(3), 1/np.sqrt(3)])\n",
    "    \n",
    "    y = (a+b)/2 + (b-a)/2 * x  # map from [-1,1] to [a,b]\n",
    "    I = (b-a)/2 * sum(c * f(y))  # estimate integral\n",
    "    return I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.8225777777777772\n"
     ]
    }
   ],
   "source": [
    "I = gauss2(f,a,b)\n",
    "print(I)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can write similar functions for 3, 4, 5, 6-point Gauss quadrature using the values from Table 22.1 of the textbook. For now, use the built-in `fixed_quad`. Specify `n=` the number of points. Docs [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.fixed_quad.html#scipy.integrate.fixed_quad)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2-point Gauss quadrature: 1.8225777778\n",
      "3-point Gauss quadrature: 1.6405333333\n",
      "4-point Gauss quadrature: 1.6405333333\n",
      "5-point Gauss quadrature: 1.6405333333\n",
      "6-point Gauss quadrature: 1.6405333333\n",
      "7-point Gauss quadrature: 1.6405333333\n",
      "8-point Gauss quadrature: 1.6405333333\n",
      "9-point Gauss quadrature: 1.6405333333\n"
     ]
    }
   ],
   "source": [
    "for n in range(2,10):\n",
    "    I,_ = integrate.fixed_quad(f, a, b, n=n)\n",
    "    print('%d-point Gauss quadrature: %0.10f' % (n,I))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our 2-point function matches this answer. Even though the 2-point estimate is not very accurate, the estimates with 3 and higher points are very accurate. \n",
    "\n",
    "Notice that 3-point Gauss should exactly integrate a $2N+1=5$th order polynomial, so in this case the answer is perfect. For a general $f(x)$ that is not a polynomial, it wouldn't work quite so well."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Compare accuracy to reach `tol=1e-4`\n",
    "- Romberg: 9 function evaluations\n",
    "- Adaptive: 21 function evaluations (with built-in function, ours was 129)\n",
    "- Gauss: 3 function evaluations (special case for 5th order polynomial)\n",
    "\n",
    "How do these compare to the trapezoid rule?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trapezoid: I=1.6405236, num_evals = 513\n"
     ]
    }
   ],
   "source": [
    "I_old = trapezoid(f, a, b, 1) # first estimate with one segment\n",
    "\n",
    "ns = 2 ** np.arange(1, 12) # 2, 4, 8, etc segments\n",
    "\n",
    "for n in ns:\n",
    "    I = trapezoid(f, a, b, n)\n",
    "    \n",
    "    if abs(I - I_old) < 1e-4:\n",
    "        break\n",
    "    else:\n",
    "        I_old = I\n",
    "\n",
    "\n",
    "print('Trapezoid: I=%0.7f, num_evals = %d' % (I, n+1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We would need over 500 function evaluations with the trapezoid rule to reach a tolerance of `1e-4`. This can be obtained with fewer than 10 function evaluations using some of these more advanced integration methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.12.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}