Cap01_Error.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\">Métodos Numéricos</h1>\n",
    "<h1 align=\"center\">Capítulo 1</h1>\n",
    "<h1 align=\"center\">Error y Representación de números en el computador</h1>\n",
    "<h1 align=\"center\">2023/02</h1>\n",
    "<h1 align=\"center\">MEDELLÍN - COLOMBIA </h1>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*** \n",
    "||![Gmail](https://img.shields.io/badge/Gmail-D14836?style=plastic&logo=gmail&logoColor=white)| <carlosalvarezh@gmail.com>| |\n",
    "|-:|:-|--:|:--|\n",
    "|[![LinkedIn](https://img.shields.io/badge/linkedin-%230077B5.svg?style=plastic&logo=linkedin&logoColor=white)](https://www.linkedin.com/in/carlosalvarez5/)|[![@alvarezhenao](https://img.shields.io/twitter/url/https/twitter.com/alvarezhenao.svg?style=social&label=Follow%20%40alvarezhenao)](https://twitter.com/alvarezhenao)|[![@carlosalvarezh](https://img.shields.io/badge/github-%23121011.svg?style=plastic&logo=github&logoColor=white)](https://github.com/carlosalvarezh)|[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/carlosalvarezh/Metodos_Numericos/blob/master/Cap01_Error.ipynb)|\n",
    "\n",
    "<table>\n",
    " <tr align=left><td><img align=left src=\"https://github.com/carlosalvarezh/Curso_CEC_EAFIT/blob/main/images/CCLogoColorPop1.gif?raw=true\" width=\"25\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license.(c) Carlos Alberto Alvarez Henao</td>\n",
    "</table>\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Comentario:*** este capítulo está basado en parte de las  notas del curso del profesor [Kyle T. Mandli](https://github.com/mandli/intro-numerical-methods) (en inglés)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C01_Img01.png?raw=true\" width=\"500\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://www.bugsnag.com/blog/bug-day-ariane-5-disaster\">Ariane 5 Explosion</a> </div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Bibliotecas a ser utilizadas en el Notebook\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sympy\n",
    "import scipy.special\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fuentes de error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Los cálculos numéricos, que involucran el uso de máquinas (análogas o digitales) presentan una serie de errores que provienen de diferentes fuentes:\n",
    "\n",
    "- del Modelo\n",
    "\n",
    "- de los datos\n",
    "\n",
    "- de truncamiento\n",
    "\n",
    "- de representación de los números (punto flotante)\n",
    "\n",
    "- $ \\ldots$\n",
    "\n",
    "***Meta:*** Categorizar y entender cada tipo de error y explorar algunas aproximaciones simples para analizarlas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Error en el modelo y los datos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Errores en la formulación fundamental\n",
    "\n",
    "- Error en los datos: imprecisiones en las mediciones o incertezas en los parámetros\n",
    "\n",
    "Infortunadamente no tenemos control de los errores en los datos y el modelo de forma directa pero podemos usar métodos que pueden ser más robustos en la presencia de estos tipos de errores."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Error de truncamiento"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Los errores surgen de la expansión de funciones con una función simple, por ejemplo, $sin(x) \\approx x$ para $|x|\\approx0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Error de  representación de punto fotante"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Los errores surgen de aproximar números reales con la representación  en precisión finita de números en el computador."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definiciones básicas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dado un valor verdadero de una función $f$ y una solución aproximada $F$, se define:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Error absoluto"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$e_a=|f-F|$$\n",
    "\n",
    "\n",
    "***Ejemplo:*** se realiza una medición y se obtiene un valor aproximado de $29.99$ *m*. Asumiendo que el valor exacto de dicha medición debería ser de $30.00$ *m*, ¿cuál es el error absoluto obtenido? Cuál sería el error absoluto si se disminuye un orden de magnitud las cantidades?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = 30.0\n",
    "F = 29.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ea = abs(f - F)\n",
    "print(\"{0:6.4f}\".format(ea)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reduciendo un órden de magnitud las cantidades\n",
    "\n",
    "f = 3.0\n",
    "F = 2.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ea = abs(f - F)\n",
    "print(\"{0:6.4f}\".format(ea)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se observa que el valor del error absoluto es igual ($\\approx 0.1$), independiente de la magnitud de las cantidades. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Error relativo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$e_r (\\%)= \\frac{e_a}{|f|}=\\frac{|f-F|}{|f|} \\times 100 \\%$$\n",
    "\n",
    "\n",
    "***Ejemplo:*** Repetir el ejemplo anterior, pero calculando el error relativo porcentual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = 30.0\n",
    "F = 29.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "er = abs(f - F) / f * 100\n",
    "print(\"{0:6.4f}%\".format(er))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reduciendo un órden de magnitud las cantidades\n",
    "\n",
    "f = 3.0\n",
    "F = 2.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "er = abs(f - F) / f * 100\n",
    "print(\"{0:6.4f}%\".format(er))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se observa que los resultados son diferentes y es mayor cuando las cantidades medidas son menores.\n",
    "\n",
    "Entre las dos formas de representar el error, la relativa es más consistente con la magnitud de lo que se está midiendo."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notación $\\text{Big}-\\mathcal{O}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "sea $$f(x)= \\mathcal{O}(g(x)) \\text{ cuando } x \\rightarrow a$$\n",
    "\n",
    "si y solo si\n",
    "\n",
    "$$|f(x)|\\leq M|g(x)| \\text{ cuando } |x-a| < \\delta \\text{ donde } M, a > 0$$\n",
    "\n",
    "\n",
    "En la práctica, usamos la notación $\\text{Big}-\\mathcal{O}$ para decir algo sobre cómo se pueden comportar los términos que podemos haber dejado fuera de una serie. Veamos el siguiente ejemplo de la aproximación de la serie de Taylor:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Ejemplo:***\n",
    "\n",
    "sea $f(x) = \\sin(x)$ con $x_0 = 0$ entonces\n",
    "\n",
    "$$T_N(x) = \\sum^N_{n=0} (-1)^{n} \\frac{x^{2n+1}}{(2n+1)!}$$\n",
    "\n",
    "Podemos escribir $f(x)$ como\n",
    "\n",
    "$$f(x) = x - \\frac{x^3}{6} + \\frac{x^5}{120} + \\mathcal{O}(x^7)$$\n",
    "\n",
    "Esto se vuelve más útil cuando lo vemos como lo hicimos antes con $\\Delta x$:\n",
    "\n",
    "$$f(x) = \\Delta x - \\frac{\\Delta x^3}{6} + \\frac{\\Delta x^5}{120} + \\mathcal{O}(\\Delta x^7)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reglas para el error de propagación basado en la notación $\\text{Big}-\\mathcal{O}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En general, existen dos teoremas que no necesitan prueba y se mantienen cuando el valor de $x$ es grande:\n",
    "\n",
    "Sea\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    f(x) &= p(x) + \\mathcal{O}(x^n) \\\\\n",
    "    g(x) &= q(x) + \\mathcal{O}(x^m) \\\\\n",
    "    k &= \\max(n, m)\n",
    "\\end{aligned}$$\n",
    "\n",
    "Entonces\n",
    "\n",
    "$$\n",
    "    f+g = p + q + \\mathcal{O}(x^k)\n",
    "$$\n",
    "\n",
    "y\n",
    "\n",
    "\\begin{align}\n",
    "    f \\cdot g &= p \\cdot q + p \\mathcal{O}(x^m) + q \\mathcal{O}(x^n) + O(x^{n + m}) \\\\\n",
    "    &= p \\cdot q + \\mathcal{O}(x^{n+m})\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "De otra forma, si estamos interesados en valores pequeños de $x$, $\\Delta x$, la expresión puede ser modificada como sigue:\n",
    "\n",
    "\\begin{align}\n",
    "    f(\\Delta x) &= p(\\Delta x) + \\mathcal{O}(\\Delta x^n) \\\\\n",
    "    g(\\Delta x) &= q(\\Delta x) + \\mathcal{O}(\\Delta x^m) \\\\\n",
    "    r &= \\min(n, m)\n",
    "\\end{align}\n",
    "\n",
    "entonces\n",
    "\n",
    "$$\n",
    "    f+g = p + q + O(\\Delta x^r)\n",
    "$$\n",
    "\n",
    "y\n",
    "\n",
    "\\begin{align}\n",
    "    f \\cdot g &= p \\cdot q + p \\cdot \\mathcal{O}(\\Delta x^m) + q \\cdot \\mathcal{O}(\\Delta x^n) + \\mathcal{O}(\\Delta x^{n+m}) \\\\\n",
    "    &= p \\cdot q + \\mathcal{O}(\\Delta x^r)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Nota:*** En este caso, supongamos que al menos el polinomio con $k=max(n,m)$ tiene la siguiente forma:\n",
    "\n",
    "$$\n",
    "    p(\\Delta x) = 1 + p_1 \\Delta x + p_2 \\Delta x^2 + \\ldots\n",
    "$$\n",
    "\n",
    "o\n",
    "\n",
    "$$\n",
    "    q(\\Delta x) = 1 + q_1 \\Delta x + q_2 \\Delta x^2 + \\ldots\n",
    "$$\n",
    "\n",
    "para que $\\mathcal{O}(1)$ \n",
    "\n",
    "\n",
    "de modo que hay un término $\\mathcal{O}(1)$ que garantiza la existencia de $\\mathcal{O}(\\Delta x^r)$ en el producto final."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para tener una idea de por qué importa más la potencia en $\\Delta x$ al considerar la convergencia, la siguiente figura muestra cómo las diferentes potencias en la tasa de convergencia pueden afectar la rapidez con la que converge nuestra solución. Tenga en cuenta que aquí estamos dibujando los mismos datos de dos maneras diferentes. Graficar el error como una función de $\\Delta x$ es una forma común de mostrar que un método numérico está haciendo lo que esperamos y muestra el comportamiento de convergencia correcto. Dado que los errores pueden reducirse rápidamente, es muy común trazar este tipo de gráficos en una escala log-log para visualizar fácilmente los resultados. Tenga en cuenta que si un método fuera realmente del orden $n$, será una función lineal en el espacio log-log con pendiente $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = np.linspace(1.0, 1e-4, 100)\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2.0)\n",
    "axes = []\n",
    "axes.append(fig.add_subplot(1, 2, 1))\n",
    "axes.append(fig.add_subplot(1, 2, 2))\n",
    "\n",
    "for n in range(1, 5):\n",
    "    axes[0].plot(dx, dx**n, label=\"$\\Delta x^%s$\" % n)\n",
    "    axes[1].loglog(dx, dx**n, label=\"$\\Delta x^%s$\" % n)\n",
    "\n",
    "axes[0].legend(loc=2)\n",
    "axes[1].set_xticks([10.0**(-n) for n in range(5)])\n",
    "axes[1].set_yticks([10.0**(-n) for n in range(16)])\n",
    "axes[1].legend(loc=4)\n",
    "for n in range(2):\n",
    "    axes[n].set_title(\"Crecimiento del Error vs. $\\Delta x^n$\")\n",
    "    axes[n].set_xlabel(\"$\\Delta x$\")\n",
    "    axes[n].set_ylabel(\"Error Estimado\")\n",
    "    axes[n].set_title(\"Crecimiento de las diferencias\")\n",
    "    axes[n].set_xlabel(\"$\\Delta x$\")\n",
    "    axes[n].set_ylabel(\"Error Estimado\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error de truncamiento"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Teorema de Taylor:*** Sea $f(x) \\in C^{m+1}[a,b]$ y $x_0 \\in [a,b]$, para todo $x \\in (a,b)$ existe un número $c = c(x)$ que se encuentra entre $x_0$ y $x$ tal que\n",
    "\n",
    "$$ f(x) = T_N(x) + R_N(x)$$\n",
    "\n",
    "donde $T_N(x)$ es la aproximación del polinomio de Taylor\n",
    "\n",
    "$$T_N(x) = \\sum^N_{n=0} \\frac{f^{(n)}(x_0)\\times(x-x_0)^n}{n!}$$\n",
    "\n",
    "y $R_N(x)$ es el residuo (la parte de la serie que obviamos)\n",
    "\n",
    "$$R_N(x) = \\frac{f^{(n+1)}(c) \\times (x - x_0)^{n+1}}{(n+1)!}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Otra forma de pensar acerca de estos resultados consiste en reemplazar $x - x_0$ con $\\Delta x$. La idea principal es que el residuo $R_N(x)$ se vuelve mas pequeño cuando $\\Delta x \\rightarrow 0$.\n",
    "\n",
    "$$T_N(x) = \\sum^N_{n=0} \\frac{f^{(n)}(x_0)\\times \\Delta x^n}{n!}$$\n",
    "\n",
    "y $R_N(x)$ es el residuo (la parte de la serie que obviamos)\n",
    "\n",
    "$$ R_N(x) = \\frac{f^{(n+1)}(c) \\times \\Delta x^{n+1}}{(n+1)!} \\leq M \\Delta x^{n+1}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Ejemplo 1:***\n",
    "\n",
    "$f(x) = e^x$ con $x_0 = 0$\n",
    "\n",
    "Usando esto podemos encontrar expresiones para el error relativo y absoluto en función de $x$ asumiendo $N=2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Derivadas:\n",
    "$$\\begin{aligned}\n",
    "    f'(x) &= e^x \\\\\n",
    "    f''(x) &= e^x \\\\ \n",
    "    f^{(n)}(x) &= e^x\n",
    "\\end{aligned}$$\n",
    "\n",
    "Polinomio de Taylor:\n",
    "$$\\begin{aligned}\n",
    "    T_N(x) &= \\sum^N_{n=0} e^0 \\frac{x^n}{n!} \\Rightarrow \\\\\n",
    "    T_2(x) &= 1 + x + \\frac{x^2}{2}\n",
    "\\end{aligned}$$\n",
    "\n",
    "Restos:\n",
    "$$\\begin{aligned}\n",
    "    R_N(x) &= e^c \\frac{x^{n+1}}{(n+1)!} = e^c \\times \\frac{x^3}{6} \\quad \\Rightarrow \\\\\n",
    "    R_2(x) &\\leq \\frac{e^1}{6} \\approx 0.5\n",
    "\\end{aligned}$$\n",
    "\n",
    "Precisión:\n",
    "$$\n",
    "    e^1 = 2.718\\ldots \\\\\n",
    "    T_2(1) = 2.5 \\Rightarrow e \\approx 0.2 ~~ r \\approx 0.1\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "¡También podemos usar el paquete `sympy` que tiene la capacidad de calcular el polinomio de *Taylor* integrado!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = sympy.symbols('x')\n",
    "f = sympy.symbols('f', cls=sympy.Function)\n",
    "\n",
    "f = sympy.exp(x)\n",
    "f.series(x0=0, n=11)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 1.1**500\n",
    "\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Graficando"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-5, 5, 100)\n",
    "T_N = 1.0 + x + x**2 / 2.0 + x**3 / 6.0 + x**4 / 24.0 + x**5 / 120.0 + x**6 / 720.0 + x**7 / 5040.0 + x**8 / 40320.0 + x**9 / 362880\n",
    "R_N = np.exp(1) * x**10 / 3628800.0\n",
    "\n",
    "plt.plot(x, T_N, 'r', x, np.exp(x), 'k', x, R_N, 'b')\n",
    "plt.plot(0.0, 1.0, 'o', markersize=10)\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"$f(x)$, $T_N(x)$, $R_N(x)$\")\n",
    "plt.legend([\"$T_N(x)$\", \"$f(x)$\", \"$R_N(x)$\"], loc=2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R_N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Ejemplo 2:***\n",
    "\n",
    "Aproximar\n",
    "\n",
    "$$ f(x) = \\frac{1}{x} \\quad x_0  = 1,$$\n",
    "\n",
    "usando $x_0 = 1$ para el tercer termino de la serie de Taylor."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{aligned}\n",
    "    f'(x) &= -\\frac{1}{x^2} \\\\\n",
    "    f''(x) &= \\frac{2}{x^3} \\\\\n",
    "    f^{(n)}(x) &= \\frac{(-1)^n n!}{x^{n+1}}\n",
    "\\end{aligned}$$\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    T_N(x) &= \\sum^N_{n=0} (-1)^n (x-1)^n \\Rightarrow \\\\\n",
    "    T_2(x) &= 1 - (x - 1) + (x - 1)^2\n",
    "\\end{aligned}$$\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    R_N(x) &= \\frac{(-1)^{n+1}(x - 1)^{n+1}}{c^{n+2}} \\Rightarrow \\\\\n",
    "    R_2(x) &= \\frac{-(x - 1)^{3}}{c^{4}}\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "x = np.linspace(0.8, 2, 100)\n",
    "T_N = 1.0 - (x-1) + (x-1)**2\n",
    "R_N = -(x-1.0)**3 / (1.1**4)\n",
    "\n",
    "plt.plot(x, T_N, 'r', x, 1.0 / x, 'k', x, R_N, 'b')\n",
    "plt.plot(1.0, 1.0, 'o', markersize=10)\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"$f(x)$, $T_N(x)$, $R_N(x)$\")\n",
    "\n",
    "plt.legend([\"$T_N(x)$\", \"$f(x)$\", \"$R_N(x)$\"], loc=8)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Laboratorio Numérico 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert alert-success\">\n",
    "$\\color{red}{\\textbf{Ejercicio:}}$ Realice la expansión de la serie de Taylor para los dos ejemplos anteriores con 3, 4 y 5 términos. Cuál es el error que se tiene a medida que se adicionan más términos? Realice una gráfica comparativa del Residuo que se obtiene para cada término adicional. Haga un análisis de lo que sucede. Si extendemos hasta el \"infinito\" dicho residuo qué pueden concluir?\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error de punto flotante"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Errores surgen de aproximar números reales con números de precisión finita\n",
    "\n",
    "$$\\pi \\approx 3.14$$\n",
    "\n",
    "o $\\frac{1}{3} \\approx 0.333333333$ en decimal, los resultados forman un número finito de registros para representar cada número."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Ej.:*** considere la representación de $\\sqrt{2}=1.4142 \\ldots$. Como sabemos, éste es un número irracional, es decir, tiene una cantidad infinita de dígitos decimales. El computador almacena de forma incompleta la represenación de ese valor empleando cierta cantidad de números decimales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$2 - (\\sqrt{2})^2$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "a = np.sqrt(2)\n",
    "print(\"a: \", a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "b = abs(2-a**2)\n",
    "print(\"|2-a^2| = \", b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aritmética de punto flotante"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Los números en sistemas de [punto flotante](https://en.wikipedia.org/wiki/Floating-point_arithmetic \"Floating point arithmetic\") se representan como una serie de bits que representan diferentes partes de un número. En los sistemas de punto flotante normalizados, existen algunas convenciones estándar para el uso de estos bits. En general, los números se almacenan dividiéndolos en la forma\n",
    "\n",
    "$$fl(x) = \\pm (0.d_1 d_2 d_3 \\ldots d_p)_\\beta \\times \\beta^E$$\n",
    "\n",
    "donde los digitos $\\{d_i\\}_{i=1}^p$ son enteros tales que $0\\leq d_i \\leq \\beta-1$ y $d_1 \\neq 0$\n",
    "\n",
    "El sistema se caracteriza por cuatro números enteros:\n",
    "\n",
    "- la *base* $\\beta>1$. Para el sistema binario $\\beta = 2$, para decimal $\\beta = 10$, etc.\n",
    "\n",
    "\n",
    "- La precisión $p \\geq 1$, que representa la cantidad de dígitos significativos, y\n",
    "\n",
    "\n",
    "- el *exponente* $E$, que es un entero en el rango $[E_{\\min}, E_{\\max}]$\n",
    "\n",
    "\n",
    "$\\pm$ es un bit único y representa el signo del número.\n",
    "\n",
    "Los puntos importantes en cualquier sistema de punto flotante son:\n",
    "\n",
    "1. Existe un conjunto discreto y finito de números representables.\n",
    "\n",
    "\n",
    "2. Estos números representables no están distribuidos uniformemente en la línea real\n",
    "\n",
    "\n",
    "3. La aritmética en sistemas de punto flotante produce resultados diferentes de la aritmética de precisión infinita (es decir, matemática \"real\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notación de punto flotante"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Es común encontrar la siguiente notación para representar un conjunto de números de punto flotante:\n",
    "\n",
    "$$F(\\beta, p, E_{min}, E_{max})$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Propiedades de los sistemas de punto flotante"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Todos los sistemas de punto flotante se caracterizan por varios números importantes\n",
    "\n",
    "- Número normalizado reducido ([*underflow*](https://en.wikipedia.org/wiki/Arithmetic_underflow) si está por debajo, relacionado con números sub-normales alrededor de cero)\n",
    "\n",
    "\n",
    "- Número normalizado más grande ([*overflow*](https://en.wikipedia.org/wiki/Integer_overflow))\n",
    "\n",
    "\n",
    "- Cero\n",
    "\n",
    "\n",
    "- $\\epsilon$ o $\\epsilon_{mach}$\n",
    "\n",
    "\n",
    "- `Inf` y `nan`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Ejemplo: Sistema de juguete***\n",
    "\n",
    "Considere el sistema decimal de 2 digitos de  precisión (normalizado)\n",
    "\n",
    "$$F(10,2,-2,0)$$\n",
    "\n",
    "$$f = \\pm 0.d_1d_2 \\times 10^E$$\n",
    "\n",
    "con $E \\in [-2, 0]$.\n",
    "\n",
    "**Numero y distribución de números**\n",
    "\n",
    "\n",
    "1. Cuántos números pueden representarse con este sistema?\n",
    "\n",
    "\n",
    "2. Cuál es la distribución en la línea real?\n",
    "\n",
    "\n",
    "3. Cuáles son los límites underflow y overflow?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cuántos números pueden representarse con este sistema?\n",
    "\n",
    "$$f = \\pm 0.d_1d_2 \\times 10^E ~~~ \\text{con} ~~~ E \\in [-2, 0]$$\n",
    "\n",
    "$$2 \\times 9 \\times 10 \\times 3 + 1 = 541$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cuál es la distribución en la recta \"real\"?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAABzCAYAAABaW0yzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAQVElEQVR4nO3de5AlZXnH8e8PkCCysOKK7nLZFaUsL0VKgkgMCZtoGUEJXmLUGJCIJsQyYqlRVFQqaDQaL2USoxEoERWSEC/rrQJeiBoUXQjXIAgIsrBAuIPxgvjkj+5hm9kzM2d2embOGb6fqlPTp99+337e9/TpeebtnnNSVUiSJGnutlrsACRJkpYKEytJkqSemFhJkiT1xMRKkiSpJyZWkiRJPTGxkiRJ6omJlbQEJPlIkrf21NYeSe5OsnX7/KwkL++j7Un7uTvJnpPWbZXk80le1uN+Pp7kHX21N+Q+1ySpJNss5H4lLT4TK2nEJbk6yU+T3JXk9iRnJzkqyX3v36o6qqqOH7Ktp0+3TVX9uKp2qKp7+4h/mv3sUFVXTVr9TuBrVXXSfO57GEl+MCjBS3J0kvWLEZOk0WdiJY2HQ6pqGbAaeDfwRuDEvney2DMsVfWmqvrQYsbQcTJw+ID1h7Vl82KxXwNJc2NiJY2RqrqjqtYBLwRemuSJcP/LXUlWJPliO7t1a5JvtZfYTgH2AL7QXoZ7Q+eS1ZFJfgx8fYrLWI9O8r0kd7SX6nZu97U2yYZujN1ZsSRbJ3lzkivbGbdzk+zellWSx7TLOyX5RJL/TXJNkmMnZuSSHJHk20n+LsltSX6U5KCpxijJk5Kc1+7vX4DtJpU/O8n5ndm/vado6hTggCSrO3UfB+wNnJrkWUn+O8mdSa5Nctw0Ma1Ksq59Pa5I8opO2XFJTk/yySR3Ake043Fiko1Jrkvyjs6l2cck+c/2tbi57aOkEWFiJY2hqvoesAH47QHFr2vLHg48AnhzU6UOA35MM/u1Q1W9p1PnQOBxwO9PscvDgZcBq4BfAsPOKr0WeDFwMLBj28b/Ddju74GdgD3bWA4H/rRT/hTgMmAF8B7gxCSZ3EiSbYHP0SRFOwP/Bjy/U74PcBLw58DDgI8C65L82uS2qmoD8A2aGaoJhwNfrqqbgZ+0z5cDzwL+IslzphiHU2lek1XAHwJ/k+RpnfJDgdPbtj5FMyP2S+AxwJOAZwAT97kdD5wBPBTYjWbsJI0IEytpfF1PkzxMdg+wElhdVfdU1bdq5i8FPa6qflJVP52i/JSquriqfgK8FfijiRmUGbwcOLaqLqvGBVV1S3eDtp0XAm+qqruq6mrgfdw/obmmqj7W3vd1ctu/RwzY3/7Ag4APtn0/Hfh+p/wVwEer6pyqureqTgZ+3tYb5OSJONoZtJe066iqs6rqoqr6VVVdSJM8HTi5gXaG7gDgjVX1s6o6HzhhUv++U1Wfq6pf0SSgBwGvaV+Tm4APAC9qt72H5pLwqra9b08Ru6RFYGIlja9dgVsHrH8vcAVwRpKrkhwzRFvXzqL8GprkZcUQ7e4OXDnDNiuAbdt2u/vYtfP8homFqpqY8dphQFurgOsmJZLddlcDr2svA96e5PY2xlVTxPYZYGWS/YG1wPbAlwCSPCXJN9rLl3cARzF4TFYBt1bVXdP0rzu+q2nGd2Mnxo8Cu7TlbwACfC/JJYNusJe0eEyspDGU5Mk0v5g3m61oZ31eV1V7AocAr+1cdppq5mqmGa3dO8t70MyaTFwO274T19Y0lyAnXAs8eoa2b2bTLEx3H9fNUG+QjcCuky4T7jEpnndW1fLOY/uqOnVQY20SdzrNJb/DgNOq6hdt8aeBdcDuVbUT8BGahGey64Gdkyybpn/d8b+WZhZtRSfGHavqCW1MN1TVK6pqFc0lzQ9P3KsmafGZWEljJMmOSZ4NnAZ8sqouGrDNs9sbnAPcCdzbPgBupLmPabb+JMnjk2wP/DVwentZ7nJgu/ZG7gcBxwLd+5VOAI5Pslcaeyd5WLfhtp1/Bd6ZZFl7s/hrgU9uQZzfobk36dVJtknyPGC/TvnHgKPa2aYkeUgb+7KBrTVOprlU+Xzu/9+Ay2hmon6WZD/gjwdVrqprgbOBdyXZrr1Z/kiae6kGbb+R5h6q97Wv91ZJHp3kQIAkL0iyW7v5bTRJ2bx+NIak4ZlYSePhC0nuopnNeAvwfu5/c3fXXsBXgbtpEo0PV9VZbdm7gGPbS0yvn8X+TwE+TnNJbjvg1dD8lyLwSpoE6jqaGazufwm+nyZpOoMmyTsRePCA9v+yrXsVzSzcp2luMp+VdjbpecARNEnHC2ku502Ur6e5z+of2vIr2m2n803gDppLjN37tV4J/HX7uryNpp9TeTGwhmb26rPA26vqzGm2P5zm8uj/tHGeTnNfGcCTgXOS3E0zY3Z0Vf1ohj5IWiCZ+Z5WSZIkDcMZK0mSpJ6YWEmSJPXExEqSJKknJlaSJEk9mdcv+1yxYkWtWbNmPnchSZLUi3PPPffmqnr4zFtObV4TqzVr1rB+/fr53IUkSVIvklwz81bT81KgJElST0ysJEmSemJiJUmS1BMTK0mSpJ6YWEmSJPXExEqSJKknJlaSJEk9MbGSJEnqiYmVJElST0ysJEmSemJiJUmS1BMTK0mSpJ6YWEmSJPXExEqSJKknJlaSJEk9MbGSJEnqyfwmVtdf3/w87rhN62azPF/bjku9UY5tS+uNcmwLXW+UY5ttvTVrYPny5jGxPN26Ldl+7dr7dveBMy/fFENf7U+3/bBj1V0exddpHGLb0nqjHNtC1xvl2La03gLFthusYo5SVXNtY0r7JrW+ChKY2M9sludr23GpN8qx2SfHYvLyQmj3t+aYL3H1u5+14PsdqfH22Buf2ByLsenTvgnrq+Z0YtlqLpUlSZK0iYmVJElST0ysJEmSemJiJUmS1BMTK0mSpJ6YWEmSJPXExEqSJKkn28xn45cBa9eubZ5M/Jzt8nxtOy71Rjm2La03yrEtdL1Rjm229RZCu78brrqFtd9974Lvd9bLo1xvlGPb0nqjHNtC1xvl2La03kLFNkfz+gGhy5L6jXlrXdJSs2HHXdiw0y69tbfbHTex25039daepKXtbpjzB4T6yeujXG+UY7NPjsXk5YXQ7s9PXvfYG6vYHIux6ZOfvC5JkjRCTKwkSZJ6YmIlSZLUExMrSZKknphYSZIk9cTESpIkqScmVpIkST0xsZIkSerJvH6lDStXNj/f/vZN62azPF/bjku9UY5tS+uNcmwLXW+UY5ttvdWr4fbbm+Xly5vliZ+D1k1XNtX2a9bct7ujn7bXphg++MF+2p9u+5n6Py6v0zjEtqX1Rjm2ha43yrFtab0Fiu1G2Mgcze8nr++7b61fv37e2pckSepLknOrat+5tOGlQEmSpJ6YWEmSJPXExEqSJKknJlaSJEk9MbGSJEnqiYmVJElST0ysJEmSemJiJUmS1BMTK0mSpJ6YWEmSJPXExEqSJKknJlaSJEk9MbGSJEnqiYmVJElST0ysJEmSepKqmr/Gk7uAy+ZtB+NpBXDzYgcxghyXwRyXwRyXzTkmgzkugzkugz22qpbNpYFt+opkCpdV1b7zvI+xkmS9Y7I5x2Uwx2Uwx2VzjslgjstgjstgSdbPtQ0vBUqSJPXExEqSJKkn851Y/fM8tz+OHJPBHJfBHJfBHJfNOSaDOS6DOS6DzXlc5vXmdUmSpAcSLwVKkiT1xMRKkiSpJ3NKrJK8IMklSX6VZMp/20zyzCSXJbkiyTGd9TsnOTPJD9ufD51LPKNimH4leWyS8zuPO5O8pi07Lsl1nbKDF7wT82DY1zvJ1Ukuavu+frb1x82Qx8vuSb6R5NL2PXd0p2zJHC9TnSs65Unyobb8wiT7DFt3nA0xLi9px+PCJGcn+fVO2cD301IwxLisTXJH573xtmHrjqshxuSvOuNxcZJ7k+zcli3lY+WkJDcluXiK8v7OLVW1xQ/gccBjgbOAfafYZmvgSmBPYFvgAuDxbdl7gGPa5WOAv51LPKPymG2/2jG6AVjdPj8OeP1i92OxxgW4Glgx13Edl8cw/QJWAvu0y8uAyzvvoyVxvEx3ruhsczDwFSDA/sA5w9Yd18eQ4/JU4KHt8kET49I+H/h+GvfHkOOyFvjiltQdx8ds+wUcAnx9qR8rbd9+B9gHuHiK8t7OLXOasaqqS6tqpk9W3w+4oqquqqpfAKcBh7ZlhwInt8snA8+ZSzwjZLb9ehpwZVVdM59BjYC5vt4P2OOlqjZW1Xnt8l3ApcCuCxXgApnuXDHhUOAT1fgusDzJyiHrjqsZ+1ZVZ1fVbe3T7wK7LXCMi2Eur/lSPV5m268XA6cuSGSLrKq+Cdw6zSa9nVsW4h6rXYFrO883sOkXwiOqaiM0vziAXRYgnoUw2369iM0P7le105EnLZVLXgw/LgWckeTcJH+2BfXHzaz6lWQN8CTgnM7qpXC8THeumGmbYeqOq9n27Uiav7wnTPV+GnfDjstvJrkgyVeSPGGWdcfN0P1Ksj3wTODfO6uX6rEyjN7OLTN+pU2SrwKPHFD0lqr6/IyhNtNqk439ZzxMNy6zbGdb4A+AN3VW/xNwPM04HQ+8D3jZlkW6sHoal9+qquuT7AKcmeQH7V8bY6vH42UHmhPha6rqznb12B4vkwxzrphqmyV5nmkN3bckv0uTWB3QWb3k3k+tYcblPJpbLO5u7z38HLDXkHXH0Wz6dQjwX1XVncVZqsfKMHo7t8yYWFXV04cMaiobgN07z3cDrm+Xb0yysqo2tlNuN81xXwtmunFJMpt+HQScV1U3dtq+bznJx4Av9hHzQuhjXKrq+vbnTUk+SzMV+00e4MdLkgfRJFWfqqrPdNoe2+NlkunOFTNts+0QdcfVMONCkr2BE4CDquqWifXTvJ/G3Yzj0vnjg6r6cpIPJ1kxTN0xNZt+bXalZAkfK8Po7dyyEJcCvw/sleRR7ezMi4B1bdk64KXt8kuBYWbAxsFs+rXZNe72l+uE5wID/4thDM04LkkekmTZxDLwDDb1/wF7vCQJcCJwaVW9f1LZUjlepjtXTFgHHN7+B8/+wB3t5dNh6o6rGfuWZA/gM8BhVXV5Z/1076dxN8y4PLJ975BkP5rfebcMU3dMDdWvJDsBB9I51yzxY2UY/Z1b5niX/XNpsryfAzcC/9GuXwV8edLd9pfT3Fn/ls76hwFfA37Y/tx5LvGMymOqfg0Yl+1p3uQ7Tap/CnARcGH7Aq5c7D4t1LjQ/OfFBe3jEo+X+8blAJrp5wuB89vHwUvteBl0rgCOAo5qlwP8Y1t+EZ3/Rp7qPLMUHkOMywnAbZ1jY327fsr301J4DDEur2r7fQHNTf1PXerHy0xj0j4/AjhtUr2lfqycCmwE7qHJW46cr3OLX2kjSZLUEz95XZIkqScmVpIkST0xsZIkSeqJiZUkSVJPTKwkSZJ6YmIlSZLUExMrSZKknphYSRoLSZ7cftH0du2nRF+S5ImLHZckdfkBoZLGRpJ3ANsBDwY2VNW7FjkkSbofEytJY6P9rq7vAz+j+XqSexc5JEm6Hy8FShonOwM7AMtoZq4kaaQ4YyVpbCRZB5wGPIrmy6ZftcghSdL9bLPYAUjSMJIcDvyyqj6dZGvg7CS/V1VfX+zYJGmCM1aSJEk98R4rSZKknphYSZIk9cTESpIkqScmVpIkST0xsZIkSeqJiZUkSVJPTKwkSZJ68v85mWE7RUfpkAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x72 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "d_1_values = [1, 2, 3, 4, 5, 6, 7, 8, 9]\n",
    "d_2_values = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]\n",
    "E_values = [0, -1, -2]\n",
    "\n",
    "fig = plt.figure(figsize=(10.0, 1.0))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "for E in E_values:\n",
    "    for d1 in d_1_values:\n",
    "        for d2 in d_2_values:\n",
    "            axes.plot( (d1 + d2 * 0.1) * 10**E, 0.0, 'r+', markersize=20)\n",
    "            axes.plot(-(d1 + d2 * 0.1) * 10**E, 0.0, 'r+', markersize=20)\n",
    "            \n",
    "axes.plot(0.0, 0.0, '+', markersize=20)\n",
    "axes.plot([-10.0, 10.0], [0.0, 0.0], 'k')\n",
    "\n",
    "axes.set_title(\"Distribución de Valores\")\n",
    "axes.set_yticks([])\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"\")\n",
    "axes.set_xlim([-1., 1.0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cuáles son los límites superior (overflow) e inferior (underflow)?\n",
    "\n",
    "- El menor número que puede ser representado (underflow) es: $0.10 \\times 10^{-2} = 0.01$\n",
    "\n",
    "\n",
    "- El mayor número que puede ser representado (overflow) es:  $0.99 \\times 10^0 = 0.99$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sistema Binario"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considere el sistema en base 2 de 2 dígitos de precisión\n",
    "\n",
    "$$F(2,2,-1,1)$$\n",
    "\n",
    "$$F=\\pm 0.d_1d_2 \\times 2^E \\quad \\text{con} \\quad E \\in [-1, 1]$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Numero y distribución de números"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Cuántos números pueden representarse con este sistema?\n",
    "\n",
    "\n",
    "2. Cuál es la distribución en la línea real?\n",
    "\n",
    "\n",
    "3. Cuáles son los límites underflow y overflow?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cuántos números pueden representarse en este sistema?\n",
    "\n",
    "\n",
    "$$f=\\pm 0.d_1d_2 \\times 2^E ~~~~ \\text{con} ~~~~ E \\in [-1, 1]$$\n",
    "\n",
    "$$ 2 \\times 1 \\times 2 \\times 3 + 1 = 13$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cuál es la distribución en la línea real?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAABzCAYAAAB3qWMJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAANlklEQVR4nO3de5CddX3H8feHAEZIwGKoQyAkXhhG6zC1RHQ6tGYq9YqXaq1aa7Qq09Sx4ogjiiiUS2m91akdR6phjMilFtTGamfQShBHuWwsXhDiIBITRBGRQLxUAt/+8TzbHtZNdpfd5Tlnn/drZmf2nOc5v+d7vic755Pf73d2U1VIkiT1xV5dFyBJkvRQMvxIkqReMfxIkqReMfxIkqReMfxIkqReMfxIkqReMfxIQyjJh5O8c47GOjzJziSL2tubkrxuLsaecJ2dSR4z4b69kvx7ktfM4XU+luSsuRpvmtdclaSS7P1QXlfS/DD8SA+xJLck+WWSe5LcleSrSdYl+b+fx6paV1VnTnOs4/Z0TlX9oKqWVNV9c1H/Hq6zpKpunnD32cB/VdV583nt6Uhy42QhLMmJSca6qElSNww/UjeeV1VLgZXA3wMnA+vn+iJdz1RU1dur6p+6rGHABmDtJPe/sj02L7p+DST9JsOP1KGq2lFVG4GXAq9K8kR44NJOkmVJ/qOdJbozyZXtctL5wOHAZ9slp7cOLM+8NskPgC/tZsnmsUmuSbKjXZY6qL3WmiTbB2scnF1KsijJKUm+185cbU6yoj1WSR7Xfn9gko8n+UmSrUlOHZ/ZSvLqJF9J8t4kP0vy/STP3l2Pkjwpydfb6/0rsHjC8eOTXDcwi3bUboY6Hzg2ycqBxz4eOAq4KMlzk/x3kruTbEty+h5qWp5kY/t63JTkhIFjpye5JMknktwNvLrtx/oktyW5NclZA8uQj0tyRfta3NE+R0nzyPAjDYGqugbYDvzBJIdPao8dDDwKOKV5SL0S+AHNLNKSqnr3wGOeBjweeOZuLrkWeA2wHNgFTHd25s3Ay4HnAAe0Y/xikvM+CBwIPKatZS3wlwPHnwJsAZYB7wbWJ8nEQZLsC3yGJrgcBPwb8OKB478HnAf8FfBI4FxgY5KHTRyrqrYDl9PM9IxbC3y+qu4Aft7efgTwXOCvk7xwN324iOY1WQ78KfB3SZ4+cPwFwCXtWBfQzCztAh4HPAl4BjC+7+pM4DLgt4DDaHonaR4ZfqTh8UOaN/iJ7gUOAVZW1b1VdWVN/Uf5Tq+qn1fVL3dz/Pyq+nZV/Rx4J/Bn4zMRU3gdcGpVbanGN6rqp4MntOO8FHh7Vd1TVbcA7+OBoWNrVX2k3Ye0oX1+j5rkek8F9gE+0D73S4BrB46fAJxbVVdX1X1VtQH4n/Zxk9kwXkc7E/WK9j6qalNVfauq7q+qb9IEnKdNHKCd6ToWOLmqflVV1wEfnfD8vlZVn6mq+2lC4rOBN7Wvye3APwIva8+9l2b5c3k73ld2U7ukOWL4kYbHocCdk9z/HuAm4LIkNyd52zTG2jaD41tpAsayaYy7AvjeFOcsA/Ztxx28xqEDt380/k1Vjc8cLZlkrOXArRPC3uC4K4GT2iWvu5Lc1da4fDe1fQo4JMlTgTXAfsDnAJI8Jcnl7VLdDmAdk/dkOXBnVd2zh+c32N+VNP29baDGc4Hfbo+/FQhwTZLrJ9uULWluGX6kIZDkyTRvnr/xv/529uSkqnoM8DzgzQNLLLubAZpqZmjFwPeH08w+jC/97DdQ1yKa5bZx24DHTjH2Hfz/bMbgNW6d4nGTuQ04dMKS2OET6jm7qh4x8LVfVV002WBt0LqEZnnrlcDFVfXr9vCFwEZgRVUdCHyYJpRM9EPgoCRL9/D8Bvu/jWY2atlAjQdU1e+0Nf2oqk6oquU0y3cfGt87JWl+GH6kDiU5IMnxwMXAJ6rqW5Occ3y7KTbA3cB97RfAj2n21czUXyR5QpL9gDOAS9olqO8Ci9vNv/sApwKD+2c+CpyZ5Ig0jkryyMGB23E+CZydZGm7wfjNwCceRJ1fo9kr88Ykeyd5EXDMwPGPAOvaWZsk2b+tfemkozU20CzLvZgHfsprKc2Mzq+SHAP8+WQPrqptwFeBc5IsbjdYv5Zmb89k599Gs6fnfe3rvVeSxyZ5GkCSlyQ5rD39ZzTBaV5/LYHUd4YfqRufTXIPzazAO4D388ANwYOOAL4I7KQJAx+qqk3tsXOAU9vllLfM4PrnAx+jWX5aDLwRmk+fAa+nCTm30swEDX766/00weYymiC2Hnj4JOP/TfvYm2lmsy6k2Zg8I+2szIuAV9MEg5fSLF2NHx+j2ffzz+3xm9pz9+TLwA6a5bTB/UOvB85oX5d30TzP3Xk5sIpmFujTwGlV9YU9nL+WZinwO22dl9DscwJ4MnB1kp00M08nVtX3p3gOkmYhU++blCRJWjic+ZEkSb1i+JEkSb1i+JEkSb1i+JEkSb0yoz+4t2zZslq1atU8lSJJkjR3Nm/efEdVHTzx/hmFn1WrVjE2NjZ3VUmSJM2TJFsnu99lL0mS1CuGH0mS1CuGH0mS1CuGH0mS1CuGH0mS1CuGH0mS1CuGH0mS1CuGH0mS1CuGH0mS1CuGH0mS1CuGH0mS1CuGH0mS1CuGH0mS1CuGH0mS1CuGH0mS1CuGH0mS1CsLK/ycfnrXFYyuYe9d3+ob9vHmWp/q87UdLsNe37Ab0f6lqqZ98urVq2tsbGwey5mlBGbwfDRg2HvXt/qGfby51qf6fG2Hy7DXN+yGvH9JNlfV6on3L6yZH0mSpCkYfiRJUq8YfiRJUq8YfiRJUq8YfiRJUq8YfiRJUq8YfiRJUq/sPZOTt2zZwpo1a+aplDky7PUNs2HvXd/qG/bx5lqf6vO1HS7DXt+wG8H+zeiXHC5durSOPvroeSxnN265BbZunbvxVq6EVavmbrxhNuy961t9Bx4IO3YM73jD3r9hr28u+drOzrDXN+wWSP+uuOKKSX/J4Yxmfo488kg2bdo0Z0XNuSH/TZNDbdh717f6hn28udan+nxth8uw1zfshrx/SSa93z0/kiSpVww/kiSpVww/kiSpVww/kiSpVww/kiSpVww/kiSpVww/kiSpVww/kiSpVxZW+DnttK4rGF3D3ru+1Tfs4821PtXnaztchr2+YTei/ZvRn7dYvXp1jY2NzWM5kiRJcyPJpH/eYmHN/EiSJE3B8CNJknrF8CNJknrF8CNJknrF8CNJknrF8CNJknrF8CNJknrF8CNJknrF8CNJknrF8CNJknrF8CNJknrF8CNJknrF8CNJknrF8CNJknrF8CNJknolVTX9k5OfAFvnr5xZWwbc0XURI8rezY79mx379+DZu9mxf7Mz7P1bWVUHT7xzRuFn2CUZq6rVXdcxiuzd7Ni/2bF/D569mx37Nzuj2j+XvSRJUq8YfiRJUq8stPDzL10XMMLs3ezYv9mxfw+evZsd+zc7I9m/BbXnR5IkaSoLbeZHkiRpjww/kiSpVxZU+ElyZpJvJrkuyWVJlndd0yhJ8p4kN7Y9/HSSR3Rd0yhJ8pIk1ye5P8nIffSzC0melWRLkpuSvK3rekZJkvOS3J7k213XMoqSrEhyeZIb2p/bE7uuaVQkWZzkmiTfaHv3t13XNFMLas9PkgOq6u72+zcCT6iqdR2XNTKSPAP4UlXtSvIPAFV1csdljYwkjwfuB84F3lJVYx2XNNSSLAK+C/wxsB24Fnh5VX2n08JGRJI/BHYCH6+qJ3Zdz6hJcghwSFV9PclSYDPwQv/9TS1JgP2rameSfYCvACdW1VUdlzZtC2rmZzz4tPYHFk6yewhU1WVVtau9eRVwWJf1jJqquqGqtnRdxwg5Bripqm6uql8DFwMv6LimkVFVXwbu7LqOUVVVt1XV19vv7wFuAA7ttqrRUI2d7c192q+Rer9dUOEHIMnZSbYBrwDe1XU9I+w1wH92XYQWtEOBbQO3t+ObjzqQZBXwJODqjksZGUkWJbkOuB34QlWNVO9GLvwk+WKSb0/y9QKAqnpHVa0ALgDe0G21w2eq/rXnvAPYRdNDDZhO/zRtmeS+kfrfo0ZfkiXApcCbJqweaA+q6r6q+l2aFYJjkozU0uveXRcwU1V13DRPvRD4HHDaPJYzcqbqX5JXAccDT6+FtCFsjszg35+mth1YMXD7MOCHHdWiHmr3q1wKXFBVn+q6nlFUVXcl2QQ8CxiZzfcjN/OzJ0mOGLj5fODGrmoZRUmeBZwMPL+qftF1PVrwrgWOSPLoJPsCLwM2dlyTeqLdtLseuKGq3t91PaMkycHjnwZO8nDgOEbs/XahfdrrUuBImk/cbAXWVdWt3VY1OpLcBDwM+Gl711V+Wm76kvwJ8EHgYOAu4LqqemanRQ25JM8BPgAsAs6rqrO7rWh0JLkIWAMsA34MnFZV6zstaoQkORa4EvgWzXsGwClV9fnuqhoNSY4CNtD83O4FfLKqzui2qplZUOFHkiRpKgtq2UuSJGkqhh9JktQrhh9JktQrhh9JktQrhh9JktQrhh9JktQrhh9JktQrhh9JnUjy5CTfTLI4yf5Jrh+1vw8kaTT5Sw4ldSbJWcBi4OHA9qo6p+OSJPWA4UdSZ9q/6XUt8Cvg96vqvo5LktQDLntJ6tJBwBJgKc0MkCTNO2d+JHUmyUbgYuDRwCFV9YaOS5LUA3t3XYCkfkqyFthVVRcmWQR8NckfVdWXuq5N0sLmzI8kSeoV9/xIkqReMfxIkqReMfxIkqReMfxIkqReMfxIkqReMfxIkqReMfxIkqRe+V/tvWUGSosGfQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x72 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "d_1_values = [1]\n",
    "d_2_values = [0, 1]\n",
    "E_values = [1, 0, -1]\n",
    "\n",
    "fig = plt.figure(figsize=(10.0, 1.0))\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "\n",
    "for E in E_values:\n",
    "    for d1 in d_1_values:\n",
    "        for d2 in d_2_values:\n",
    "            axes.plot( (d1 + d2 * 0.5) * 2**E, 0.0, 'r+', markersize=20)\n",
    "            axes.plot(-(d1 + d2 * 0.5) * 2**E, 0.0, 'r+', markersize=20)\n",
    "            \n",
    "axes.plot(0.0, 0.0, 'r+', markersize=20)\n",
    "axes.plot([-4.5, 4.5], [0.0, 0.0], 'k')\n",
    "\n",
    "axes.set_title(\"Distribución de Valores\")\n",
    "axes.set_yticks([])\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"\")\n",
    "axes.set_xlim([-3.5, 3.5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cuáles son los límites superior (*overflow*) e inferior (*underflow*)?\n",
    "\n",
    "- El menor número que puede ser representado (*underflow*) es: $1.0 \\times 2^{-1} = 0.5$\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "- El mayor número que puede ser representado (*overflow*) es:  $1.1 \\times 2^1 = 3$\n",
    "\n",
    "Observe que estos números son en sistema binario. \n",
    "\n",
    "Una rápida regla de oro:\n",
    "\n",
    "$$2^3 2^2 2^1 2^0 . 2^{-1} 2^{-2} 2^{-3}$$\n",
    "\n",
    "corresponde a\n",
    "\n",
    "8s, 4s, 2s, 1s . mitades, cuartos, octavos, $\\ldots$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Ejercicio:*** Cuál sería la representación en punto flotante del siguiente conjunto de números:\n",
    "\n",
    "$$F(2,3,-1,3)$$\n",
    "\n",
    "- Cuántos números se pueden representar?\n",
    "\n",
    "\n",
    "- Cuál sería el menor número representable?\n",
    "\n",
    "\n",
    "- Cuál sería el mayor número?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sistema real - [IEEE 754](https://en.wikipedia.org/wiki/IEEE_754) sistema binario de punto flotante"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Precisión simple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "- Almacenamiento total es de 32 bits\n",
    "\n",
    "\n",
    "- Exponente de 8 bits $\\Rightarrow E \\in [-126, 127]$\n",
    "\n",
    "\n",
    "- Fracción 23 bits ($p = 24$)\n",
    "\n",
    "\n",
    "```\n",
    "s EEEEEEEE FFFFFFFFFFFFFFFFFFFFFFF\n",
    "0 1      8 9                     31\n",
    "```\n",
    "\n",
    "Overflow $= 2^{127} \\approx 3.4 \\times 10^{38}$\n",
    "\n",
    "Underflow $= 2^{-126} \\approx 1.2 \\times 10^{-38}$\n",
    "\n",
    "$\\epsilon_{\\text{machine}} = 2^{-23} \\approx 1.2 \\times 10^{-7}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Precisión doble"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Almacenamiento total asignado es 64 bits\n",
    "\n",
    "- Exponenete de 11 bits $\\Rightarrow E \\in [-1022, 1024]$\n",
    "\n",
    "- Fracción de 52 bits ($p = 53$)\n",
    "\n",
    "```\n",
    "s EEEEEEEEEE FFFFFFFFFF FFFFFFFFFF FFFFFFFFFF FFFFFFFFFF FFFFFFFFFF FF\n",
    "0 1       11 12                                                      63\n",
    "```\n",
    "Overflow $= 2^{1024} \\approx 1.8 \\times 10^{308}$\n",
    "\n",
    "Underflow $= 2^{-1022} \\approx 2.2 \\times 10^{-308}$\n",
    "\n",
    "$\\epsilon_{\\text{machine}} = 2^{-52} \\approx 2.2 \\times 10^{-16}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Acceso de Python a números de la IEEE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Accede a muchos parámetros importantes, como el epsilon de la máquina\n",
    "\n",
    "```python\n",
    "import numpy as np\n",
    "np.finfo(float).eps\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "np.finfo(float).eps\n",
    "\n",
    "print(np.finfo(np.float16))\n",
    "print(np.finfo(np.float32))\n",
    "print(np.finfo(float))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculo \"manual\" del $\\epsilon_{mach}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La determinación \"manual\" del $\\epsilon_{mach}$ es muy simple. Veamos el siguiente algoritmo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "17\n",
      "1.0000000000000001e-16\n"
     ]
    }
   ],
   "source": [
    "eps = 1.0\n",
    "count = 1\n",
    "\n",
    "while 1.0 + eps > 1.0:\n",
    "    eps = eps / 10.0\n",
    "    count += 1\n",
    "    \n",
    "print(count)\n",
    "print(eps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si lo comparamos con el valor obtenido en el numeral anterior para `float64`, `2.2204460492503131e-16`, se observa que es del orden de dos veces menor, por qué?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Por qué debería importarnos esto?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C01_Img02.png?raw=true\" width=\"500\" />\n",
    "</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Aritmética de punto flotante no es conmutativa o asociativa\n",
    "\n",
    "\n",
    "- Errores de punto flotante compuestos, No asuma que la precisión doble es suficiente\n",
    "\n",
    "\n",
    "- Mezclar precisión es muy peligroso"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***EL ORDEN DE LOS FACTORES NO ALTERA EL PRODUCTO???***\n",
    "\n",
    "$$2 \\times 3 = 3 \\times 2 = 6$$\n",
    "\n",
    "\n",
    "$$ 10^{300} \\times 10^{50} \\times 10^{-60} = 10^{300} \\times 10^{-60} \\times 10^{50} ??$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 10**300\n",
    "b = 10**50\n",
    "c = 10**-60\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "OverflowError",
     "evalue": "int too large to convert to float",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mOverflowError\u001b[0m                             Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-9-8b26f45c3c82>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0md1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0ma\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mb\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mc\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0mprint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"d1: \"\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0md1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mOverflowError\u001b[0m: int too large to convert to float"
     ]
    }
   ],
   "source": [
    "d1 = a * b * c\n",
    "print(\"d1: \", d1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "d2:  1e+290\n"
     ]
    }
   ],
   "source": [
    "d2 = b * c * a\n",
    "print(\"d2: \", d2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ejemplo 1: Aritmética simple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aritmética simple $\\delta < \\epsilon_{\\text{machine}}$\n",
    "\n",
    "   $$(1+\\delta) - 1 = 1 - 1 = 0$$\n",
    "\n",
    "   $$1 - 1 + \\delta = \\delta$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n"
     ]
    }
   ],
   "source": [
    "delta = 1.000000000000001 * eps\n",
    "\n",
    "value = (1 + delta) - 1\n",
    "print(value)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ejemplo 2: Cancelación catastrófica"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Miremos qué sucede cuando sumamos dos números $x$ y $y$ cuando $x+y \\neq 0$. De hecho, podemos estimar estos límites haciendo un análisis de error. Aquí necesitamos presentar la idea de que cada operación de punto flotante introduce un error tal que\n",
    "\n",
    "$$\n",
    "    \\text{fl}(x ~\\text{op}~ y) = (x ~\\text{op}~ y) (1 + \\delta)\n",
    "$$\n",
    "\n",
    "donde $\\text{fl}(\\cdot)$ es una función que devuelve la representación de punto flotante de la expresión encerrada, $\\text{op}$ es alguna operación (ex. $+, -, \\times, /$), y $\\delta$ es el error de punto flotante debido a  $\\text{op}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "De vuelta a nuestro problema en cuestión. El error de coma flotante debido a la suma es\n",
    "\n",
    "$$\\text{fl}(x + y) = (x + y) (1 + \\delta).$$\n",
    "\n",
    "\n",
    "Comparando esto con la solución verdadera usando un error relativo tenemos\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\frac{(x + y) - \\text{fl}(x + y)}{x + y} &= \\frac{(x + y) - (x + y) (1 + \\delta)}{x + y} = \\delta.\n",
    "\\end{aligned}$$\n",
    "\n",
    "entonces si $\\delta = \\mathcal{O}(\\epsilon_{\\text{machine}})$ no estaremos muy preocupados."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Que pasa si consideramos un error de punto flotante en la representación de $x$ y $y$, $x \\neq y$, y decimos que $\\delta_x$ y $\\delta_y$ son la magnitud de los errores en su representación. Asumiremos que esto constituye el error de punto flotante en lugar de estar asociado con la operación en sí.\n",
    "\n",
    "Dado todo esto, tendríamos\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\text{fl}(x + y) &= x (1 + \\delta_x) + y (1 + \\delta_y) \\\\\n",
    "    &= x + y + x \\delta_x + y \\delta_y \\\\\n",
    "    &= (x + y) \\left(1 + \\frac{x \\delta_x + y \\delta_y}{x + y}\\right)\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculando nuevamente el error relativo, tendremos\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\frac{x + y - (x + y) \\left(1 + \\frac{x \\delta_x + y \\delta_y}{x + y}\\right)}{x + y} &= 1 - \\left(1 + \\frac{x \\delta_x + y \\delta_y}{x + y}\\right) \\\\\n",
    "    &= \\frac{x}{x + y} \\delta_x + \\frac{y}{x + y} \\delta_y \\\\\n",
    "    &= \\frac{1}{x + y} (x \\delta_x + y \\delta_y)\n",
    "\\end{aligned}$$\n",
    "\n",
    "Lo importante aquí es que ahora el error depende de los valores de $x$ y $y$, y más importante aún, su suma. De particular preocupación es el tamaño relativo de $x + y$. A medida que se acerca a cero en relación con las magnitudes de $x$ y $y$, el error podría ser arbitrariamente grande. Esto se conoce como ***cancelación catastrófica***."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = np.array([10**(-n) for n in range(1, 16)])\n",
    "x = 1.0 + dx\n",
    "y = -np.ones(x.shape)\n",
    "error = np.abs(x + y - dx) / (dx)\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.loglog(dx, x + y, 'o-')\n",
    "axes.set_xlabel(\"$\\Delta x$\")\n",
    "axes.set_ylabel(\"$x + y$\")\n",
    "axes.set_title(\"$\\Delta x$ vs. $x+y$\")\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.loglog(dx, error, 'o-')\n",
    "axes.set_xlabel(\"$\\Delta x$\")\n",
    "axes.set_ylabel(\"$|x + y - \\Delta x| / \\Delta x$\")\n",
    "axes.set_title(\"Diferencia entre $x$ y $y$ vs. Error relativo\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ejemplo 3: Evaluación de una función"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considere la función\n",
    "\n",
    "$$\n",
    "    f(x) = \\frac{1 - \\cos x}{x^2}\n",
    "$$\n",
    "\n",
    "con $x\\in[-10^{-4}, 10^{-4}]$. \n",
    "\n",
    "Tomando el límite cuando $x \\rightarrow 0$ podemos ver qué comportamiento esperaríamos ver al evaluar esta función:\n",
    "\n",
    "$$\n",
    "    \\lim_{x \\rightarrow 0} \\frac{1 - \\cos x}{x^2} = \\lim_{x \\rightarrow 0} \\frac{\\sin x}{2 x} = \\lim_{x \\rightarrow 0} \\frac{\\cos x}{2} = \\frac{1}{2}.\n",
    "$$\n",
    "\n",
    "¿Qué hace la representación de punto flotante?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = (1-np.cos(0))/0**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1e-3, 1e-3, 100, dtype=np.float32)\n",
    "error = (0.5 - (1.0 - np.cos(x)) / x**2) / 0.5\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, error, 'o')\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"Error Relativo\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ejemplo 4: Evaluación de un Polinomio"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " $$f(x) = x^7 - 7x^6 + 21 x^5 - 35 x^4 + 35x^3-21x^2 + 7x - 1$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "x = np.linspace(0.988, 1.012, 1000, dtype=np.float16)\n",
    "y = x**7 - 7.0 * x**6 + 21.0 * x**5 - 35.0 * x**4 + 35.0 * x**3 - 21.0 * x**2 + 7.0 * x - 1.0\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, y, 'r')\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"y\")\n",
    "axes.set_ylim((-0.1, 0.1))\n",
    "axes.set_xlim((x[0], x[-1]))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ejemplo 5: Evaluación de una función racional"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calcule $f(x) = x + 1$ por la función $$F(x) = \\frac{x^2 - 1}{x - 1}$$\n",
    "\n",
    "¿Cuál comportamiento esperarías encontrar?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0.5, 1.5, 101, dtype=np.float16)\n",
    "f_hat = (x**2 - 1.0) / (x - 1.0)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, np.abs(f_hat - (x + 1.0)))\n",
    "axes.set_xlabel(\"$x$\")\n",
    "axes.set_ylabel(\"Error Absoluto\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combinación de error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En general, nos debemos ocupar de la combinación de error de truncamiento con el error de punto flotante.\n",
    "\n",
    "- Error de Truncamiento: errores que surgen de la aproximación de una función, truncamiento de una serie.\n",
    "\n",
    "$$\\sin x \\approx x - \\frac{x^3}{3!} + \\frac{x^5}{5!} + O(x^7)$$\n",
    "\n",
    "\n",
    "- Error de punto flotante: errores derivados de la aproximación de números reales con números de precisión finita\n",
    "\n",
    "$$\\pi \\approx 3.14$$\n",
    "\n",
    "o $\\frac{1}{3} \\approx 0.333333333$ en decimal, los resultados forman un número finito de registros para representar cada número."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ejemplo 1:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considere la aproximación de diferencias finitas donde $f(x) = e^x$ y estamos evaluando en $x=1$\n",
    "\n",
    "$$f'(x) \\approx \\frac{f(x + \\Delta x) - f(x)}{\\Delta x}$$\n",
    "\n",
    "Compare el error entre disminuir $\\Delta x$ y la verdadera solucion $f'(1) = e$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta_x = np.linspace(1e-20, 5.0, 100)\n",
    "delta_x = np.array([2.0**(-n) for n in range(1, 60)])\n",
    "x = 1.0\n",
    "f_hat_1 = (np.exp(x + delta_x) - np.exp(x)) / (delta_x)\n",
    "f_hat_2 = (np.exp(x + delta_x) - np.exp(x - delta_x)) / (2.0 * delta_x)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.loglog(delta_x, np.abs(f_hat_1 - np.exp(1)), 'o-', label=\"Unilateral\")\n",
    "axes.loglog(delta_x, np.abs(f_hat_2 - np.exp(1)), 's-', label=\"Centrado\")\n",
    "axes.legend(loc=3)\n",
    "axes.set_xlabel(\"$\\Delta x$\")\n",
    "axes.set_ylabel(\"Error Absoluto\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ejemplo 2:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evalúe $e^x$ con la serie de *Taylor*\n",
    "\n",
    "$$e^x = \\sum^\\infty_{n=0} \\frac{x^n}{n!}$$\n",
    "\n",
    "podemos elegir $n< \\infty$ que puede aproximarse $e^x$ en un rango dado $x \\in [a,b]$ tal que el error relativo $E$ satisfaga $E<8 \\cdot \\varepsilon_{\\text{machine}}$?\n",
    "\n",
    "¿Cuál podría ser una mejor manera de simplemente evaluar el polinomio de Taylor directamente por varios $N$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def my_exp(x, N=10):\n",
    "    value = 0.0\n",
    "    for n in range(N + 1):\n",
    "        value += x**n / scipy.special.factorial(n)\n",
    "        \n",
    "    return value\n",
    "\n",
    "x = np.linspace(-2, 2, 100, dtype=np.float32)\n",
    "for N in range(1, 50):\n",
    "    error = np.abs((np.exp(x) - my_exp(x, N=N)) / np.exp(x))\n",
    "    if np.all(error < 8.0 * np.finfo(float).eps):\n",
    "        break\n",
    "\n",
    "print(N)\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.plot(x, error)\n",
    "axes.set_xlabel(\"x\")\n",
    "axes.set_ylabel(\"Error Relativo\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ejemplo 3: Error relativo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Digamos que queremos calcular el error relativo de dos valores $x$ y $y$ usando $x$ como valor de normalización\n",
    "\n",
    "$$\n",
    "    E = \\frac{x - y}{x}\n",
    "$$\n",
    "y\n",
    "$$\n",
    "    E = 1 - \\frac{y}{x}\n",
    "$$\n",
    "\n",
    "son equivalentes. En precisión finita, ¿qué forma pidría esperarse que sea más precisa y por qué?\n",
    "\n",
    "Ejemplo tomado de [blog](https://nickhigham.wordpress.com/2017/08/14/how-and-how-not-to-compute-a-relative-error/) posteado por Nick Higham*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usando este modelo, la definición original contiene dos operaciones de punto flotante de manera que\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    E_1 = \\text{fl}\\left(\\frac{x - y}{x}\\right) &= \\text{fl}(\\text{fl}(x - y) / x) \\\\\n",
    "    &= \\left[ \\frac{(x - y) (1 + \\delta_+)}{x} \\right ] (1 + \\delta_/) \\\\\n",
    "    &= \\frac{x - y}{x}  (1 + \\delta_+) (1 + \\delta_/)\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para la otra formulación tenemos\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    E_2 = \\text{fl}\\left( 1 - \\frac{y}{x} \\right ) &= \\text{fl}\\left(1 - \\text{fl}\\left(\\frac{y}{x}\\right) \\right) \\\\\n",
    "    &= \\left(1 - \\frac{y}{x} (1 + \\delta_/) \\right) (1 + \\delta_-)\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si suponemos que todos las $\\text{op}$s tienen magnitudes de error similares, entonces podemos simplificar las cosas dejando que \n",
    "\n",
    "$$\n",
    "    |\\delta_\\ast| \\le \\epsilon.\n",
    "$$\n",
    "\n",
    "Para comparar las dos formulaciones, nuevamente usamos el error relativo entre el error relativo verdadero $e_i$ y nuestras versiones calculadas $E_i$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Definición original\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\frac{e - E_1}{e} &= \\frac{\\frac{x - y}{x} - \\frac{x - y}{x}  (1 + \\delta_+) (1 + \\delta_/)}{\\frac{x - y}{x}} \\\\\n",
    "    &\\le 1 - (1 + \\epsilon) (1 + \\epsilon) = 2 \\epsilon + \\epsilon^2\n",
    "\\end{aligned}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Definición manipulada:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "    \\frac{e - E_2}{e} &= \\frac{e - \\left[1 - \\frac{y}{x}(1 + \\delta_/) \\right] (1 + \\delta_-)}{e} \\\\\n",
    "    &= \\frac{e - \\left[e - \\frac{y}{x} \\delta_/) \\right] (1 + \\delta_-)}{e} \\\\\n",
    "    &= \\frac{e - \\left[e + e\\delta_- - \\frac{y}{x} \\delta_/ - \\frac{y}{x} \\delta_/ \\delta_-)) \\right] }{e} \\\\\n",
    "    &= - \\delta_- + \\frac{1}{e} \\frac{y}{x} \\left(\\delta_/ + \\delta_/ \\delta_- \\right) \\\\\n",
    "    &= - \\delta_- + \\frac{1 -e}{e} \\left(\\delta_/ + \\delta_/ \\delta_- \\right) \\\\\n",
    "    &\\le \\epsilon + \\left |\\frac{1 - e}{e}\\right | (\\epsilon + \\epsilon^2)\n",
    "\\end{aligned}$$\n",
    "\n",
    "Vemos entonces que nuestro error de punto flotante dependerá de la magnitud relativa de $e$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Based on the code by Nick Higham\n",
    "# https://gist.github.com/higham/6f2ce1cdde0aae83697bca8577d22a6e\n",
    "# Compares relative error formulations using single precision and compared to double precision\n",
    "\n",
    "N = 501    # Note: Use 501 instead of 500 to avoid the zero value\n",
    "d = numpy.finfo(numpy.float32).eps * 1e4\n",
    "a = 3.0\n",
    "x = a * numpy.ones(N, dtype=numpy.float32)\n",
    "y = [x[i] + numpy.multiply((i - numpy.divide(N, 2.0, dtype=numpy.float32)), d, dtype=numpy.float32) for i in range(N)]\n",
    "\n",
    "# Compute errors and \"true\" error\n",
    "relative_error = numpy.empty((2, N), dtype=numpy.float32)\n",
    "relative_error[0, :] = numpy.abs(x - y) / x\n",
    "relative_error[1, :] = numpy.abs(1.0 - y / x)\n",
    "exact = numpy.abs( (numpy.float64(x) - numpy.float64(y)) / numpy.float64(x))\n",
    "\n",
    "# Compute differences between error calculations\n",
    "error = numpy.empty((2, N))\n",
    "for i in range(2):\n",
    "    error[i, :] = numpy.abs((relative_error[i, :] - exact) / numpy.abs(exact))\n",
    "\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "axes.semilogy(y, error[0, :], '.', markersize=10, label=\"$|x-y|/|x|$\")\n",
    "axes.semilogy(y, error[1, :], '.', markersize=10, label=\"$|1-y/x|$\")\n",
    "\n",
    "axes.grid(True)\n",
    "axes.set_xlabel(\"y\")\n",
    "axes.set_ylabel(\"Error Relativo\")\n",
    "axes.set_xlim((numpy.min(y), numpy.max(y)))\n",
    "axes.set_ylim((5e-9, numpy.max(error[1, :])))\n",
    "axes.set_title(\"Comparasión Error Relativo\")\n",
    "axes.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Algunos enlaces de utilidad con respecto al punto flotante IEEE:\n",
    "\n",
    "- [What Every Computer Scientist Should Know About Floating-Point Arithmetic](http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html)\n",
    "\n",
    "\n",
    "- [IEEE 754 Floating Point Calculator](http://babbage.cs.qc.edu/courses/cs341/IEEE-754.html)\n",
    "\n",
    "\n",
    "- [Numerical Computing with IEEE Floating Point Arithmetic](http://epubs.siam.org/doi/book/10.1137/1.9780898718072)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Operaciones de conteo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ***Error de truncamiento:*** *¿Por qué no usar más términos en la serie de Taylor?*\n",
    "\n",
    "\n",
    "- ***Error de punto flotante:*** *¿Por qué no utilizar la mayor precisión posible?*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ejemplo 1: Multiplicación matriz - vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sea $A, B \\in \\mathbb{R}^{N \\times N}$ y $x \\in \\mathbb{R}^N$.\n",
    "\n",
    "1. Cuenta el número aproximado de operaciones que tomará para calcular $Ax$\n",
    "\n",
    "2. Hacer lo mismo para $AB$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Producto Matriz-vector:*** Definiendo $[A]_i$ como la $i$-ésima fila de $A$ y $A_{ij}$ como la $i$,$j$-ésima entrada entonces\n",
    "\n",
    "$$\n",
    "    A x = \\sum^N_{i=1} [A]_i \\cdot x = \\sum^N_{i=1} \\sum^N_{j=1} A_{ij} x_j\n",
    "$$\n",
    "\n",
    "Tomando un caso en particular, siendo $N=3$, entonces la operación de conteo es\n",
    "\n",
    "$$\n",
    "    A x = [A]_1 \\cdot v + [A]_2 \\cdot v + [A]_3 \\cdot v = \\begin{bmatrix}\n",
    "        A_{11} \\times v_1 + A_{12} \\times v_2 + A_{13} \\times v_3 \\\\\n",
    "        A_{21} \\times v_1 + A_{22} \\times v_2 + A_{23} \\times v_3 \\\\\n",
    "        A_{31} \\times v_1 + A_{32} \\times v_2 + A_{33} \\times v_3\n",
    "    \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Esto son 15 operaciones (6 sumas y 9 multiplicaciones)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tomando otro caso, siendo $N=4$, entonces el conteo de operaciones es:\n",
    "\n",
    "$$\n",
    "    A x = [A]_1 \\cdot v + [A]_2 \\cdot v + [A]_3 \\cdot v = \\begin{bmatrix}\n",
    "        A_{11} \\times v_1 + A_{12} \\times v_2 + A_{13} \\times v_3 + A_{14} \\times v_4 \\\\\n",
    "        A_{21} \\times v_1 + A_{22} \\times v_2 + A_{23} \\times v_3 + A_{24} \\times v_4 \\\\\n",
    "        A_{31} \\times v_1 + A_{32} \\times v_2 + A_{33} \\times v_3 + A_{34} \\times v_4 \\\\\n",
    "        A_{41} \\times v_1 + A_{42} \\times v_2 + A_{43} \\times v_3 + A_{44} \\times v_4 \\\\\n",
    "    \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Esto lleva a 28 operaciones (12 sumas y 16 multiplicaciones).\n",
    "\n",
    "Generalizando, hay $N^2$ mutiplicaciones y $N(N-1)$ sumas para un total de \n",
    "\n",
    "$$\n",
    "    \\text{operaciones} = N (N - 1) + N^2 = \\mathcal{O}(N^2).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Producto Matriz-Matriz ($AB$):*** Definiendo $[B]_j$ como la $j$-ésima columna de $B$ entonces\n",
    "\n",
    "$$\n",
    "    (A B)_{ij} = \\sum^N_{i=1} \\sum^N_{j=1} [A]_i \\cdot [B]_j\n",
    "$$\n",
    "\n",
    "El producto interno de dos vectores es representado por \n",
    "\n",
    "$$\n",
    "    a \\cdot b = \\sum^N_{i=1} a_i b_i\n",
    "$$\n",
    "\n",
    "conduce a $\\mathcal{O}(3N)$ operaciones. Como hay $N^2$ entradas en la matriz resultante, tendríamos $\\mathcal{O}(N^3)$ operaciones"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Existen métodos para realizar la multiplicación matriz - matriz más rápido. En la siguiente figura vemos una colección de algoritmos a lo largo del tiempo que han podido limitar el número de operaciones en ciertas circunstancias\n",
    "$$\n",
    "    \\mathcal{O}(N^\\omega)\n",
    "$$\n",
    "![matrix multiplication operation bound](./images/bound_matrix_multiply.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ejemplo 2: Método de Horner para evaluar polinomios"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dado\n",
    "\n",
    "$$P_N(x) = a_0 + a_1 x + a_2 x^2 + \\ldots + a_N x^N$$ \n",
    "\n",
    "o\n",
    "\n",
    "\n",
    "$$P_N(x) = p_1 x^N + p_2 x^{N-1} + p_3 x^{N-2} + \\ldots + p_{N+1}$$\n",
    "\n",
    "queremos encontrar la mejor vía para evaluar $P_N(x)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Primero considere dos vías para escribir $P_3$\n",
    "\n",
    "$$ P_3(x) = p_1 x^3 + p_2 x^2 + p_3 x + p_4$$\n",
    "\n",
    "y usando multiplicación anidada\n",
    "\n",
    "$$ P_3(x) = ((p_1 x + p_2) x + p_3) x + p_4$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considere cuántas operaciones se necesitan para cada...\n",
    "\n",
    "$$ P_3(x) = p_1 x^3 + p_2 x^2 + p_3 x + p_4$$\n",
    "\n",
    "$$P_3(x) = \\overbrace{p_1 \\cdot x \\cdot x \\cdot x}^3 + \\overbrace{p_2 \\cdot x \\cdot x}^2 + \\overbrace{p_3 \\cdot x}^1 + p_4$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sumando todas las operaciones, en general podemos pensar en esto como una pirámide\n",
    "\n",
    "![Original Count](./images/horners_method_big_count.png)\n",
    "\n",
    "podemos estimar de esta manera que el algoritmo escrito de esta manera tomará aproximadamente $\\mathcal{O}(N^2/2)$ operaciones para completar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mirando nuetros otros medios de evaluación\n",
    "\n",
    "$$ P_3(x) = ((p_1 x + p_2) x + p_3) x + p_4$$\n",
    "\n",
    "Aquí encontramos que el método es $\\mathcal{O}(N)$ (el 2 generalmente se ignora en estos casos). Lo importante es que la primera evaluación es $\\mathcal{O}(N^2)$ y la segunda $\\mathcal{O}(N)$!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Algoritmo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Complete la función e implemente el método de *Horner*\n",
    "\n",
    "```python\n",
    "def eval_poly(p, x):\n",
    "    \"\"\"Evaluates polynomial given coefficients p at x\n",
    "    \n",
    "    Function to evaluate a polynomial in order N operations.  The polynomial is defined as\n",
    "    \n",
    "    P(x) = p[0] x**n + p[1] x**(n-1) + ... + p[n-1] x + p[n]\n",
    "    \n",
    "    The value x should be a float.\n",
    "    \"\"\"\n",
    "    pass\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scalar version\n",
    "def eval_poly(p, x):\n",
    "    \"\"\"Evaluates polynomial given coefficients p at x\n",
    "    \n",
    "    Function to evaluate a polynomial in order N operations.  The polynomial is defined as\n",
    "    \n",
    "    P(x) = p[0] x**n + p[1] x**(n-1) + ... + p[n-1] x + p[n]\n",
    "    \n",
    "    The value x should be a float.\n",
    "    \"\"\"\n",
    "    \n",
    "    y = p[0]\n",
    "    for coefficient in p[1:]:\n",
    "        y = y * x + coefficient\n",
    "    \n",
    "    return y\n",
    "\n",
    "# Vectorized version\n",
    "def eval_poly(p, x):\n",
    "    \"\"\"Evaluates polynomial given coefficients p at x\n",
    "    \n",
    "    Function to evaluate a polynomial in order N operations.  The polynomial is defined as\n",
    "    \n",
    "    P(x) = p[0] x**n + p[1] x**(n-1) + ... + p[n-1] x + p[n]\n",
    "    \n",
    "    The value x can by a NumPy ndarray.\n",
    "    \"\"\"\n",
    "    \n",
    "    y = numpy.ones(x.shape) * p[0]\n",
    "    for coefficient in p[1:]:\n",
    "        y = y * x + coefficient\n",
    "    \n",
    "    return y\n",
    "\n",
    "p = [1, -3, 10, 4, 5, 5]\n",
    "x = numpy.linspace(-10, 10, 100)\n",
    "plt.plot(x, eval_poly(p, x))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open('./nb_style.css', 'r').read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Tabla de Contenidos",
   "title_sidebar": "Contenidos",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "165px"
   },
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}