{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# Numerical integration \n", "\n", " \n", "**Anne Kværnø**\n", "\n", "Date: **Oct 14, 2018**\n", "\n", "# Introduction\n", "Given the finite integral" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "I[f](a,b) = \\int_a^b f(x) dx.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A *numerical quadrature* or a *quadrature rule* is a formula for approximating such integrals. Quadratures are usually of the form" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "Q[f](a,b) = \\sum_{i=0}^n w_i f(x_i),\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $x_i$, $w_i$ for $i=0,1,\\dotsc,n$ are respectively the *nodes* and the *weights* of the quadrature rule. \n", "If the function $f$ is given from the context, we will for simplicity denote the integral and the quadrature simply as $I(a,b)$ and $Q(a,b)$.\n", "\n", "The [trapezoidal rule, the midpoint rule and Simpson's rule](https://wiki.math.ntnu.no/tma4100/tema/numerics?&#numerisk_integrasjon) known from previous courses are all examples of numerical quadratures. \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "In this note we will see how quadrature rules can be constructed from integration of interpolation polynomials.\n", "We will demonstrate how to do error analysis and how to find error estimates. We will also demonstrate how the integration interval can be automatically partitioned in subintervals, so called *adaptive* integration. \n", "\n", "In the sequel, we will use material from *Preliminaries*, section 3.2, 4 and 5. \n", "\n", "# Quadrature based on polynomial interpolation.\n", "This section relies on the content of the note on polynomial interpolation, in particular the section on Lagrange polynomials. \n", "\n", "Choose $n+1$ distinct nodes $x_i$, $i=0,\\dotsc,n$ in the interval $[a,b]$, and let $p_n(x)$ be the interpolation polynomial satisfying\n", "the interpolation condition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p_n(x_i) = f(x_i), \\qquad i=0,1,\\dotsc.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will then use $\\int_a^b p_n(x)dx$ as an approximation to $\\int_a^b f(x)dx$. By using the Lagrange form of the polynomial" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "p_n(x) = \\sum_{i=0}^n f(x_i) \\ell_i(x)\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with the cardinal functions $\\ell_i(x)$ given by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\ell_i(x) = \\prod_{j=0,j\\not=i}^n \\frac{x-x_j}{x_i-x_j},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "the following quadrature formula is obtained" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "Q[f](a,b) = \\int_a^b p_n(x)dx\n", " = \\sum_{i=0}^n f(x_i) \\int_a^b \\ell_i(x) dx = \\sum_{i=0}^n w_i f(x_i) = Q(a,b).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The weights in the quadrature is simply the integral of the cardinal functions over the interval.\n", "\n", "Let us derive two schemes for integration over the interval $[0,1]$, and apply them to the integral" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "I(0,1) = \\int_0^1 \\cos\\left(\\frac{\\pi}{2}x\\right) = \\frac{2}{\\pi} = 0.636619\\dotsc.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Example 1:**\n", "Let $x_0=0$ and $x_1=1$. The cardinal functions and thus the weights are given by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\\begin{align*}\n", "\\ell_0(x) &= 1-x, & w_0 &= \\int_0^1(1-x)dx = 1/2 \\\\ \n", "\\ell_1(x) &= x, & w_1 &= \\int_0^1 xdx = 1/2\n", "\\end{align*}\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and the corresponding quadrature rule, better known as the trapezoidal rule and usually denoted by $T$, is given by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "T(0,1) = \\frac{1}{2} \\left[ f(0) + f(1) \\right].\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This formula applied to the function $f(x)=\\cos(\\pi x/2)$ gives" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "T(0,1) = \\frac{1}{2}\\left[ \\cos(0) + \\cos\\left(\\frac{\\pi}{2}\\right)\\right] = \\frac{1}{2},\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and the error is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "I(0,1) - T(0,1) = \\frac{2}{\\pi}-\\frac{1}{2} = 0.138\\dotsc\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Example 2:**\n", "Let $x_0=1/2 + \\sqrt{3}/6$ and $x_1 = 1/2 - \\sqrt{3}/6$. Then" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\\begin{align*}\n", "\\ell_0(x) &= -\\sqrt{3}x + \\frac{1+\\sqrt{3}}{2}, & w_0 &= \\int_0^1 \\ell_0(x)dx= 1/2, \\\\ \n", "\\ell_1(x) &= \\sqrt{3}x + \\frac{1-\\sqrt{3}}{2}, & w_1 &= \\int_0^1 \\ell_1(x)dx = 1/2.\n", "\\end{align*}\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The quadrature rule is" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "Q(0,1) = \\frac{1}{2}\\left[f\\left(\\frac{1}{2}-\\frac{\\sqrt{3}}{6}\\right) + \n", "f\\left(\\frac{1}{2}+\\frac{\\sqrt{3}}{6}\\right) \\right].\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And this quadrature applied to $f(x)=\\cos(\\pi x/2)$ is given by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "Q(0,1) = \\frac{1}{2}\\left[\\cos\\left(\\frac{\\pi}{2}x_0\\right) + \\cos\\left(\\frac{\\pi}{2}x_1\\right) \\right] = 0.635647\\dotsc\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with an error" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "I(0,1)-Q(0,1) = 9.72\\dotsc \\cdot 10^{-4}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So the choice of nodes clearly matters!\n", "\n", "Before concluding this section, let us present simple indication on the quality of a method: \n", "\n", "**Definition: The degree of precision.**\n", "\n", "A numerical quadrature has degree of precision $d$ if \n", "$Q[p](a,b) = I[p](a,b)$ for alle $p \\in \\mathbb{P}_d$.\n", "\n", "\n", "\n", "\n", "Since both integrals and quadratures are linear in the integrand $f$, the degree of precision is $d$ if" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\\begin{align*}\n", "I[x^j](a,b) &= Q[x^j](a,b), \\qquad j=0,1,\\dotsc, d, \\\\ \n", "I[x^{d+1}](a,b) &\\not= Q[x^{d+1}](a,b)\n", "\\end{align*}\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All quadratures constructed from Lagrange interpolation polynomials in $n+1$ distinct nodes will automatically be of precision at least $n$. \n", "\n", "It is left to the reader to show that\n", "the trapezoidal rule from Example 1 is of precision 1, the formula from Example 2 of precision 3.\n", "\n", "# Construction of numerical quadratures\n", "In the following, you will learn the steps on how to construct realistic\n", "algorithms for numerical integration, similar to those used in software like\n", "Matlab of SciPy. The steps are: \n", "**Construction.**\n", "\n", "1. Choose $n+1$ distinct nodes on a standard interval $[-1,1]$. \n", "\n", "2. Let $p_n(x)$ be the polynomial interpolating some general function $f$ in the nodes, and let the $Q[f](-1,1)=I[p_n](-1,1)$. \n", "\n", "3. Transfer the formula $Q$ from $[-1,1]$ to some interval $[a,b]$.\n", "\n", "4. Find the composite formula, by dividing the interval $[a,b]$ into subintervals and applying the quadrature formula on each subinterval.\n", "\n", "5. Find an expression for the error $E[f](a,b) = I[f](a,b)-Q[f](a,b)$. \n", "\n", "6. Find an expression for an estimate of the error, and use this to create an adaptive algorithm.\n", "\n", "\n", "\n", "## Simpson's rule\n", "We will go through the steps above for one method, Simpson's formula. The strategy is quite generic, so it is more important to understand and remember how results are derived, not exactly what they are. The different algorithms will be implemented and tested, and theoretical results will be verified by numerical experiments. \n", "\n", "We will adopt the standard notation and denote this particular quadrature by $S[f](a,b)$. \n", "\n", "### The quadrature formula on the standard interval [-1,1]\n", "\n", "The quadrature rule is defined by the choice of nodes on a standard interval $[-1,1]$. For Simpson's rule, choose\n", "the nodes $t_0=-1$, $t_1=0$ and $t_2=1$. \n", "The corresponding cardinal\n", "functions are" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\ell_0 = \\frac{1}{2}(t^2-t), \\qquad\n", "\\ell_1(t) = 1-t^2, \\qquad\n", "\\ell_2(t) = \\frac{1}{2}(t^2+t).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "which gives the weights" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "w_0 = \\int_{-1}^1 \\ell_0(t)dt = \\frac{1}{3}, \\qquad\n", "w_1 = \\int_{-1}^1 \\ell_1(t)dt = \\frac{4}{3}, \\qquad\n", "w_2 = \\int_{-1}^1 \\ell_2(t)dt = \\frac{1}{3}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "such that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\int_{-1}^1 f(t) dt \\approx \\int_{-1}^1 p_2(t) dt = \\sum_{i=0}^2 w_i f(t_i) =\n", "\\frac{1}{3} \\left[\\; f(-1) + 4 f(0) + f(1) \\; \\right].\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Simpson's rule has degree of precision 3 (check it yourself). \n", "\n", "**Example 3:**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\int_{-1}^1 \\cos \\left( \\frac{\\pi t}{2}\\right)dt = \\frac{4}{\\pi}\n", "\\approx \\frac{1}{3}\\left[\\cos(-\\pi/2) + 4 \\cos(0)\n", "+ \\cos(\\pi/2) \\right]= \\frac{4}{3}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Transfer the integral and the quadrature to the interval $[a,b]$\n", "\n", "The integral and the quadrature is transferred to some arbitrary interval\n", "$[a,b]$ by the transformation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "x = \\frac{b-a}{2}t + \\frac{b+a}{2}, \\qquad \\text{so} \\qquad dx = \\frac{b-a}{2}dt.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Thus" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\int_a^b f(x)dx = \\frac{b-a}{2} \\int_{-1}^1 f\\left(\\frac{b-a}{2}t + \\frac{b+a}{2}\\right) dt\n", "\\approx \\frac{b-a}{6} \\left[\\;f(a)+4f\\left(\\frac{b+a}{2}\\right)+f(b)\\;\\right].\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Simpson's rule over the interval $[a,b]$ becomes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "S(a,b) = \\frac{b-a}{6}\\left[\\; f(a)+4f(c)+f(b)\\; \\right], \\qquad c=\\frac{b+a}{2}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Composite Simpson's rule\n", "\n", "Divide $[a,b]$ into $2m$ equal intervals of length \n", "$h = (b-a)/(2m)$. Let $x_j = a+jh$, $i=0,\\cdots,2m$, and apply Simpson's rule\n", "on each subinterval $[x_{2j}, x_{2j+2}]$. The result is:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\int_a^b f(x)dx = \\sum_{j=0}^{m-1} \\int_{x_{2j}}^{x_{2j+2}} f(x) dx\n", "\\approx \\sum_{j=0}^{m-1} S(x_{2j},x_{2j+2}) \n", "\\label{_auto1} \\tag{1}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", " = \\sum_{j=0}^{m-1} \\frac{h}{3}\n", "\\left[ f(x_{2j}) + 4 f(x_{2j+1})+ f(x_{2j+2}) \\right] \n", "\\label{_auto2} \\tag{2}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} =\n", "\\frac{h}{3} \\left[ f(x_0) + 4\\sum_{j=0}^{m-1}f(x_{2j+1}) + 2 \\sum_{j=1}^{m-1}f(x_{2j}) + f(x_{2m}) \\right]\n", "\\label{_auto3} \\tag{3}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will use the the notation $S_m(a,b)$ for the composite Simpson's rule on $m$ subintervals. \n", "\n", "### Implementation and testing\n", "\n", "It is now time to implement the composite Simpson's method, and see how well it\n", "works. \n", "Start by calling the necessary modules etc:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "from numpy import *\n", "from matplotlib.pyplot import *\n", "from math import factorial\n", "newparams = {'figure.figsize': (8.0, 4.0), 'axes.grid': True,\n", " 'lines.markersize': 8, 'lines.linewidth': 2,\n", " 'font.size': 14}\n", "rcParams.update(newparams)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def simpson(f, a, b, m=10):\n", "# Find an approximation to an integral by the composite Simpson's method:\n", "# Input: \n", "# f: integrand\n", "# a, b: integration interval\n", "# m: number of subintervals\n", "# Output: The approximation to the integral\n", " n = 2*m\n", " x_noder = linspace(a, b, n+1) # equidistributed nodes from a to b \n", " h = (b-a)/n # stepsize\n", " S1 = f(x_noder) + f(x_noder[n]) # S1 = f(x_0)+f(x_n)\n", " S2 = sum(f(x_noder[1:n:2])) # S2 = f(x_1)+f(x_3)+...+f(x_m)\n", " S3 = sum(f(x_noder[2:n-1:2])) # S3 = f(x_2)+f(x_4)+...+f(x_{m-1})\n", " S = h*(S1 + 4*S2 + 2*S3)/3\n", " return S" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Test if the code is correct. We know that Simpson's rule has precision 3, thus\n", "all third degree polynomials can be integrated exactly. Choose one such\n", "polynomial, find the exact integral, and compare. \n", "\n", "**Numerical experiment 1:**\n", "Apply the code on the integral, and compare with the exact result." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\int_{-1}^2(4x^3+x^2+2x-1)dx = 18.\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Numerical experiment 1\n", "def f(x): # Integrand\n", " return 4*x**3+x**2+2*x-1 \n", "a, b = -1, 2 # Integration interval\n", "I_exact = 18.0 # Exact value of the integral (for comparision)\n", "S = simpson(f, a, b, m=1) # Numerical solution, using m subintervals \n", "err = I_exact-S # Error\n", "print('I = {:.8f}, S = {:.8f}, error = {:.3e}'.format(I_exact, S, err))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Numerical experiment 2:**\n", "We will assume that the error decreases when the number of subintervals $m$\n", "increases. But how much? \n", "\n", "Apply the composite method on the integral (again with a known solution):" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\int_0^1 \\cos\\left(\\frac{\\pi x}{2}\\right )dx = \\frac{2}{\\pi}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use the function 'simpson' with $m=1,2,4,8,16$ and see how the error changes\n", "with $m$. Comment on the result." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Numerical experiment 2\n", "def f(x):\n", " return cos(0.5*pi*x)\n", "a, b = 0, 1\n", "I_exact = 2/pi\n", "for m in [1,2,4,8,16]:\n", " S = simpson(f, a, b, m=m) # Numerical solution, using m subintervals \n", " err = I_exact-S # Error\n", " if m == 1:\n", " print('m = {:3d}, error = {:.3e}'.format(m, err))\n", " else:\n", " print('m = {:3d}, error = {:.3e}, reduction factor = {:.3e}'.format(m, err, err/err_prev))\n", " err_prev=err" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "From the experiment we observe that the error is reduced by a factor\n", "approximately $0.0625 = 1/16$ whenever the number of subintervals increases with\n", "a factor 2. In the following, we will prove that this is in fact what can be\n", "expected. \n", "\n", "## Error analysis\n", "First we will find an expression for the error $E(a,b)=I(a,b)-S(a,b)$ over one\n", "interval $(a,b)$. This will then be used to find an expression for the composite\n", "formula. \n", "\n", "Let $c=(a+b)/2$ be the midpoint of the interval, and $h=(b-a)/2$ be the distance between $c$\n", "and the endpoints $a$ and $b$. Do a Taylor series expansion of the integrand $f$\n", "around the midpoint, and integrate each term in the series." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\int_a^b f(x)dx = \\int_{-h}^{h} f(c+s)ds =\n", "\\int_{-h}^h \\left( f(c) + sf'(c) + \\frac{1}{2}s^2 f''(c) + \\frac{1}{6} s^3 f'''(c) + \\frac{1}{24}s^4 f^{(4)}(c) + \\dotsm \\right) ds \n", "\\label{_auto4} \\tag{4}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "= 2h f(c) + \\frac{h^3}{3} f''(c) + \\frac{h^5}{60} f^{(4)}(c) + \\dotsm.\n", "\\label{_auto5} \\tag{5}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly, do a Taylor series expansion of the quadrature $S(a,b)$ around c:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "S(a,b) = \\frac{h}{3}\\left( f(c-h)+4f(c)+f(c+h) \\right) \n", "\\label{_auto6} \\tag{6}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", " = \\frac{h}{3}\\left( f(c) - hf'(c) + \\frac{1}{2}h^2 f''(c) - \\frac{1}{6} h^3 f'''(c) + \\frac{1}{24}h^4 f^{(4)}(c) + \\dotsm \\right. \n", "\\label{_auto7} \\tag{7}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", " \\qquad + 4f(c) \n", "\\label{_auto8} \\tag{8}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", " \\qquad + \\left. f(c) + hf'(c) + \\frac{1}{2}h^2 f''(c) + \\frac{1}{6} h^3 f'''(c) + \\frac{1}{24}h^5 f^{(4)}(c) + \\dotsm \\right) \n", "\\label{_auto9} \\tag{9}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", " = 2h f(c) + \\frac{h^3}{3} f''(c) + \\frac{h^5}{36} f^{(4)}(c) + \\dotsm\n", "\\label{_auto10} \\tag{10}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The series expansion of the error becomes:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E(a,b) = \\int_a^b f(x) dx - S(a,b) = -\\frac{h^5}{90} f^{(4)}(c) + \\cdots\n", " = - \\frac{(b-a)^5}{2^5 \\cdot 90} f^{(4)}(c) + \\dotsm,\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "using $h=(b-a)/2$. \n", "\n", "**NB!** By choosing to do the Taylor-expansions around the midpoint, every second\n", "term disappear thanks to symmetry. Choosing another point $\\hat{c}$ in the interval will\n", "give the same dominant error term (with $c$ replaced by $\\hat{c}$), but the\n", "calculations will be much more cumbersome.\n", "\n", "Usually, we will assume $h$ to be small, such that the first nonzero term in the\n", "series dominates the error, and the rest of the series can be\n", "ignored. It is however possible, but not trivial, to prove the following result:\n", "\n", "**Theorem: Error in Simpson's method.**\n", "\n", "Let $f(x) \\in C^{4}[a,b]$. There exist a $\\xi \\in (a,b)$ such that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E(a,b) = \\int_a^b f(x)dx - \\frac{b-a}{6} \\left[\\;f(a)+4f\\left(\\frac{b+a}{2}\\right)+f(b)\\;\\right] = -\\frac{(b-a)^5}{2880}f^{(4)}(\\xi).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**NB!**: Since $p^{(4)}(x)=0$ for all $p \\in \\mathbb{P}_3$ the degree of precisision is 3. \n", "\n", "Use the theorem to find an expression for the error in the composite Simpson's\n", "formula $S_m(a,b)$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "\\int_a^b f(x)dx - S_{m}(a,b) =\n", "\\sum_{j=0}^{m-1} \\left( \\int_{x_{2j}}^{x_{2j+2}} f(x)dx - \\frac{h}{3}\n", "\\left[ f(x_{2j}) + 4 f(x_{2j+1})+ f(x_{2j+2}) \\right] \\right) \n", "\\label{_auto11} \\tag{11}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", " = \\sum_{j=0}^{m-1} -\\frac{(2h)^5}{2880} f^{(4)}(\\xi_j)\n", "\\label{_auto12} \\tag{12}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $\\xi_{j} \\in (x_{2j}, x_{2j+2})$. We can then use the generalized mean\n", "value theorem, see *Preliminaries*, section 5, result 2. According to this, there is a $\\xi \\in (a,b)$ such that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\sum_{j=0}^{m-1} f^{(4)}(\\xi_j) = m f^{(4)}(\\xi).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Use $2mh = b-a$, and the following theorem has been proved:\n", "\n", "**Theorem: Error in composite Simpson's method.**\n", "\n", "Let $f(x) \\in C^{4}[a,b]$. There exist a $\\xi \\in (a,b)$ such that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\int_a^b f(x)dx - S_{m}(a,b) = -\\frac{(b-a)h^4}{180} f^{(4)}(\\xi).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Example 4:** \n", "Find the upper bound for the error when the composite Simpson's rule is applied\n", "to the integral $\\int_0^1 \\cos(\\pi x/2)dx$. \n", "\n", "In this case $f^{(4)}(x) = (\\pi^4/16) \\cos(\\pi\n", "x/2)$, so that $|f^{4)}(x)| \\leq (\\pi/2)^4$. The error bound becomes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "|I(a,b)-S_m(a,b)| \\leq \\frac{1}{180} \\left(\\frac{1}{2m}\\right)^4 \n", "\\left(\\frac{\\pi}{2}\\right)^4 = \\frac{\\pi^4}{46080}\\frac{1}{m^4}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If $m$ is increased by a factor 2, the error will be reduced by a factor of\n", "1/16, as indicated by Numerical experiment 2. \n", "\n", "**Numerical exercise:** \n", "Include the error bound in the output of Numerical experiment 2, and confirm that it really holds.\n", "\n", "\n", "## Error estimate\n", "From a practical point of view, the error expression derived above has some\n", "limitations. First, the bounds are often much too high. Second, we do not always\n", "know (or want to find) $|f^{(4)}(x)|$, even less its upper bound. So the question arises:\n", "How can we find an estimate of the error, without any extra analytical\n", "calculations? \n", "\n", "This is the idea: \n", "Let the interval $(a,b)$ chosen small, such that $f^{(4)}(x)$ can be\n", "assumed to be almost constant over the interval. Let $H=b-a$ be the length of the interval. Let $S_1(a,b)$ and $S_2(a,b)$ be the results from Simpson's formula over one and two subintervals respectively. Further, let $C = -f^{(4)}(x)/2880$ for some $x\\in [a,b]$, which $x$ does not matter since $f^{(4)}$ is assumed almost constant anyway. The errors of the two approximations are then given by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\\begin{align*}\n", "I(a,b) - S_1(a,b) &\\approx C H^5, \\\\ \n", "I(a,b) - S_2(a,b) &\\approx 2 C \\left(\\frac{H}{2}\\right)^5.\n", "\\end{align*}\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Subtract the two expressions to eliminate $I(a,b)$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "S_2(a,b) - S_1(a,b) \\approx \\frac{15}{16}C H^5\n", " \\qquad \\Rightarrow \\qquad\n", " CH^5 \\approx \\frac{16}{15}(S_2(a,b) - S_1(a,b)).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Insert this in the expression for the error:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "E_1(a,b) = I(a,b) - S_1(a,b) \\approx \\frac{16}{15} (\\,S_2(a,b) - S_1(a,b)\\, ) = \\mathcal{E}_1(a,b), \n", "\\label{_auto13} \\tag{13}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "E_2(a,b) = I(a,b) - S_2(a,b) \\approx \\frac{1}{15} (\\,S_2(a,b) - S_1(a,b)\\, ) = \\mathcal{E}_2(a,b).\n", "\\label{_auto14} \\tag{14}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This gives us a computable estimate for the error, both in $S_1$ and $S_2$. As\n", "the error in $S_2(a,b)$ is about 1/16 of the error in $S_1(a,b)$, and we anyway\n", "need to compute both, we will use $S_2(a,b)$ as our approximation. An even better\n", "approximation to the integral is given for free by just adding the error\n", "estimate.\n", "\n", "** Example 5:**\n", "Find an approximation to the integral $\\int_0^1\\cos(x)dx = \\sin(1)$ by Simpson's\n", "rule over one and two subintervals. Find the error estimates $\\mathcal{E}_m$,\n", "$m=1,2$ and compare with the exact error. \n", "\n", "*Solution:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "S_1(0,1) = \\frac{1}{6} \\big[ \\cos(0.0) + 4\\cos(0.5) + \\cos(1.0) \\big] = 0.8417720923 \n", "\\label{_auto15} \\tag{15}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "S_2(0,1) = \\frac{1}{12} \\big[ \\cos(0.0) + 4 \\cos(0.25) +2 \\cos(0.5) + 4 \\cos(0.75) + \\cos(1.0) \\big] = 0.8414893826\n", "\\label{_auto16} \\tag{16}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The exact error and the error estimate become:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation}\n", "E_1(0,1) = \\sin(1) - S_1(0,1) = -3.011 \\cdot 10^{-4}, \\quad\n", "\\mathcal{E}_1(0,1) = \\frac{16}{15}(S_2-S_1) = -3.016\\cdot 10^{-4}, \n", "\\label{_auto17} \\tag{17}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", "\n", "$$\n", "\\begin{equation} \n", "E_2(0,1) = \\sin(1)-S_2(0,1) = -1.840 \\cdot 10^{-5}, \\quad \n", "\\mathcal{E}_2(0,1) = \\frac{1}{16} (S_2-S_1) = -1.885 \\cdot 10^{-5}.\n", "\\label{_auto18} \\tag{18}\n", "\\end{equation}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case, it is a very good correspondence between the error estimate and\n", "the exact error. An even better approximation is obtained by adding the error\n", "estimate to $S_2$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "Q = S_{2}(0,1) + \\mathcal{E}_2(0,1) = 0.8414705353607151\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "with an error $\\sin(1)-Q = 4.4945 \\cdot 10^{-7}$. This gives a lot of additional\n", "accuracy without any extra work. \n", "\n", "### Implementation of Simpson's method with an error estimate\n", "\n", "The function simpson_basic returns" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "S_2(a,b) \\approx \\int_{a}^b f(x)dx\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "including an error estimate." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def simpson_basic(f, a, b):\n", " # Simpson's method with error estimate\n", " # Input: \n", " # f: integrand\n", " # a, b: integration interval\n", " # Output:\n", " # S_2(a,b) and the error estimate.\n", " \n", " # The nodes \n", " c = 0.5*(a+b)\n", " d = 0.5*(a+c)\n", " e = 0.5*(c+b)\n", " \n", " # Calculate S1=S_1(a,b), S2=S_2(a,b) \n", " H = b-a\n", " S1 = H*(f(a)+4*f(c)+f(b))/6\n", " S2 = 0.5*H*(f(a)+4*f(d)+2*f(c)+4*f(e)+f(b))/6\n", "\n", " error_estimate = (S2-S1)/15 # Error estimate for S2\n", " return S2, error_estimate" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Test:**\n", "As a first check of the implementation, use the example above, and make sure\n", "that the results are the same:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Test of simpson_basic\n", "\n", "def f(x): # Integrand\n", " return cos(x)\n", "\n", "a, b = 0, 1 # Integration interval\n", " \n", "I_exact = sin(1) # Exact solution for comparision\n", "\n", "# Simpson's method over two intervals, with error estimate\n", "S, error_estimate = simpson_basic(f, a, b)\n", "\n", "# Print the result and the exact solution \n", "print('Numerical solution = {:.8f}, exact solution = {:.8f}'.format(S, I_exact))\n", "\n", "# Compare the error and the error estimate \n", "print('Error in S2 = {:.3e}, error estimate for S2 = {:.3e}'.format(I_exact-S, error_estimate))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, let us see how reliable the quadrature and the error estimates are for\n", "another example, which you have to do yourself: \n", "\n", "**Numerical experiment 3:**\n", "Given the integral (with solution)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "I(a,b) = \\int_a^b \\frac{1}{1+16x^2} dx = \\left. \\frac{\\arctan(4x)}{4}\n", "\\right|_a^b\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Use simson_basic to find an approximation to the integral over the interval $[0,8]$. Print out $S_2(0,8)$, the error estimate $\\mathcal{E}_2(0,8)$ and the real error $E_2(0,8)$. How reliable are the error estimates?\n", "\n", "2. Repeat the experiment over the intervals $[0,1]$ and $[4, 8]$. Notice the difference between exact error of the two intervals.\n", "\n", "3. Repeat the experiment over the interval $[0,0.1]$.\n", "\n", "This is what you should observe from the experiment:\n", "1. Interval $[0,8]$: The error is large, and the error estimate is significantly smaller than the real error (the error is *under-estimated*).\n", "\n", "2. Interval $[0,1]$: As for the interval $[0,8]$. \n", "\n", "3. Interval $[4,8]$: Small error, and a reasonable error estimate.\n", "\n", "4. Interval $[0,0.1]$: Small error, reasonable error estimate.\n", "\n", "Why is it so, and how can we deal with it? Obviously, we need small subintervals\n", "near $x=0$, while large subintervals are acceptable in the last half of the\n", "interval. \n", "\n", "**Explanation:**\n", "The error in Simpson's method is given by" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "E(a,b) = -\\frac{(b-a)^5}{2880}f^{(4)}(\\xi).\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So let us take a look at $f^{(4)}(x)$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "f(x)=\\frac{1}{1+16x^2} \\qquad \\Rightarrow \\qquad\n", " f^{(4)}(x) = 6144 \\frac{1280 x^4 - 160x^2 +1}{(1-16x^2)^5}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot the 4th derivate of Runge's function:\n", "def df4(x):\n", " return 6144*(1280*x**4-160*x**2+1)/((1+16*x**2)**5)\n", "x = linspace(0, 8, 1001)\n", "plot(x, df4(x))\n", "title('The 4th derivative of Runges function');\n", "xlabel('x')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is no surprise that the error is large and the error estimates fail (we have assumed $f^{(4)}$ almost constant for the estimates) over the interval $[0,1]$. The part of the interval where $f^{(4)}(x)$ is large has to be partitioned in significantly smaller subintervals to get an acceptable result. But how, as $f^{(4)}$ is in general not known? This is the topic of the next section. \n", "\n", "## Adaptive integration\n", "\n", "Given a basic function, for example simpson_basic, returning an approximation $Q(a,b)$ to the integral, as well as an error estimate $\\mathcal{E}(a,b)$. Based on this, we want to find a partitioning of the interval:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "a = X_0 < X_1 \\cdots < X_m = b\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "such that" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "|\\mathcal{E}(X_j, X_{j+1})| \\approx \\frac{X_{k+1}-X_k}{b-a} \\cdot Tol\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $Tol$ is a tolerance given by the user. In this case" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\text{Accumulated error over (a,b)} \\approx \\sum_{j=0}^{m-1} \\mathcal{E}(X_k, X_{k+1})\n", " \\leq \\text{Tol}.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Such a partitioning can be done by an recursive algorithm:\n", "\n", "**Algorithm: Adaptive quadrature.**\n", "\n", "Given $f$, $a$, $b$ and a user defined tolerance Tol.\n", "* Calculate $Q(a,b)$ and $\\mathcal{E}(a,b)$.\n", "\n", "* **if** $|\\mathcal{E}(a,b)| \\leq \\text{Tol}$:\n", "\n", " * Accept the result, return $Q(a,b) + \\mathcal{E}(a,b)$ as an approximation to $I(a,b)$.\n", "\n", "\n", "* **else**:\n", "\n", " * Let $c=(a+b)/2$, and repeat the process on each of the subintervals $[a,c]$ and $[c,b]$, with tolerance $\\text{Tol}/2$.\n", "\n", "\n", "* Sum up the accepted results from each subinterval.\n", "\n", "\n", "\n", "### Implementation\n", "\n", "The adaptive algorithm is implemented below with simpson_basic as the basic quadrature routine. \n", "The function simpson_adaptive is a recursive function, that is a function that calls itself.\n", "To avoid it to do so infinitely many times, an extra variable level is introduced, this will increase by one for each time the function calls itself. If level is over some maximum value, the result is returned, and a warning printed." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def simpson_basic(f, a, b):\n", " # Simpson's method with error estimate\n", " # Input: \n", " # f: integrand\n", " # a, b: integration interval\n", " # Output:\n", " # S_2(a,b) and the error estimate.\n", " \n", " # The nodes \n", " c = 0.5*(a+b)\n", " d = 0.5*(a+c)\n", " e = 0.5*(c+b)\n", " \n", " # Calculate S1=S_1(a,b), S2=S_2(a,b) \n", " H = b-a\n", " S1 = H*(f(a)+4*f(c)+f(b))/6\n", " S2 = 0.5*H*(f(a)+4*f(d)+2*f(c)+4*f(e)+f(b))/6\n", "\n", " error_estimate = (S2-S1)/15 # Error estimate for S2\n", " return S2, error_estimate" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def simpson_adaptive(f, a, b, tol = 1.e-6, level = 0, maks_level=15):\n", " # Simpson's adaptive method\n", " # Input: \n", " # f: integrand\n", " # a, b: integration interval\n", " # tol: tolerance\n", " # level, maks_level: For the recursion. Just ignore them. \n", " # Output:\n", " # The approximation to the integral\n", " \n", " \n", " Q, error_estimate = simpson_basic(f, a, b) # The quadrature and the error estimate \n", " \n", " # -------------------------------------------------\n", " # Write the output, and plot the nodes. \n", " # This part is only for illustration. \n", " if level == 0:\n", " print(' l a b feil_est tol')\n", " print('==============================================') \n", " print('{:2d} {:.6f} {:.6f} {:.2e} {:.2e}'.format(\n", " level, a, b, abs(error_estimate), tol))\n", " \n", " x = linspace(a, b, 101)\n", " plot(x, f(x), [a, b], [f(a), f(b)], '.r')\n", " title('The integrand and the subintervals')\n", " # -------------------------------------------------\n", " \n", " if level >= maks_level:\n", " print('Warning: Maximum number of levels used.')\n", " return Q\n", " \n", " if abs(error_estimate) < tol: # Accept the result, and return\n", " result = Q + error_estimate \n", " else:\n", " # Divide the interval in two, and apply the algorithm to each interval.\n", " c = 0.5*(b+a)\n", " result_left = simpson_adaptive(f, a, c, tol = 0.5*tol, level = level+1)\n", " result_right = simpson_adaptive(f, c, b, tol = 0.5*tol, level = level+1)\n", " result = result_right + result_left\n", " return result" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Numerical experiment 4:**\n", "Use adaptive Simpson to find an approximation to the integral $\\int_0^5 1/(1+16x^2)dx$ using the tolerances Tol=$10^{-3}, 10^{-5}, 10^{-7}$. Compare the numerical result with the exact one." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Test: The adaptive Simpson's method\n", "def f(x): # Integrand \n", " return 1/(1+(4*x)**2)\n", "a, b = 0, 8 # Integration interval\n", "I_exact = 0.25*(arctan(4*b)-arctan(4*a)) # Exact integral\n", "tol = 1.e-3 # Tolerance\n", "# Apply the algorithm\n", "result = simpson_adaptive(f, a, b, tol=tol)\n", "# Print the result and the exact solution \n", "print('\\nNumerical solution = {:8f}, exact solution = {:8f}'\n", " .format(result, I_exact))\n", "# Compare the measured error and the tolerance\n", "err = I_exact - result\n", "print('\\nTolerance = {:.1e}, error = {:.3e}'.format(tol, abs(err)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Other quadrature formulas\n", "Simpson's rule is only one example of quadrature rules derived from polynomial interpolations. There are many others, and the whole process of deriving the methods, do error analysis, develope error estimates and adaptive algorithms can be repeated. \n", "\n", "Let us just conclude with a few other popular classes of methods: \n", "\n", "[Newton-Cotes formulas](https://en.wikipedia.org/wikTi/Newton–Cotes_formulas). These are based on equidistributed nodes. Both the Trapezoidal rule and Simpson's rule are examples of closed Newton-Cotes methods, the midpoint rule is an example of an open one. \n", "\n", "*Gauss-Legendre quadrature*: For the standard interval $[-1,1]$ choose the nodes as the zeros of the polynomial of degree $m$:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "L_m(t) = \\frac{d^m}{dt^m}(t^2-1)^m.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These are of precision $d=2m-1$, which is the best result possible with $m$ nodes. The midpoint rule is an example of a Gauss-Legendre quadrature with $m=1$." ] } ], "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.6.4" } }, "nbformat": 4, "nbformat_minor": 2 }