{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Definite Integrals, Part 4: Romberg Integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**References:**\n",
    "\n",
    "- Section 5.3 *Romberg Integration* of {cite}`Sauer`.\n",
    "- Section 4.5 *Romberg Integration* of {cite}`Burden-Faires`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "Romberg Integration is based on repeated Richardson extrapolalation from the composite trapezoidal rule, starting with one interval and repeatedly doubling.\n",
    "Our notation starts with\n",
    "\n",
    "$$ R_{i,0} = T_{2^i}, \\quad i=0, 1, 2, \\dots $$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "T_n = \\left(\\frac{f(a)}{2} + \\sum_{k=1}^{n-1} f(a + kh) + \\frac{f(b)}{2}\\right) h, \\quad h = \\frac{b-a}{n}\n",
    "$$\n",
    "\n",
    "and the second index will indicate the number of extapolation steps done (none so far!)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually we only need this $T_n$ formula for the single trapezoidal rule, to get\n",
    "\n",
    "$$ R_{0,0} = T_1 = \\frac{f(a) + f(b)}{2}(b-a), $$\n",
    "\n",
    "because the most efficient way to get the other values is recursively, with\n",
    "\n",
    "$$ T_{2n} = \\frac{T_n + M_n}{2} $$\n",
    "\n",
    "where $M_n$ is the composite midpoint rule,\n",
    "\n",
    "$$ M_n = h \\sum_{k=1}^n f(a + (k-1/2)h), \\quad h = \\frac{b-a}{n} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extrapolation is then done with the formula\n",
    "\n",
    "$$ R_{i,j} = \\frac{4^j R_{i,j-1} - R_{i-1,j-1}}{4^j - 1}, \\quad j = 1, 2, \\dots, i $$\n",
    "\n",
    "which can also be expressed as\n",
    "\n",
    "$$ R_{i,j} = R_{i,j-1} + E_{i,j-1}, \\;\\text{where}\\; E_{i,j-1} = \\frac{R_{i,j-1} - R_{i-1,j-1}}{4^j - 1}\\quad\\text{is an error estimate.} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "R_{i,1} = S_{2n} = \\frac{4 T_{2 n} - T_{n}}{4 - 1}, \\; n=2^{i-1}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An algorithm, in pseudocode\n",
    "\n",
    "The above can now be arranged into a basic algorithm.\n",
    "It does a fixed number $M$ of levels of extrapolation so using $2^M$ intervals;\n",
    "a refinement would be to use the above error estimate $E_{i,j-1}$ as the basis for a stopping condition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{prf:algorithm} Romberg Integration\n",
    ":label: romberg-integration\n",
    "\n",
    "$n \\leftarrow 1$\n",
    "<br>\n",
    "$h \\leftarrow b-a$\n",
    "<br>\n",
    "$\\displaystyle R_{0,0} = \\frac{f(a) + f(b)}{2} h$\n",
    "<br>\n",
    "`for i from 1 to M:`\n",
    "<br>\n",
    "$\\quad$ $R_{i,0} = \\left( R_{i-1,0} + h \\sum_{k=1}^n f(a + (i - 1/2)h) \\right)/2$\n",
    "<br>\n",
    "$\\quad$ `for j from 1 to i:`\n",
    "<br>\n",
    "$\\quad\\quad \\displaystyle R_{i,j} = \\frac{4^j R_{i,j-1} - R_{i-1,j-1}}{4^j - 1}$\n",
    "<br>\n",
    "$\\quad$ `end for`\n",
    "<br>\n",
    "$\\quad n \\leftarrow 2 n$\n",
    "<br>\n",
    "$\\quad h \\leftarrow h/2$\n",
    "<br>\n",
    "`end for`\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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}