{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 02 - Error, accuracy, stability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Floating point representation of numbers, roundoff error, truncation error, numerical stability and condition number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Main sources of numerical errors: \n",
    "- rounding and cancellation (due to the use of finite-precision arithmetic)\n",
    "- truncation or approximation errors (due to approximation of infinite sequences or continuous functions by a finite number of samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic definitions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Absolute and relative error:**\n",
    "\n",
    "If $ \\tilde{x} $ is an approximation for $ x $, then:\n",
    "\n",
    "  * absolute error $ A(x) = | \\tilde{x} - x | $\n",
    "\n",
    "  * relative error $ R(x) = | \\tilde{x} - x | \\, / \\, | x | $\n",
    "  \n",
    "**Decimal precision:**\n",
    "\n",
    "Given a relative error $ R $, the decimal precision $ p $ is the largest integer such that $ R \\leq 5 \\times 10^{-p} $\n",
    "\n",
    "**Big-O notation:**\n",
    "\n",
    "The error term in an approximation to a mathematical function can be described by the [big-O](https://en.wikipedia.org/wiki/Big_O_notation) notation:\n",
    "\n",
    "$$\n",
    "f(x) =  O(g(x)) \\quad \\text{as} \\quad x \\rightarrow a\n",
    "$$ \n",
    "\n",
    "if and only if \n",
    "\n",
    "$$\n",
    "|f(x)| \\leq M |g(x)| \\quad \\text{as}\\quad  |x - a| < \\delta \\quad \\text{where} \\quad M, a > 0.\n",
    "$$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Roundoff error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Representation of real numbers in computer\n",
    "\n",
    "Most widely used representation of real numbers in computers is the floating-point representation. Floating-point representations have a base $ \\beta $, exponent $ E $ and precision $ p $. In general, a floating-point number is represented as \n",
    "\n",
    "$$ f = \\pm \\, d_1.d_2d_3 \\dots d_p \\times \\beta^E, $$\n",
    "\n",
    "where $ d_1.d_2d_3 \\dots d_p $ is called significand (also mantissa)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Properties of floating-point systems:**\n",
    "\n",
    "- Smallest positive number ([underflow](https://en.wikipedia.org/wiki/Arithmetic_underflow) if below)\n",
    "- Largest number ([overflow](https://en.wikipedia.org/wiki/Integer_overflow) if above)\n",
    "- [Machine epsilon](https://en.wikipedia.org/wiki/Machine_epsilon), $ \\varepsilon $, defined as the difference between 1 and the next larger floating point number (upper bound on the relative error due to rounding in floating point arithmetic)\n",
    "- Special values: zero (`0`), infinities (`+Inf`, `-Inf`), [not a number](https://en.wikipedia.org/wiki/NaN) (`NaN`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 02.1:** Consider a decimal system of floating-point numbers defined as $ f = \\pm \\, d_1.d_2 \\cdot \\beta^E, $ where $ \\beta = 10 $ and $ E \\in \\{-2, -1, 0 \\} $.  \n",
    "\n",
    "1) How many numbers such system represents?  \n",
    "2) What is the smallest positive number and the largest number?  \n",
    "3) What is the machine epsilon?  \n",
    "4) What is the distribution of numbers on the real line?\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 02.2:** Consider a binary system of floating-point numbers defined as $ f = \\pm \\, d_1.d_2 \\cdot \\beta^E, $ where $ \\beta = 2 $ and $ E = \\{-1, 0, 1 \\} $.\n",
    "    \n",
    "1) How many numbers such system represents?  \n",
    "2) What is the smallest positive number and the largest number?  \n",
    "3) What is the machine epsilon?  \n",
    "4) What is the distribution of numbers on the real line?\n",
    "    \n",
    "</div>    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that unlike the real number system (which is continuous), a floating-point systems have always gaps between numbers. If a number is not exactly representable, then it must be approximated by one of the nearest representable values. Furthermore, the distribution of numbers on the real line is not uniform."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation of floating-point systems\n",
    "\n",
    "Over the years, a variety of floating-point representations have been used in computers. In 1985, the [IEEE 754](https://en.wikipedia.org/wiki/IEEE_754) standard for floating-point arithmetic was established. Since then, most implementations of the floating-point systems in computers conform to the rules defined by IEEE."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Half precision** ([IEEE 754-2008](https://en.wikipedia.org/wiki/IEEE_754-2008_revision)): \n",
    "\n",
    "Total storage alloted is 16 bits (10 bits for mantissa, 5 bits for exponent, 1 bit for sign)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.finfo(np.float16))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Single precision:**\n",
    "\n",
    "Total storage alloted is 32 bits (23 bits for mantissa, 8 bits for exponent, 1 bit for sign)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.finfo(np.float32))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Double precision:**\n",
    "\n",
    "Total storage allocated is 64 bits (52 bits for mantissa, 11 bits for exponent, 1 bit for sign)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print(np.finfo(np.float64))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The IEEE 754 standard alse defines algorithm for addition, subtraction, multiplication, division and square root, and requires that implementations produce the same result as that algorithm. Thus, when a program is moved from one machine to another, the results of the basic operations will be the same in every bit if both machines support the IEEE 754 standard. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Example:** Using the double precision floating-point format, compare that $ 0.1 + 0.2 = 0.3 $.\n",
    "    \n",
    "</div>  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print((0.1 + 0.2) == 0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result above is not a Python bug (read more e.g. [here](https://docs.python.org/3/faq/design.html#why-are-floating-point-calculations-so-inaccurate) or [here](https://docs.python.org/3/tutorial/floatingpoint.html)). Representing infinitely many real numbers by a finite number of bits requires an approximate representation. Given any fixed number of bits, most calculations with real numbers will produce quantities that cannot be exactly represented using that many bits. None of the numbers 0.1, 0.2, and 0.3 (actually most of the decimal fractions) has an exact representation as a binary floating-point number, no matter how many base 2 digits you’re willing to use. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "**Note:** The decimal numbers could be represented exactly by the base-10 floating-point systems. However, base-10 implementations are rare because base-2 (binary) arithmetic is so much faster on computers.\n",
    "    \n",
    "</div>    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Python by default displays always a rounded value (to keep the number of digits manageable):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "but the actual stored value is the nearest representable binary fraction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(format(0.1, \".55f\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 02.3:** Using the double precision floating-point format, calculate $ x = 0.1 + 0.2 - 0.3 $. Afterwards, perform $ 100 \\ \\times $ the operation $ x := x + x $ and display the value of $ x $.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 02.4:** Consider the following function, \n",
    "    \n",
    "$$\n",
    "f(x) = \\frac{1 - \\cos{x}}{x^2}.\n",
    "$$\n",
    "\n",
    "The behaviour of this function as $ x $ approaches zero can be determined by evaluating the limit \n",
    "    \n",
    "$$ \n",
    "\\lim_{x \\to 0} f(x) = \\frac{1}{2}. \n",
    "$$\n",
    "\n",
    "Can be found the same result using the floating-point representation?\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 02.5:** Calculate the sum of the following series,\n",
    "    \n",
    "$$\n",
    "0.9^0, \\ 0.9^1, \\ \\dots, \\ 0.9^n,\n",
    "$$\n",
    "    \n",
    "for $ n = 400 $ in the forward and backward directions. Compare the results.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "**Note:** When summing a series of numbers using floating-point systems, always sum the smaller numbers first. To improve accuracy when summing numbers, consider the [Kahan summation algorithm](https://en.wikipedia.org/wiki/Kahan_summation_algorithm).\n",
    "    \n",
    "</div>    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary\n",
    "\n",
    "- Floating point arithmetic is not commutative, associative, and not necessarily distributive:\n",
    "\n",
    "valid operations:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "& 1 \\cdot x = x, \\\\\n",
    "& x \\cdot y = y \\cdot x, \\\\\n",
    "& x + x = 2 \\cdot x \n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "not necessarily valid operations:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "& x \\cdot x^{-1} = 1, \\\\\n",
    "& (1 + x) - 1 = x, \\\\\n",
    "& (x + y) + z = x + (y + z) \n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "- Errors propagate between calculations. A small error in the input may result in a large error in the output.\n",
    "\n",
    "**Which operations should I avoid in order to minimize the roundoff error?**\n",
    "\n",
    "- subtractions of numbers that are nearly equal\n",
    "- additions and subtractions of numbers that differ greatly in magnitude\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Representation of integers in computer\n",
    "\n",
    "Although standard arithmetic operations are safe, some care may be necessary when working with large numbers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**16 bit unsigned**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.iinfo(np.uint16))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**16 bit signed**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.iinfo(np.int16))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**32 bit signed**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.iinfo(np.int32))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Truncation error\n",
    "\n",
    "The discrepancy between the true answer and the answer obtained by a numerical method regardless the roundoff error. Truncation error would persist even if we would have an infinitely accurate representation of numbers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 02.6:** Construct a Taylor series of the function \n",
    "\n",
    "$$\n",
    "f(x) = e^x\n",
    "$$\n",
    "\n",
    "using the first 3 elements. \n",
    "    \n",
    "1) Plot the Taylor series in $ x \\in [-1, 1] $ and compare the approximation with $ f(x) $.   \n",
    "2) Calculate and plot the absolute and relative errors of the approximation. Find the maximum values.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combination of errors\n",
    "\n",
    "In general, the roundoff errors and the truncation errors may combine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 02.7:** Calculate the derivative of the function\n",
    "\n",
    "$$\n",
    "f(x) = \\sin{x}\n",
    "$$\n",
    "\n",
    "in $ x \\in [-2 \\pi, 2 \\pi] $ using finite differences with [forward and central](https://en.wikipedia.org/wiki/Finite_difference#Basic_types) schemes, respectively, i.e.,\n",
    "\n",
    "$$\n",
    "f_{f}^{\\prime}(x) = \\frac{f(x + h) - f(x)}{h}, \\qquad f_{c}^{\\prime}(x) = \\frac{f(x + h) - f(x - h)}{2 h},\n",
    "$$\n",
    "\n",
    "where $ h $ is a finite step. \n",
    "    \n",
    "1) Compare the approximations with $ f^{\\prime}(x) = \\cos{x}.$\n",
    "    \n",
    "2) Calculate the relative error of the approximations above for $ x = 1 $ and $ h = 2^{-n} $, where $ n \\in \\{ 1, \\dots, 60 \\} $. Find the optimal value of $ h $ for each approximation.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The relative error drops faster in the case of central scheme, because it is second order of accuracy, i.e.,\n",
    "\n",
    "$$\n",
    "f(x + h) = f(x) + h f^{\\prime}(x) + \\frac{h^2}{2} f^{\\prime\\prime}(x) + \\mathcal{O}(h^3),\n",
    "$$\n",
    "$$\n",
    "f(x - h) = f(x) - h f^{\\prime}(x) + \\frac{h^2}{2} f^{\\prime\\prime}(x) + \\mathcal{O}(h^3),\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\to f^{\\prime}(x) = \\frac{f(x + h) - f(x - h)}{2 h} + \\mathcal{O}(h^2),\n",
    "$$\n",
    "\n",
    "while the forward scheme is first order of accuracy, i.e.,\n",
    "\n",
    "$$\n",
    "f(x + h) = f(x) + h f^{\\prime}(x) + \\mathcal{O}(h^2) \\ \\to \\ f^{\\prime}(x) = \\frac{f(x + h) - f(x)}{h} + \\mathcal{O}(h).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical stability\n",
    "\n",
    "Numerical method may magnify the roundoff error introduced into the computation at an early stage. Such a method is called unstable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 02.8:** Solve the 1st order ordinary differential equation\n",
    "\n",
    "$$\n",
    "y^{\\prime}(x) + y(x) = 0, \\quad y(0) = 1\n",
    "$$\n",
    "\n",
    "in $ x \\in [0, 10] $ using forward and central differencing schemes with $ h = 0.1 $. Compare both results with analytical solution $ y(x) = \\exp(-x) $.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "**Note:** Even though the central differencing scheme has higher order of accuracy than the forward scheme, it it not always better. The suitability has to be examined for each problem separately.\n",
    "    \n",
    "</div>    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 02.9:** Using the half-precision floating-point format, calculate the first 20 integer powers of the number called \"golden mean\",\n",
    "\n",
    "$$\n",
    "\\phi = \\frac{\\sqrt{5} - 1}{2}.\n",
    "$$\n",
    "    \n",
    "with the following recursive algorithm,\n",
    "    \n",
    "$$\n",
    "\\phi^n = \\phi^{n-1} - \\phi^{n}.\n",
    "$$\n",
    "\n",
    "Compare the calculated values with the ones obtained by the simple multiplication.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Condition number\n",
    "\n",
    "[Condition number](https://en.wikipedia.org/wiki/Condition_number) measures how much the output value of the function can change for a small change in the input argument (i.e., sensitivity to input error). Given the problem $ f $ and the input $ x $, the condition number $ C_p $ is defined as \n",
    "\n",
    "$$\n",
    "C_p = \\frac{\\| \\delta f(x) \\| \\, / \\, \\| f(x) \\|}{\\| \\delta x \\| \\, / \\, \\| x \\|}.\n",
    "$$\n",
    "\n",
    "- A problem with $ C_p \\approx 1 $ is said to be well-conditioned\n",
    "- A problem with $ C_p > 100 $ is said to be ill-conditioned\n",
    "- A problem with $ C_p > \\varepsilon^{-1} $ is not solvable within the specified precision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Exercise 02.10:** Find the condition number of the following system of equations,\n",
    "\n",
    "$$ \n",
    "x + \\alpha y = 1, \n",
    "$$\n",
    "\n",
    "$$ \n",
    "\\alpha x + y = 0, \n",
    "$$\n",
    "\n",
    "depending on the parameter $ \\alpha \\in [-2, 2] $. Plot the result and decide for which $ \\alpha $ the system is ill-conditioned.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "# add your code here\n"
   ]
  }
 ],
 "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": 2
}