{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Definite Integrals, Part 4: Romberg Integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Updated on Thursday, April 1** adding a note about an error estimate that can be used in an algorithm that computes to a specified error tolerance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**References:**\n",
    "\n",
    "- Section 5.3 *Romberg Integration* of [Sauer](../references.html#Sauer)\n",
    "- Section 4.5 *Romberg Integration* of [Burden&Faires](../references.html#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",
    "$$\n",
    "R_{0,0} = T_1 = \\frac{f(a) + f(b)}{2}(b-a),\n",
    "$$\n",
    "\n",
    "because the most efficient way to get the other values is recursively, with\n",
    "\n",
    "$$\n",
    "T_{2n} = \\frac{T_n + M_n}{2}\n",
    "$$\n",
    "\n",
    "where $M_n$ is the composite midpoint rule,\n",
    "\n",
    "$$\n",
    "M_n = h \\sum_{k=1}^n f(a + (k-1/2)h), \\quad h = \\frac{b-a}{n}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the \"Romberg notation\" $R_{i,0} = T_{2^i}$ so that $n=2^i$, this becomes\n",
    "\n",
    "$$R_{i+1,0} = \\frac{R_{i,0} + M_{2^i}}{2}$$\n",
    "\n",
    "or more conveniently for below,\n",
    "\n",
    "$$R_{i,0} = \\frac{R_{i-1,0} + M_{2^{i-1}}}{2}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extrapolation is then done with the formula\n",
    "\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",
    "\n",
    "which can also be expressed as\n",
    "\n",
    "$$R_{i,j} = R_{i,j-1} + E_{i,j-1}$$\n",
    "\n",
    "where\n",
    "\n",
    "$$E_{i,j-1} = \\frac{R_{i,j-1} - R_{i-1,j-1}}{4^j - 1}$$\n",
    "\n",
    "is an error estimate.\n",
    "This can be used as a stopping condition in an algorithm that computes to a specified error tolerance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Aside:** the first level of extrapolation gives the Composite Simpson's Rule:\n",
    "\n",
    "$$\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": [
    "## A first 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": [
    "$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`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "This work is licensed under [Creative Commons Attribution-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-sa/4.0/)"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}