{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Definite Integrals, Part 3: The (Composite) Simpson's Rule and Richardson Extrapolation\n", "\n", "**References:**\n", "\n", "- Sections 5.2.2 and 5.2.3 of Chapter 5 *Numerical Differentiation and Integration* in {cite}`Sauer`.\n", "- Sections 4.3 and 4.4 of Chapter 5 *Numerical Differentiation and Integration* in {cite}`Burden-Faires`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "The Composite Simpson's Rule can be be derived in several ways.\n", "The traditional approach is to devise Simpson's Rule by approximating the integrand function with a colocating quadratic (using three equally spaced nodes) and then \"compounding\", as seen with the Trapezoid and Midpoint Rules.\n", "\n", "We have already seen another approach: using a 2:1 weighted average of the Trapezoid and Midpoint Rules with th goal of cancelling their $O(h^2)$ error terms.\n", "\n", "This section will show a third approach, based on Richardson extrapolation:\n", "this will set us up for Romberg Integration." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Basic Simpson's Rule by Richardson Extrapolation\n", "\n", "From the section on\n", "{doc}`The Composite Trapezoid and Midpoint Rules `,\n", "we have\n", "\n", "$$\n", "T_n = \\int_a^b f(x) \\,dx + \\frac{Df(b) - Df(a)}{12} h^2 + O(h^4),\n", "= I + c_2 h^2 + O(h^4)\n", "$$\n", "\n", "where $I$ is the integral to be approximated (the \"Q\" in the section on\n", "{doc}`richardson-extrapolation`,\n", "and $c_2 = (Df(b) - Df(a))/12$.\n", "\n", "Thus the \"n form\" of Richardson Extrapolation with $p=2$ gives a new approximation that I will call $S_{2n}$:\n", "\n", "$$ S_{2n} = \\frac{4T_{2n} - T_n}{4-1} $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To start, look at the simplest case of this:\n", "\n", "$$\n", "S_{2} = \\frac{4 T_{2} - T_1}{3}\n", "$$\n", "\n", "Definfing $h = (b-a)/2$, the ingredients are\n", "\n", "$$\n", "T_1 = \\frac{f(a) + f(b)}{2} (b-a) = \\frac{f(a) + f(b)}{2} 2 h = (f(a) + f(b)) h\n", "$$\n", "\n", "and\n", "\n", "$$\n", "T_2 = \\left[ \\frac{f(a)}{2} + f(a+h) + \\frac{f(b)}{2} \\right] h\n", "$$\n", "\n", "so\n", "\n", "$$\n", "S_{2} = \\frac{\\left[ 2f(a) + 4f(a+h) + 2f(b) \\right] - [f(a) + f(b)]}{3} h,\n", "= \\frac{f(a) + 4f(a+h) + f(b)}{3} h\n", "$$\n", "\n", "which is the basic Simpson's Rule.\n", "The subscript \"2\" is because this uses two intervals, with $h=(b-a)/2$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Accuracy and Order of Precision of Simpson's Rule\n", "\n", "Rather than derive this the traditional way — by fitting a quadratic to the function values at $x=a$, $a+h$ and $b$ —\n", "this can be confirmed \"a postiori\" by showing that the degree of precision is at least 2,\n", "so that it is exact for all quadratics. And actually we get a bonus, thanks to some symmetry." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**For $f(x) = 1$,** the exact integral is $I = b-a, = 2h$, and also\n", "\n", "$$ S_2 = \\frac{1 + 4 \\times 1 + 1}{3} h, = 2h $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**For $f(x) = x$,** the exact integral is\n", "$I = \\int_a^b x \\, dx = [x^2/2]_a^b = (b^2 - a^2)/2 = (b-a)(b+a)/2 = (a+b)h$\n", "\n", "and\n", "\n", "$$\n", "S_2 = \\frac{a + 4(a+b)/2 + b}{3} h\n", "= \\frac{a + 2(a+b) + b}{3} h\n", "= (a+b)h\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "However, it is sufficient to traslate the domain to the symmetric interval $[-h, h]$, so redo the $f(x)=x$ case this easier way:\n", "\n", "The exact integral is $\\int_{-h}^h x \\, dx = 0$ (because the function is odd)\n", "\n", "$$\n", "S_2 = \\frac{-h + 4 \\times 0 + h}{3} h = 0\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**For $f(x) = x^2$,** again do it just on the symmetric interval $[-h, h]$:\n", "the exact integral is\n", "$\\int_{-h}^h x^2 \\, dx = [x^3/3]_{-h}^h = 2h^3/3$\n", "and\n", "\n", "$$ S_2 = \\frac{(-h)^2 + 4\\times 0^2 + h^2}{3} h = 2h^3/3 $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So the degree of precision is *at least* 2, as expected." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What about cubics? Check with $f(x)=x^3$, again on interval $[-h, h]$.\n", "\n", "Almost no calculation is needed: symmetry does it all for us:\n", "\n", "- on one hand, the exact integral is zero due to the function being odd on a symmetric interval:\n", "$\\int_{-h}^h x^3 \\, dx = [x^4/4]_{-h}^h = 0$\n", "\n", "- on the other hand,\n", "\n", "$$ S_2 = \\frac{(-h)^3 + 4 \\times 0^3 + h^3}{3} h = 0 $$\n", "\n", "The degree of precision is at least 3.\n", "\n", "Our luck ends here, but looking at $f(x)=x^4$ is informative:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**For $f(x) = x^4$,**\n", "\n", "- the exact integral is $\\int_{-h}^h x^4 \\, dx = [x^5/5]_{-h}^h = 2h^5/5$\n", "\n", "- on the other hand\n", "\n", "$$ S_2 = \\frac{(-h)^4 + 4\\times 0^4 + h^4}{3} h = 2h^5/3 $$\n", "\n", "So there is a discrepancy of $(4/15) h^5, = O(h^5)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This Simpson's Rule has degree of precision 3: it is exact for all cubics, but not for all quartics.\n", "\n", "The last result also indicate the order of error:\n", "\n", "$$ S_2 - I = O(h^5) $$\n", "\n", "Just as for the composite Trapezoid and Midpoint Rules, when we combine multiple simple Simpson's Rule approximations with $2n$ intervals each of width $h= (b-a)/(2n)$, the error is roughly multiplied by $n$, so $h^5$ goes to $n h^5, = (b-a)h^4$,\n", "leading to\n", "\n", "$$ S_{2n} - I = O(h^4) $$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Appendix: Deriving Simpson's Rule by the Method of Undetermined Coefficients\n", "\n", "We wish the determine the most accurate approximation of the form\n", "\n", "$$ \\int_a^b f(x) \\, dx \\approx \\left[ C_1 f(a) + C_2 f(c) + C_3 f(b) \\right] h $$\n", "\n", "where $c$ is the midpoint, $c = (a+b)/2$\n", "\n", "This wilk be done by the first, \"hardest\" method: inserting Taylor polynomial and error terms,\n", "but to make it a bit less hard, we can consider just the symmetric case $a=-h$, $b=h$, $h= (b-a)/2$ by making the change of variables $x \\to x-c$.\n", "\n", "As we now know that this will be exact for cubics, use third order Tayloe polynomials:\n", "\n", "$$\n", "f(\\pm h) = f(0) \\pm f'(0) h + \\frac{f''(0)}{2} h^2 \\pm \\frac{f'''(0)}{6} h^3 + \\frac{f''''(\\xi_\\pm)}{24} h^4\n", "$$\n", "\n", "(Note that the special values $\\xi_\\pm$ are in general different for the \"$+h$\" and \"$-h$\" cases." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As usual, gather terms with the same power of $h$:\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{split}\n", "S_2 &= h f(0) (C_1 + C_2 + C_3)\n", "\\\\&+ h^2 f^{(1)}(0) (-C_1 + C_3)\n", "\\\\&+ h^3 f^{(2)}(0) (C_1/2 + C_3/2)\n", "\\\\&+ h^4 f^{(3)}(0) (-C_1/6 + C_3/6)\n", "\\\\&+ h^5 (C_1 f^{(4)}(\\xi_-) + C_3 f^{(4)}(\\xi_+))/24 \n", "\\end{split}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The exact integral can also be computed with Taylor's formula:\n", "\n", "$$\n", "\\begin{split}\n", "I = \\int_{-h}^h f(x)\\, dx\n", "&= \\int_{-h}^h \\left[ f(0) + Df(0) x + \\frac{D^2f(0)}{2}x^2 + \\frac{D^3f(0)}{6}x^3 + \\frac{D^4f(24)}{2}x^4 + \\frac{D^5f(\\xi_x)}{120}x^5 \\right]\\, dx\n", "\\\\&= 2 h f(0) + \\frac{D^2f(0)}{3} h^3 + \\frac{D^3f(0)}{12} h^5 + O(h^6)\n", "\\end{split}\n", "$$\n", "\n", "(Symmetry causes all the odd power integrals to valish.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "so the error is\n", "\n", "$$\n", "\\begin{split}\n", "S_2 - I\n", "&= h f(0) (C_1 + C_2 + C_3 - 2)\n", "\\\\&+ h^2 Df(0) (-C_1 + C_3)\n", "\\\\&+ h^3 D^2f(0) (C_1/2 + C_3/2 - 1/3)\n", "\\\\&+ O(h^5)\n", "\\end{split}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The best possibility is setting the coeficients of $h$, $h^2$ and $h^3$ to zero:\n", "\n", "$$\n", "\\begin{split}\n", "C_1 + C_2 + C_3 &= 2\n", "\\\\-C_1 + C_3 &= 0\n", "\\\\C_1/2 + C_3/2 &= 1/3\n", "\\end{split}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Symmetry helps, as the \"$h^2\"$ equation $-C_1 + C_3 = 0$ gives $C_3 = C_1$, leaving\n", "\n", "$$ C_1 = 1/3, \\quad 2C_1 + C_2 = 2 $$\n", "\n", "and thus\n", "\n", "$$ C_1 = C_3 = 1/3, \\quad C_2 = 4/3 $$\n", "\n", "as claimed above." ] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.8.1", "language": "julia", "name": "julia-1.8" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.8.1" } }, "nbformat": 4, "nbformat_minor": 4 }