chapter5-1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a40cea43",
   "metadata": {},
   "source": [
    "# 5.1 Newton-Cotes formulas"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27dde733",
   "metadata": {},
   "source": [
    "The idea is to construct the polynomial interpolant $P(x)$ and compute\n",
    "$\\int_{a}^{b}P(x)dx$ as an approximation to $\\int_{a}^{b}f(x)dx$. Given nodes $x_{0},x_{1},...,x_{n}$,\n",
    "the Lagrange form of the interpolant is\n",
    "\\begin{equation*}\n",
    "P_{n}(x)=\\sum_{i=0}^{n}f(x_{i})l_{i}(x)\n",
    "\\end{equation*}\n",
    "and from the interpolation error formula Theorem 3 Chapter 4-1, we have\n",
    "\\begin{equation*}\n",
    "f(x)=P_{n}(x)+(x-x_{0})\\cdots(x-x_{n})\\frac{f^{(n+1)}(\\xi(x))}{(n+1)!},\n",
    "\\end{equation*}\n",
    "where $\\xi(x)\\in[a,b]$. (We have written $\\xi(x)$ instead of $\\xi$ to emphasize that $\\xi$ depends on the value of $x$.)\n",
    "\n",
    "Taking the integral of both sides yields\n",
    "\\begin{equation}\\label{eq:newtoncotes}\n",
    "\\int_{a}^{b}f(x)dx  = \\underbrace{\\int_{a}^{b}P_{n}(x)dx}_{\\text{quadrature rule}}+\\underbrace{\\frac{1}{(n+1)!}\\int_{a}^{b}\\prod_{i=0}^{n}(x-x_{i})f^{(n+1)}(\\xi(x))dx.}_{\\text{error term}}\n",
    "\\end{equation}\n",
    "\n",
    "The first integral on the right-hand side gives the quadrature rule:\n",
    "\\begin{equation*}\n",
    "\\int_{a}^{b}P_{n}(x)dx=\\int_a^b \\left( \\sum_{i=0}^{n}f(x_{i})l_{i}(x) \\right) dx =\\sum_{i=0}^{n}\\left(\\underbrace{\\int_{a}^{b}l_{i}(x)dx}_{w_i}\\right) f(x_{i})=\\sum_{i=0}^n w_i f(x_i),\n",
    "\\end{equation*}\n",
    "and the second integral gives the error term.\n",
    "\n",
    "We obtain different quadrature rules by taking different nodes, or number of\n",
    "nodes. The following result is useful in the theoretical analysis\n",
    "of Newton-Cotes formulas."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f3f9a14",
   "metadata": {},
   "source": [
    "\\begin{theorem}[Weighted mean value theorem for integrals]\\label{theo:chap4_thm65}\n",
    "Suppose $f\\in C^0[a,b],$ the Riemann integral of $g(x)$ exists on\n",
    "$[a,b],$ and $g(x)$ does not change sign on $[a,b].$ Then there\n",
    "exists $\\xi\\in(a,b)$ with $\\int_{a}^{b}f(x)g(x)dx=f(\\xi)\\int_{a}^{b}g(x)dx.$\n",
    "\\end{theorem}\n",
    "\n",
    "\n",
    "Two well-known numerical quadrature rules, trapezoidal rule and Simpson's rule, are examples of Newton-Cotes formulas:\n",
    "\n",
    "- **Trapezoidal rule**\n",
    "\n",
    "Let $f\\in C^2[a,b]$. Take two nodes, $x_{0}=a, x_{1}=b,$ and use the linear\n",
    "Lagrange polynomial\n",
    "\\begin{equation*}\n",
    "P_{1}(x)=\\frac{x-x_{1}}{x_{0}-x_{1}}f(x_{0})+\\frac{x-x_{0}}{x_{1}-x_{0}}f(x_{1})\n",
    "\\end{equation*}\n",
    "to estimate $f(x)$. Substitute $n=1$ in Equation \\eqref{eq:newtoncotes} to get\n",
    "\\begin{equation*}\n",
    "\\int_{a}^{b}f(x)dx=\\int_{a}^{b}P_{1}(x)dx+\\frac{1}{2}\\int_{a}^{b}\\prod_{i=0}^{1}(x-x_{i})f''(\\xi(x))dx,\n",
    "\\end{equation*}\n",
    "and then substitute for $P_1$ to obtain\n",
    "\\begin{align*}\n",
    "\\int_{a}^{b}f(x)dx & =\\int_{a}^{b}\\frac{x-x_{1}}{x_{0}-x_{1}}f(x_{0})dx+\\int_{a}^{b}\\frac{x-x_{0}}{x_{1}-x_{0}}f(x_{1})dx +\\frac{1}{2}\\int_{a}^{b}(x-x_{0})(x-x_{1})f''(\\xi(x))dx.\n",
    "\\end{align*}\n",
    "The first two integrals on the right-hand side can be evaluated easily: the first integral is $\\frac{h}{2}f(x_{0})$ and the second one is\n",
    "$\\frac{h}{2}f(x_{1})$, where $h=x_1-x_0=b-a$. Let's evaluate\n",
    "\\begin{equation*}\n",
    "\\int_{a}^{b}(x-x_{0})(x-x_{1})f''(\\xi(x))dx.\n",
    "\\end{equation*}\n",
    "We will use Theorem \\ref{theo:chap4_thm65} for this computation. Note that the function $(x-x_0)(x-x_1)=(x-a)(x-b)$ does not change sign on the interval $[a,b]$ and it is integrable: so this function serves the role of $g(x)$ in Theorem \\ref{theo:chap4_thm65}. The other term, $f''(\\xi(x))$, serves the role of $f(x)$. Applying the theorem, we get\n",
    "\\begin{equation*}\n",
    "\\int_{a}^{b}(x-a)(x-b)f''(\\xi(x))dx=f''(\\xi) \\int_a^b (x-a)(x-b)dx\n",
    "\\end{equation*}\n",
    "where we kept the same notation $\\xi$, somewhat inappropriately, as we moved $f''(\\xi(x))$ from inside to the outside of the integral. Finally, observe that\n",
    "\\begin{equation*}\n",
    "\\int_a^b (x-a)(x-b)dx=\\frac{(a-b)^3}{6}=\\frac{-h^3}{6},\n",
    "\\end{equation*}\n",
    "where $h=b-a$. Putting all the pieces together, we have obtained\n",
    "\\begin{equation*}\n",
    "\\int_{a}^{b}f(x)dx=\\frac{h}{2}\\left[f(x_{0})+f(x_{1})\\right]-\\frac{h^{3}}{12}f''(\\xi).\n",
    "\\end{equation*}\n",
    "\n",
    "- **Simpson's rule**\n",
    "\n",
    "Let $f\\in C^4[a,b]$. Take three equally-spaced nodes, $x_{0}=a$, $x_{1}=a+h$, $x_{2}=b,$\n",
    "where $h=\\frac{b-a}{2},$ and use the second degree Lagrange polynomial\n",
    "\\begin{equation*}\n",
    "P_{2}(x)=\\frac{(x-x_{1})(x-x_{2})}{(x_{0}-x_{1})(x_{0}-x_{2})}f(x_{0})+\\frac{(x-x_{0})(x-x_{2})}{(x_{1}-x_{0})(x_{1}-x_{2})}f(x_{1})+\\frac{(x-x_{0})(x-x_{1})}{(x_{2}-x_{0})(x_{2}-x_{1})}f(x_{2})\n",
    "\\end{equation*}\n",
    "to estimate $f(x)$. Substitute $n=2$ in Equation \\eqref{eq:newtoncotes} to get\n",
    "\\begin{equation*}\n",
    "\\int_{a}^{b}f(x)dx=\\int_{a}^{b}P_{2}(x)dx+\\frac{1}{3!}\\int_{a}^{b}\\prod_{i=0}^{2}(x-x_{i})f^{(3)}(\\xi(x))dx,\n",
    "\\end{equation*}\n",
    "and then substitute for $P_2$ to obtain\n",
    "\\begin{align*}\n",
    "\\int_{a}^{b}f(x)dx & =\\int_{a}^{b}\\frac{(x-x_{1})(x-x_{2})}{(x_{0}-x_{1})(x_{0}-x_{2})}f(x_{0})dx+\\int_{a}^{b}\\frac{(x-x_{0})(x-x_{2})}{(x_{1}-x_{0})(x_{1}-x_{2})}f(x_{1})dx+\\\\\n",
    " & +\\int_{a}^{b}\\frac{(x-x_{0})(x-x_{1})}{(x_{2}-x_{0})(x_{2}-x_{1})}f(x_{2})dx+\\frac{1}{6}\\int_{a}^{b}(x-x_{0})(x-x_{1})(x-x_{2})f^{(3)}(\\xi(x))dx.\n",
    "\\end{align*}\n",
    "The sum of the first three integrals on the right-hand side simplify as: $\\frac{h}{3}\\left[f(x_{0})+4f(x_{1})+f(x_{2})\\right]$. The last integral cannot be evaluated using Theorem \\ref{theo:chap4_thm65} directly, like in the trapezoidal rule, since the function $(x-x_0)(x-x_1)(x-x_2)$ changes sign on $[a,b]$. However, a clever application of integration by parts transforms the integral to an integral where  Theorem \\ref{theo:chap4_thm65} is applicable (see Atkinson [1989] for details), and the integral simplifies as $-\\frac{h^{5}}{90}f^{(4)}(\\xi)$ for some $\\xi \\in (a,b)$. In summary, we obtain\n",
    "\\begin{equation*}\n",
    "\\int_{a}^{b}f(x)dx=\\frac{h}{3}\\left[f(x_{0})+4f(x_{1})+f(x_{2})\\right]-\\frac{h^{5}}{90}f^{(4)}(\\xi)\n",
    "\\end{equation*}\n",
    "where $\\xi\\in(a,b)$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86f06da6",
   "metadata": {},
   "source": [
    "**Exercise 5.1-1**\n",
    "\n",
    "Prove that the sum of the weights in Newton-Cotes rules is $b-a$, for any $n$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12cec0fd",
   "metadata": {},
   "source": [
    "\\begin{definition}\\label{def:chap4_def66}\n",
    "The degree of accuracy, or precision, of a quadrature formula is the\n",
    "largest positive integer $n$ such that the formula is exact for $f(x)=x^{k},$\n",
    "when $k=0,1,...,n,$ or equivalently, for any polynomial of degree\n",
    "less than or equal to $n$.\n",
    "\\end{definition}\n",
    "\n",
    "Observe that the trapezoidal and Simpson's rules have degrees of accuracy\n",
    "of one and three. These two rules are examples of closed Newton-Cotes\n",
    "formulas; closed refers to the fact that the end points $a,b$ of\n",
    "the interval are used as nodes in the quadrature rule. Here is the\n",
    "general definition."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab9bc9d9",
   "metadata": {},
   "source": [
    "\\begin{definition}\\label{def:chap4_def67}\n",
    "The $(n+1)$-point closed Newton-Cotes formula uses nodes $x_{i}=x_{0}+ih$,\n",
    "for $i=0,1,...,n,$ where $x_{0}=a,x_{n}=b,h=\\frac{b-a}{n},$ and\n",
    "\\begin{equation*}\n",
    "w_{i}=\\int_{x_{0}}^{x_{n}}l_{i}(x)dx=\\int_{x_{0}}^{x_{n}}\\prod_{j=0,j\\neq i}^{n}\\frac{x-x_{j}}{x_{i}-x_{j}}dx.\n",
    "\\end{equation*}\n",
    "\\end{definition}\n",
    "\n",
    "The following theorem provides an error formula for the closed Newton-Cotes formula. A proof can be found in Keller [1966]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7b03487",
   "metadata": {},
   "source": [
    "\\begin{theorem}\\label{theo:chap4_thm68}\n",
    "For the $(n+1)$-point closed Newton-Cotes formula, we have:\n",
    "\n",
    "- if $n$ is even and $f\\in C^{n+2}[a,b]$\n",
    "\\begin{equation*}\n",
    "\\int_{a}^{b}f(x)dx=\\sum_{i=0}^{n}w_{i}f(x_{i})+\\frac{h^{n+3}f^{(n+2)}(\\xi)}{(n+2)!}\\int_{0}^{n}t^{2}(t-1)\\cdots(t-n)dt,\n",
    "\\end{equation*}\n",
    "- if $n$ is odd and $f\\in C^{n+1}[a,b]$\n",
    "\\begin{equation*}\n",
    "\\int_{a}^{b}f(x)dx=\\sum_{i=0}^{n}w_{i}f(x_{i})+\\frac{h^{n+2}f^{(n+1)}(\\xi)}{(n+1)!}\\int_{0}^{n}t(t-1)\\cdots(t-n)dt\n",
    "\\end{equation*}\n",
    "where $\\xi$ is some number in $(a,b)$.\n",
    "\\end{theorem}\n",
    "\n",
    "Some well known examples of closed Newton-Cotes formulas are the trapezoidal\n",
    "rule ($n=1)$, Simpson's rule ($n=2)$, and Simpson's three-eighth\n",
    "rule ($n=3)$. Observe that in the $(n+1)$-point closed Newton-Cotes\n",
    "formula, if $n$ is even, then the degree of accuracy is $(n+1$),\n",
    "although the interpolating polynomial is of degree $n$. \n",
    "The open Newton-Cotes formulas exclude the end points of the interval."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd436d8f",
   "metadata": {},
   "source": [
    "\\begin{definition}\\label{def:chap4_def69}\n",
    "The $(n+1)$-point open Newton-Cotes formula uses nodes $x_{i}=x_{0}+ih$,\n",
    "for $i=0,1,...,n,$ where $x_{0}=a+h,x_{n}=b-h,h=\\frac{b-a}{n+2}$\n",
    "and\n",
    "\\begin{equation*}\n",
    "w_{i}=\\int_{a}^{b}l_{i}(x)dx=\\int_{a}^{b}\\prod_{j=0,j\\neq i}^{n}\\frac{x-x_{j}}{x_{i}-x_{j}}dx.\n",
    "\\end{equation*}\n",
    "We put $a=x_{-1}$ and $b=x_{n+1}$.\n",
    "\\end{definition}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "014e0458",
   "metadata": {},
   "source": [
    "The error formula for the open Newton-Cotes formula is given next; for a proof see Isaacson and Keller [1966].\n",
    "\\begin{theorem}\\label{theo:chap4_thm70}\n",
    "For the $(n+1)$-point open Newton-Cotes formula, we have:\n",
    "\n",
    "- if $n$ is even and $f\\in C^{n+2}[a,b]$\n",
    "\\begin{equation*}\n",
    "\\int_{a}^{b}f(x)dx=\\sum_{i=0}^{n}w_{i}f(x_{i})+\\frac{h^{n+3}f^{(n+2)}(\\xi)}{(n+2)!}\\int_{-1}^{n+1}t^{2}(t-1)\\cdots(t-n)dt,\n",
    "\\end{equation*}\n",
    "- if $n$ is odd and $f\\in C^{n+1}[a,b]$\n",
    "\\begin{equation*}\n",
    "\\int_{a}^{b}f(x)dx=\\sum_{i=0}^{n}w_{i}f(x_{i})+\\frac{h^{n+2}f^{(n+1)}(\\xi)}{(n+1)!}\\int_{-1}^{n+1}t(t-1)\\cdots(t-n)dt,\n",
    "\\end{equation*}\n",
    "where $\\xi$ is some number in $(a,b)$.\n",
    "\\end{theorem}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b57aa1c9",
   "metadata": {},
   "source": [
    "The well-known midpoint rule is an example of open Newton-Cotes formula:\n",
    "\n",
    "- **Midpoint rule**\n",
    "\n",
    "Take one node, $x_0=a+h$, which corresponds to $n=0$ in the above theorem to obtain\n",
    "\\begin{equation*}\n",
    "\\int_{a}^{b}f(x)dx=2hf(x_{0})+\\frac{h^{3}f''(\\xi)}{3}\n",
    "\\end{equation*}\n",
    "where $h=(b-a)/2.$ This rule interpolates $f$ by\n",
    "a constant (the values of $f$ at the midpoint), that is, a polynomial\n",
    "of degree 0, but it has degree of accuracy 1."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98752739",
   "metadata": {},
   "source": [
    "\\begin{remark}\\label{rem:chap4_rem71}\n",
    "Both closed and open Newton-Cotes formulas using odd number of nodes\n",
    "($n$ is even), gain an extra degree of accuracy beyond that of the\n",
    "polynomial interpolant on which they are based. This is due to cancellation\n",
    "of positive and negative error.\n",
    "\n",
    "There are some drawbacks of Newton-Cotes formulas:\n",
    "\n",
    "- In general, these rules are not of the highest degree of accuracy possible\n",
    "for the number of nodes used.\n",
    "- The use of large number of equally spaced nodes may incur the erratic\n",
    "behavior associated with high-degree polynomial interpolation. Weights\n",
    "for a high-order rule may be negative, potentially leading to loss\n",
    "of significance errors.\n",
    "- Let $I_{n}$ denote the Newton-Cotes estimate of an integral based\n",
    "on $n$ nodes. $I_{n}$ may not converge to the true integral as $n\\rightarrow\\infty$\n",
    "for perfectly well-behaved integrands.\n",
    "\\end{remark}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32aa1337",
   "metadata": {},
   "source": [
    "\\begin{example}\\label{example:chap4_exa72}\n",
    "Estimate $\\int_{0.5}^1 x^x dx$ using the midpoint, trapezoidal, and Simpson's rules.\n",
    "\\end{example}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "551ffc04",
   "metadata": {},
   "source": [
    "**Solution.**\n",
    "\n",
    "Let $f(x)=x^x$.\n",
    "The midpoint estimate for the integral is $2hf(x_{0})$ where $h=(b-a)/2=1/4$ and $x_{0}=0.75.$\n",
    "Then the midpoint estimate, using 6 digits, is $f(0.75)/2=0.805927/2=0.402964$.\n",
    "The trapezoidal estimate is $\\frac{h}{2}\\left[f(0.5)+f(1)\\right]$\n",
    "where $h=1/2$, which results in $1.707107/4=0.426777.$\n",
    "Finally, for Simpson's rule, $h=(b-a)/2=1/4,$ and thus the estimate is\n",
    "\\begin{equation*}\n",
    "\\frac{h}{3}\\left[f(0.5)+4f(0.75)+f(1)\\right]=\\frac{1}{12}\\left[0.707107+4(0.805927)+1\\right]=0.410901.\n",
    "\\end{equation*}\n",
    "Here is a summary of the results:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{array}{ccc}\n",
    "\\text{Midpoint} & \\text{Trapezoidal} & \\text{Simpson's} \\\\\n",
    "\\hline\n",
    "0.402964 & 0.426777 & 0.410901  \\\\\n",
    "\\end{array}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed6f954a",
   "metadata": {},
   "source": [
    "\\begin{example}\\label{example:chap4_exa73}\n",
    "Find the constants $c_{0},c_{1},x_{1}$ so that the quadrature formula\n",
    "\\begin{equation*}\n",
    "\\int_{0}^{1}f(x)dx=c_{0}f(0)+c_{1}f(x_{1})\n",
    "\\end{equation*}\n",
    "has the highest possible degree of accuracy.\n",
    "\\end{example}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb37998c",
   "metadata": {},
   "source": [
    "**Solution.**\n",
    "\n",
    "We will find how many of the polynomials $1,x,x^{2},...$ the rule\n",
    "can integrate exactly. If $p(x)=1,$ then\n",
    "\\begin{equation*}\n",
    "\\int_{0}^{1}p(x)dx=c_{0}p(0)+c_{1}p(x_{1})\\Rightarrow1=c_{0}+c_{1}.\n",
    "\\end{equation*}\n",
    "If $p(x)=x,$ we get\n",
    "\\begin{equation*}\n",
    "\\int_{0}^{1}p(x)dx=c_{0}p(0)+c_{1}p(x_{1})\\Rightarrow\\frac{1}{2}=c_{1}x_{1}\n",
    "\\end{equation*}\n",
    "and $p(x)=x^{2}$ implies\n",
    "\\begin{equation*}\n",
    "\\int_{0}^{1}p(x)dx=c_{0}p(0)+c_{1}p(x_{1})\\Rightarrow\\frac{1}{3}=c_{1}x_{1}^{2}.\n",
    "\\end{equation*}\n",
    "We have three unknowns and three equations, so we have to stop here.\n",
    "Solving the three equations we get: $c_{0}=1/4,c_{1}=3/4,x_{1}=2/3.$\n",
    "So the quadrature rule is of precision two and it is:\n",
    "\\begin{equation*}\n",
    "\\int_{0}^{1}f(x)dx=\\frac{1}{4}f(0)+\\frac{3}{4}f(\\frac{2}{3}).\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67d1a162",
   "metadata": {},
   "source": [
    "**Exercise 5.1-2**\n",
    "\n",
    "Find $c_{0},c_{1},c_{2}$ so that the quadrature rule $\\int_{-1}^{1}f(x)dx=c_{0}f(-1)+c_{1}f(0)+c_{2}f(1)$\n",
    "has degree of accuracy 2."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "013def30",
   "metadata": {},
   "source": [
    "**References**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb4eff25",
   "metadata": {},
   "source": [
    "K.E. Atkinson. An Introduction to Numerical Analysis. John Wiley & Sons, 2nd edition, 1989.\n",
    "\n",
    "Isaacson, E. and H.B. Keller. Analysis of Numerical Methods. John Wiley & Sons, 1966."
   ]
  }
 ],
 "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.7"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}