{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# 08 Errors\n", "\n", "(See also *Computational Physics* (Landau, Páez, Bordeianu), Chapter 3)\n", "\n", "\n", "These slides include material from *Computational Physics. eTextBook Python 3rd Edition.* Copyright © 2012 Landau, Rubin, Páez. Used under the Creative-Commons Attribution-NonCommerical-ShareAlike 3.0 Unported License. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "## Stupidity or Incompetence\n", "(e.g., [PEBCAK](https://en.wiktionary.org/wiki/PEBCAK))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "## Random errors \n", "- cosmic rays\n", "- random bit flips" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Approximation errors\n", "\n", "\"**algorithmic errors**\"\n", " - simplifying and adapting mathematics to the computer\n", " - should decrease as $N$ increases\n", "\n", "#### Example:\n", "Approximate $\\sin(x)$ with its truncated series expansion:\n", "\\begin{align}\n", "\\sin x &= \\sum_{n=1}^{+\\infty} \\frac{(-1)^{n-1} x^{2n-1}}{(2n - 1)!}\\\\\n", " &\\approx \\sum_{n=1}^{N} \\frac{(-1)^{n-1} x^{2n-1}}{(2n - 1)!} + \\mathcal{E}(x, N)\n", "\\end{align}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "## Round-off errors\n", "- finite precision for storing floating-point numbers (32 bit, 64 bit)\n", "- not known exactly (treat as uncertainty)\n", "- can *accumulate* and lead to *garbage*\n", "\n", "#### Example: \n", "Assume you can only store four decimals:\n", "\n", "\\begin{align}\n", " \\text{storage}:&\\quad \\frac{1}{3} = 0.3333_c \\quad\\text{and}\\quad \\frac{2}{3} = 0.6667_c\\\\\n", " \\text{exact}:&\\quad 2\\times\\frac{1}{3} - \\frac{2}{3} = 0\\\\\n", " \\text{computer}:&\\quad 2 \\times 0.3333 - 0.6667 = -0.0001 \\neq 0\n", "\\end{align}" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "slideshow": { "slide_type": "fragment" } }, "source": [ "... now imagine adding \"$2\\times\\frac{1}{3} - \\frac{2}{3}$\" in a loop 100,000 times." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## The problems with *subtractive cancelation* \n", "Model the computer representation $x_c$ of a number $x$ as\n", "\n", "$$\n", "x_c \\simeq x(1+\\epsilon_x)\n", "$$\n", "\n", "with the *relative* error $|\\epsilon_x| \\approx \\epsilon_m$ (similar to machine precision).\n", "\n", "Note: The *absolute* error is $\\Delta x = x_c - x$ and is related to the relative error by $\\epsilon_x = \\Delta x/x$." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "What happens when we subtract two numbers $b$ and $c$: \n", "\n", "$$a = b - c$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\\begin{gather}\n", "a_c = b_c - c_c = b(1+\\epsilon_b) - c(1+\\epsilon_c)\\\\\n", "\\frac{a_c}{a} = 1 + \\frac{b}{a}\\epsilon_b - \\frac{c}{a} \\epsilon_c\n", "\\end{gather}\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "No guarantee that the errors cancel, and the relative error on $a$\n", "\n", "$$\n", "\\epsilon_a = \\frac{a_c}{a} - 1 = \\frac{b}{a}\\epsilon_b - \\frac{c}{a} \\epsilon_c\n", "$$ \n", "\n", "can be huge for small $a$!" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Subtracting two nearly equal numbers\n", "\n", "$$b \\approx c$$ is bad!" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "\\begin{align}\n", "\\frac{a_c}{a} &= 1 + \\frac{b}{a}(\\epsilon_b - \\epsilon_c) \\\\\n", "\\left| \\frac{a_c}{a} \\right| &\\leq 1 + \\left| \\frac{b}{a} \\right| (|\\epsilon_b| + |\\epsilon_a|)\n", "\\end{align}\n", "\n", "i.e. the large number $b/a$ magnifies the error. " ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "# Beware of subtractions!\n", "\n", "**If you subtract two large numbers and end up with a small one, then the small one is less significant than any of the large ones.**" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Round-off errors\n", "\n", "Repeated calculations of quantities with errors beget new errors: In general, analyze with the rules of *error propagation*: function $f(x_1, x_2, \\dots, x_N)$ with absolute errors on the $x_i$ of $\\Delta x_i$ (i.e., $x_i \\pm \\Delta x_i$):\n", "$$\n", "\\Delta f(x_1, x_2, \\dots; \\Delta x_1, \\Delta x_2, \\dots) =\n", " \\sqrt{\\sum_{i=1}^N \\left(\\Delta x_i \\frac{\\partial f}{\\partial x_i}\\right)^2}\n", "$$\n", "\n", "Note: relative error $$\\epsilon_i = \\frac{\\Delta x_i}{x_i}$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Example: division $a = b/c$ (... with short cut)\n", "\\begin{align}\n", "a_c &= \\frac{b_c}{c_c} = \\frac{b(1+\\epsilon_b)}{c(1+\\epsilon_b)} \\\\\n", "\\frac{a_c}{a} &= \\frac{1+\\epsilon_b}{1+\\epsilon_c} \n", " = \\frac{(1+\\epsilon_b)(1-\\epsilon_c)}{1-\\epsilon_c^2} \\approx (1+\\epsilon_b)(1-\\epsilon_c)\\\\\n", " &\\approx 1 + |\\epsilon_b| + |\\epsilon_c|\\\\\n", "\\epsilon_a = \\frac{a_c}{a} - 1 &\\approx |\\epsilon_b| + |\\epsilon_c|\n", "\\end{align}\n", "\n", "(neglected terms of order $\\mathcal{O}(\\epsilon^2)$); and same for multiplication.\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "**Errors accumulate with every operation.**" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Model for round-off error accumulation\n", "View error in each calculation as a step in a *random walk*. The total \"distance\" (i.e. total error) $R$ over $N$ steps of length $r$ (the individual, \"random\" errors), is on average" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "$$ R \\approx \\sqrt{N} r $$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Total relative error $\\epsilon_{\\text{ro}}$ after $N$ calculations with error of the order of the machine precision $\\epsilon_m$ is\n", "\n", "$$ \\epsilon_{\\text{ro}} \\approx \\sqrt{N} \\epsilon_m $$\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "(Only a model, depending on algorithm may be less or even $N!$...)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "## Total error of an algorithm\n", "What you need to know to evaluate an algorithm:\n", "1. Does it converge? (What $N$ do I need?)\n", "2. How precise are the converged results (What is the error $\\epsilon_\\text{tot}$?)\n", "3. What is its run time? (How fast is it for a given problem size?)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "The total error contains *approximation* and *round off* errors:\n", "\n", "\\begin{gather}\n", "\\epsilon_\\text{tot} = \\epsilon_\\text{app} + \\epsilon_\\text{ro}\n", "\\end{gather}\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Model for the approximation error for an algorithm that takes $N$ steps (operations) to find a \"good\" answer:\n", "\n", "$$\n", "\\epsilon_\\text{app} \\simeq \\frac{\\alpha}{N^\\beta}\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "and round off error as\n", "\n", "$$\n", "\\epsilon_{\\text{ro}} \\approx \\sqrt{N} \\epsilon_m\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Model for total error:\n", "$$\n", "\\epsilon_\\text{tot} = \\frac{\\alpha}{N^\\beta} + \\sqrt{N} \\epsilon_m\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Analyze $\\log_{10} $ of the relative error (direct readout of number of significant decimals).\n", "\n", "\n", "\n", "Image from Computational Physics. eTextBook Python 3rd Edition. Copyright © 2012 Landau, Rubin, Páez. Used under the Creative-Commons Attribution-NonCommerical-ShareAlike 3.0 Unported License." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "### Example analysis\n", "\\begin{gather}\n", "\\epsilon_\\text{app} = \\frac{1}{N^2}, \\quad \\epsilon_\\text{ro} = \\sqrt{N}\\epsilon_m\\\\\n", "\\epsilon_\\text{tot} = \\frac{1}{N^2} + \\sqrt{N}\\epsilon_m\n", "\\end{gather}" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Total error is a *minimum* for\n", "\n", "\\begin{gather}\n", "\\frac{d\\epsilon_\\text{tot}}{dN} = -\\frac{2}{N^{3}} + \\frac{1}{2}\\frac{\\epsilon_m}{\\sqrt{N}} = 0, \\quad\\text{thus} \\quad\n", "N^{5/2} = 4 \\epsilon_m^{-1}\\\\\n", "N = \\left(\\frac{4}{\\epsilon_m}\\right)^{2/5}\n", "\\end{gather}\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "What is the best $N$ for single precision $\\epsilon_m \\approx 10^{-7}$?" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [], "source": [ "import math\n", "def N_opt(eps_m):\n", " return round(math.pow(4./eps_m, 2./5.))\n", "def eps_app(N):\n", " return 1./(N*N)\n", "def eps_ro(N, eps_m):\n", " return math.sqrt(N)*eps_m" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "slideshow": { "slide_type": "fragment" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "best N = 1099 (for eps_m=1e-07)\n", "eps_tot = 4.14e-06\n", "eps_app = 8.28e-07, eps_ro = 3.32e-06\n" ] } ], "source": [ "epsilon_m = 1e-7 # single precision\n", "\n", "N = N_opt(epsilon_m)\n", "err_app = eps_app(N)\n", "err_ro = eps_ro(N, epsilon_m)\n", "print(\"best N = {0} (for eps_m={1})\".format(N, epsilon_m))\n", "print(\"eps_tot = {0:.3g}\".format(err_app + err_ro))\n", "print(\"eps_app = {0:.3g}, eps_ro = {1:.3g}\".format(err_app, err_ro))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "Single precision $\\epsilon_m \\approx 10^{-7}$:\n", "\n", "$$\n", "N \\approx 1099\\\\\n", "\\epsilon_\\text{tot} \\approx 4 \\times 10^{-6} \\\\\n", "\\epsilon_\\text{app} = 8.28 \\times 10^{-7} \\\\\n", "\\epsilon_\\text{ro} = 3.32 \\times 10^{-6}\n", "$$\n" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Here, most of the error is round-off error! What can you do?" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "* use double precision (delay round-off error)\n", "* use a better algorithm, e.g. $\\epsilon_\\text{app}\\simeq \\frac{2}{N^4}$ (uses fewer steps)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "**Better algorithms are always a good idea :-)**\n", "\n", "Remember: trade-off between **approximation error** and **rounding error**." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "slideshow": { "slide_type": "skip" } }, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "celltoolbar": "Slideshow", "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.7" } }, "nbformat": 4, "nbformat_minor": 2 }