Cap05_IntegracionNumerica.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\">Métodos Numéricos</h1>\n",
    "<h1 align=\"center\">Capítulo 5: Diferenciación e integración numérica</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/Cap05_IntegracionNumerica.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": [
    "## Diferenciación Numérica"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introducción"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La [diferenciación numérica](https://en.wikipedia.org/wiki/Numerical_differentiation) se emplea para determinar (estimar) el valor de la derivada de una función en un punto específico. No confundir con la derivada de una función, pues lo que se obtendrá es un valor puntual y no una función. En este capítulo nos centraremos únicamente en ecuiaciones unidimensionales. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Series de Taylor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "De la [serie de Taylor](https://en.wikipedia.org/wiki/Taylor_series) \n",
    "\n",
    "<a id='Ec5_1'></a>\n",
    "\\begin{equation*}\n",
    "f(x_{i \\pm 1}) = f(x_i) \\pm f'(x_i)h + \\frac{f''(x_i)h^2}{2!} \\pm \\frac{f'''(x_i)h^3}{3!} + \\ldots\n",
    "\\label{eq:Ec5_1} \\tag{5.1}\n",
    "\\end{equation*}\n",
    "\n",
    "con $h=\\Delta x = x_{i+1}-x_i$ siendo el tamaño de paso.\n",
    "\n",
    "Dada que la serie contiene infinitos términos, partir de la ecuación ($5.1$) se pueden obtener infinitos esquemas numéricos para determinar cada una de las infinitas derivadas de dicho polinomio. En este curso usaremos la técnica de [Diferencias Finitas](https://en.wikipedia.org/wiki/Finite_difference) para desarrollarlas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Esquemas de diferencias finitas para la primera derivada"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Esquema de primer orden hacia adelante (forward)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "De la ecuación [(5.1)](#Ec5_1) tomando los valores positivos, que involucran únicamente términos hacia adelante, se trunca la serie hasta la primera derivada y se realiza un despeje algebraico para llegar a:\n",
    "\n",
    "<a id='Ec5_2'></a>\n",
    "\\begin{equation*}\n",
    "f'(x_i) = \\frac{f(x_{i+1})-f(x_i)}{h} + \\mathcal{O}(h)\n",
    "\\label{eq:Ec5_2} \\tag{5.2}\n",
    "\\end{equation*}\n",
    "\n",
    "se puede observar que el término $\\mathcal{O}(h)$ indica que el error es de orden lineal, es decir, si se reduce el tamaño de paso, $h$, a la mitad, el error se reducirá a la mitad. Si se reduc el tamaño de paso a una cuarta parte, el error se reducirá, linealmente, una cuarta parte."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Esquema de primer orden hacia atrás (backward)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "De la ecuación [(5.1)](#Ec5_1) tomando los valores negativos, que involucran únicamente términos hacia atrás (backward), se trunca la serie hasta la primera derivada y se realiza un despeje algebraico para llegar a:\n",
    "\n",
    "<a id='Ec5_3'></a>\n",
    "\\begin{equation*}\n",
    "f'(x_i) = \\frac{f(x_{i})-f(x_{i-1})}{h} + \\mathcal{O}(h)\n",
    "\\label{eq:Ec5_3} \\tag{5.3}\n",
    "\\end{equation*}\n",
    "\n",
    "se observa que se llega a una expresión similar a la de la ecuación [(5.2)](#Ec5_2), pero de esta vez, se tiene en cuenta es el valor anterior al punto $x_i$. También se observa que el error es de orden lineal, por lo que se mantiene un esquema de primer orden.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Esquema de segundo orden (central)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una forma de aumentar el orden de estos esquemas, es realizar el truncamiento de la *serie de Taylor* hasta la segunda derivada, hacia adelante y hacia atras, y realizar su resta aritmética.\n",
    "\n",
    "<a id='Ec5_4'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "f(x_{i+1}) & = f(x_i) + f'(x_i)h + \\frac{f''(x_i)h^2}{2!} \\\\\n",
    "- \\\\\n",
    "f(x_{i-1}) & = f(x_i) - f'(x_i)h + \\frac{f''(x_i)h^2}{2!} \\\\\n",
    "\\hline \\\\\n",
    "f(x_{i+1}) - f(x_{i-1}) & = 2 f'(x_i)h\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_4} \\tag{5.4}\n",
    "\\end{equation*}\n",
    " \n",
    "de la anterior ecuación, despejando el término que corresponde a la primera derivada queda:\n",
    "\n",
    "<a id='Ec5_5'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "f'(x_i) = \\frac{f(x_{i+1}) - f(x_{i-1})}{2h} + \\mathcal{O}(h^2)\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_5} \\tag{5.5}\n",
    "\\end{equation*}\n",
    "\n",
    "se llega al esquema de diferencias finitas central para la primera derivada, que es de orden dos, es decir, si se disminuye el tamaño de paso, $h$, a la mitad, el error se disminuye una cuarta partes. En principio, esta es una mejor aproximación que los dos esquemas anteriores. La selección del esquema dependerá de la disponibilidad de puntos y del fenómeno físico a tratar."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Resumen esquemas diferencias finitas para la primera derivada"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como la serie de Taylor es infinita, se podrían determinar infinitos esquemas de diferentes ordenes para la primera derivada. En la siguiente tabla se presentan algunos esquemas de diferencias finitas para la primera derivada de diferentes órdenes. Se deja al estudiante la consulta de otros esquemas.\n",
    "\n",
    "|***Esquema***|***Función***|***Error***|\n",
    "|:-----:|:-----:|:---:|\n",
    "|***Forward***|$$f´(x_0)=\\frac{f(x_0+h)-f(x_0)}{h}$$|$$\\mathcal{O}(h)$$|\n",
    "|       |$$f´(x_0)=\\frac{-3f(x_0)+4f(x_0+h)-f(x_0+2h)}{2h}$$|$$\\mathcal{O}(h^2)$$|\n",
    "|***Central***|$$f´(x_0)=\\frac{f(x_0+h)-f(x_0-h)}{2h}$$|$$\\mathcal{O}(h^2)$$|\n",
    "|       |$$f´(x_0)=\\frac{f(x_0-2h)-8f(x_0-h)+8f(x_0+h)-f(x_0+2h)}{12h}$$|$$\\mathcal{O}(h^4)$$|\n",
    "|***Backward***|$$f´(x_0)=\\frac{f(x_0)-f(x_0-h)}{h}$$|$$\\mathcal{O}(h)$$|\n",
    "|        |$$f´(x_0)=\\frac{f(x_0-2h)-4f(x_0-h)+3f(x_0)}{2h}$$|$$\\mathcal{O}(h^2)$$|\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Esquemas de diferencias finitas para la segunda derivada"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Siguiendo con la misma forma de abordar el problema para la primera derivada, si se amplian los términos en la serie de Taylor hasta la tercera derivada tanto hacia adelante como hacia atrás, y se suman, se llega a:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "f(x_{i+1}) & = f(x_i) + f'(x_i)h + \\frac{f''(x_i)h^2}{2!} + \\frac{f'''(x_i)h^3}{3!}\\\\\n",
    "+ \\\\\n",
    "f(x_{i-1}) & = f(x_i) - f'(x_i)h + \\frac{f''(x_i)h^2}{2!} - \\frac{f'''(x_i)h^3}{3!}\\\\\n",
    "\\hline \\\\\n",
    "f(x_{i+1}) + f(x_{i-1}) & = 2 f(x_i) + 2f''(x_i)\\frac{h^2}{2!} + \\mathcal{O}(h^3)\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_6} \\tag{5.6}\n",
    "\\end{equation*}\n",
    " \n",
    "Despejando para el término de la segunda derivada, se llega a:\n",
    "\n",
    "<a id='Ec5_7'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "f''(x_i) = \\frac{f(x_{i+1}) - 2f(x_i) + f(x_{i-1})}{h^2} + \\mathcal{O}(h^3)\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_7} \\tag{5.7}\n",
    "\\end{equation*}\n",
    "\n",
    "Que corresponde a un esquema de diferencias finitas de segundo orden para la segunda derivada. A este esquema también se le llama \"*molécula de tres puntos*\"\n",
    "\n",
    "Igual que para la primera derivada, se pueden determinar infinitos esquemas de diferentes órdenes para la segunda derivada, y derivadas superiores. A continuación se muestra un cuadro resumen de algunos esquemas de diferencias finitas para la segunda derivada. Se deja al estudiante la revisión de esquemas de mayor orden para la segunda derivada y derivadas superiores.\n",
    "\n",
    "|***Esquema***|***Función***|***Error***|\n",
    "|:-----:|:-----:|:---:|\n",
    "|***Forward***|$$f''(x_0)=\\frac{f(x_0)-2f(x_0+h)+f(x_0+2h)}{h^2}$$|$$\\mathcal{O}(h)$$|\n",
    "|       |$$f''(x_0)=\\frac{2f(x_0)-5f(x_0+h)+4f(x_0+2h)-f(x_0+3h)}{h^2}$$|$$\\mathcal{O}(h^2)$$|\n",
    "|***Central***|$$f''(x_0)=\\frac{f(x_0-h)-2f(x_0)+f(x_0+h)}{h^2}$$|$$\\mathcal{O}(h^2)$$|\n",
    "|       |$$f''(x_0)=\\frac{-f(x_0-2h)+16f(x_0-h)-30f(x_0)+16f(x_0+h)-f(x_0+2h)}{12h^2}$$|$$\\mathcal{O}(h^4)$$|\n",
    "|***Backward***|$$f''(x_0)=\\frac{f(x_0-2h)-2f(x_0-h)+f(x_0)}{h}$$|$$\\mathcal{O}(h^2)$$|\n",
    "|        |$$f''(x_0)=\\frac{-f(x_0-3h)+4f(x_0-2h)-5f(x_0-h)+2f(x_0)}{h^2}$$|$$\\mathcal{O}(h^2)$$|\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementación computacional de algunos esquemas de diferencias finitas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A manera de ejemplo, se implementarán algunos esquemas simples de diferencias finitas para la primera derivada. Se deja como actividad a los estudiantes la implementación de otros esquemas para las diferentes derivadas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sympy as sym\n",
    "\n",
    "sym.init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Esquemas de diferencias finitas para la primera derivada\n",
    "\n",
    "def df1df(x0, h):\n",
    "    # Esquema de diferencias finitas para la primera derivada hacia adelante (forward)\n",
    "    return (f(x0 + h) - f(x0)) / h\n",
    "\n",
    "def df1db(x0, h):\n",
    "    # Esquema de diferencias finitas para la primera derivada hacia atrás (backward)\n",
    "    return (f(x0) - f(x0 - h) ) / h\n",
    "\n",
    "def df1dc(x0,h):\n",
    "    # Esquema de diferencias finitas para la primera derivada central (central)\n",
    "    return (f(x0 + h) - f(x0 - h) ) / (2 * h)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#funcion a determinar el valor de la derivada\n",
    "def f(x):\n",
    "    return 2*x**3 - 3*x**2 + 5*x+0.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#cálculo y evaluación de la primera derivada empleando cálculo simbólico\n",
    "\n",
    "def df1de(x0):\n",
    "\n",
    "    x = sym.Symbol('x')\n",
    "    df = sym.diff(f(x), x)\n",
    "    #print(df)\n",
    "    df1 = df.evalf(subs={x:x0})\n",
    "    return df1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h  = 0.05\n",
    "x0 = 0.8\n",
    "\n",
    "print(\"1st derivative \\t Value \\t\\t Error(%)\")\n",
    "print('---------------------------------------')\n",
    "\n",
    "pde = df1de(x0)\n",
    "print(\"exacta \\t\\t {0:6.4f} \\t {1}\".format(pde, '   -'))\n",
    "\n",
    "pdf = df1df(x0, h)\n",
    "epdf = abs((pde - pdf) / pde * 100) \n",
    "print(\"forward \\t {0:6.4f} \\t {1:6.2f}\".format(pdf,epdf))\n",
    "\n",
    "pdb = df1db(x0, h)\n",
    "epdb = abs((pde - pdb) / pde * 100) \n",
    "print(\"backward \\t {0:6.4f} \\t {1:6.2f}\".format(pdb,epdb))\n",
    "\n",
    "pdc = df1dc(x0,h)\n",
    "epdc = abs((pde - pdc) / pde * 100) \n",
    "print(\"central \\t {0:6.4f} \\t {1:6.2f}\".format(pdc, epdc))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integración Numérica"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introducción"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La [integración numérica](https://en.wikipedia.org/wiki/Numerical_integration) aborda una amplia gama de algoritmos para determinar el valor numérico (aproximado) de una integral definida. En este curso nos centraremos principalmente en los métodos de cuadratura, tanto de interpolación como [gaussiana](https://en.wikipedia.org/wiki/Gaussian_quadrature), como dos ejemplos de dichos algoritmos. \n",
    "\n",
    "El problema a tratar en este capítulo es la solución aproximada de la función\n",
    "\n",
    "<a id='Ec5_8'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I = \\int_a^b f(x) dx\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_8} \\tag{5.8}\n",
    "\\end{equation*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fórmulas de integración de *Newton - Cotes*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La idea básica en la integración numérica es cambiar una función difícil de integrar, $f(x)$, dada por la ecuación [(5.8)](#Ec5_8), por una función más simple, $p_n(x)$,\n",
    "\n",
    "<a id='Ec5_9'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "\\widetilde{I} \\approx \\int_{a=x_0}^{b=x_n} p_{n}(x) dx\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_9} \\tag{5.9}\n",
    "\\end{equation*}\n",
    "\n",
    "Cabe resaltar que en integración numérica no se conocerá la función a integrar, solo se dispondrá de una serie de $n+1$ puntos $(x_i, y_i), i = 0, 1, 2, \\ldots, n$, y a partir de ellos se construye un polinomio interpolante de grado $n$, $p_n$, entre los valores de los límites de integración $a = x_0$ y $b=x_n$. $p_n(x)$ es un polinomio de interpolación de la forma\n",
    "\n",
    "<a id='Ec5_10'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "p_n(x)=a_0+a_1x+a_2x^2+\\ldots+a_{n-1}x^{n-1}+a_nx^n\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_10} \\tag{5.10}\n",
    "\\end{equation*}\n",
    "\n",
    "Las fórmulas de integración de [*Newton - Cotes*](https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas), también llamadas de <a id='Quadrature'></a>[cuadratura](https://en.wikipedia.org/wiki/Quadrature_(mathematics)), son un grupo de fórmulas de integración numérica de tipo interpolación, evaluando la función en puntos equidistantes, para determinar un valor aproximado de la integral. Si no se tienen puntos espaciados, otros métodos deben ser usados, como por ejemplo cuadratura gaussiana, que se verá al final del capítulo.\n",
    "\n",
    "La forma general de las fórmulas de Newton - Cotes está dada por la función:\n",
    "\n",
    "<a id='Ec5_11'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "p_n(x)=\\sum \\limits_{i=0}^n f(x_i)L_{in}(x)\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_11} \\tag{5.11}\n",
    "\\end{equation*}\n",
    "\n",
    "donde\n",
    "\n",
    "<a id='Ec5_12'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "L_{in}(x)=\\frac{(x-x_0)\\ldots(x-x_{i-1})(x-x_{i+1})\\ldots(x-x_n)}{(x_i-x_0)\\ldots(x_i-x_{i-1})(x_i-x_{i+1})\\ldots(x_i-x_n)}\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_12} \\tag{5.12}\n",
    "\\end{equation*}\n",
    "\n",
    "es el polinomio de Lagrange, de donde se deduce que:\n",
    "\n",
    "<a id='Ec5_13'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "\\int_a^b p(x)dx=(b-a)\\sum \\limits_{i=0}^n f(x_i) \\frac{1}{(b-a)} \\int_a^b L_{in}(x)dx\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_13} \\tag{5.13}\n",
    "\\end{equation*}\n",
    "\n",
    "entonces,\n",
    "\n",
    "<a id='Ec5_14'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "\\int_a^b f(x)dx \\approx \\int_a^b p(x)dx=(b-a)\\sum \\limits_{i=0}^n w_if(x_i) \n",
    "\\end{split}\n",
    "\\label{eq:Ec5_14} \\tag{5.14}\n",
    "\\end{equation*}\n",
    "\n",
    "donde los pesos, $w_i$ de la función son representados por\n",
    "\n",
    "<a id='Ec5_15'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "w_i=\\frac{1}{(b-a)} \\int_a^b L_{in}(x)dx\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_15} \\tag{5.15}\n",
    "\\end{equation*}\n",
    "\n",
    "A partir de esta idea se obtienen los diferentes esquemas de integración numérica de *Newton - Cotes*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Regla trapezoidal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Regla trapezoidal de aplicación simple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La [regla trapezoidal](https://en.wikipedia.org/wiki/Trapezoidal_rule) emplea una aproximación de la función mediante una línea recta\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C05_Img03_TrapezoidalRule.PNG?raw=true\" width=\"250\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://upload.wikimedia.org/wikipedia/commons/4/40/Trapezoidal_rule_illustration.svg\">wikipedia.com</a> </div>\n",
    "\n",
    "y corresponde al caso en el que el polinomio en la ecuación [(5.11)](#Ec5_11) es de primer orden\n",
    "\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I=\\int_{a}^{b}f(x)dx \\approx \\int_a^b \\left[ f(a) + \\frac{f(b)-f(a)}{b-a}(x-a)\\right]dx\n",
    "= (b-a)\\frac{f(a)+f(b)}{2}\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_16} \\tag{5.16}\n",
    "\\end{equation*}\n",
    "\n",
    "Geométricamente, es equivalente a aproximar el área del trapezoide bajo la línea recta que conecta $f(a)$ y $f(b)$. La integral se representa como:\n",
    "\n",
    "$$I ≈ \\text{ancho} \\times \\text{altura promedio}$$\n",
    "\n",
    "El error en la regla trapezoidal simple se puede determinar como:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "E_t=-\\frac{1}{12}f''(\\xi)(b-a)^3\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_17} \\tag{5.17}\n",
    "\\end{equation*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Regla trapezoidal de aplicación múltiple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una manera de mejorar la exactitud de la regla trapezoidal es dividir el intervalo de integración de $a$ a $b$ en un número $n$ de segmentos y aplicar el método a cada uno de ellos. Las ecuaciones resultantes son llamadas fórmulas de integración de múltiple aplicación o compuestas.\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C05_Img04_TrapezoidalRuleMultiple.gif?raw=true\" width=\"350\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://en.wikipedia.org/wiki/Trapezoidal_rule#/media/File:Trapezium2.gif\">wikipedia.com</a> </div>\n",
    "\n",
    "Hay $n+1$ puntos base igualmente espaciados $(x_0, x_1, x_2, \\ldots, x_n)$. En consecuencia hay $n$ segmentos de igual anchura: $h = (b–a) / n$. Si $a$ y $b$ son designados como $x_0$ y $x_n$ respectivamente, la integral total se representará como:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I=\\int_{x_0}^{x_1}f(x)dx+\\int_{x_1}^{x_2}f(x)dx+\\int_{x_2}^{x_3}f(x)dx+\\ldots+\\int_{x_{n-2}}^{x_{n-1}}f(x)dx+\\int_{x_{n-1}}^{x_n}f(x)dx\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_18} \\tag{5.18}\n",
    "\\end{equation*}\n",
    "\n",
    "Al sustituir la regla trapezoidal simple en cada integrando, se tiene\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I\\approx \\left(f(x_0)+f(x_1)\\right)\\frac{h}{2}+\\left(f(x_1)+f(x_2)\\right)\\frac{h}{2}+\\left(f(x_2)+f(x_3)\\right)\\frac{h}{2}+\\ldots\\left(f(x_{n-2})+f(x_{n-1})\\right)\\frac{h}{2}+\\left(f(x_{n-1})+f(x_n)\\right)\\frac{h}{2}\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_19} \\tag{5.19}\n",
    "\\end{equation*}\n",
    "\n",
    "ahora, agrupando términos\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I\\approx \\frac{h}{2}\\left[ f(x_0) + 2\\sum_{i=1}^{n-1}f(x_i)+f(x_n) \\right]\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_20} \\tag{5.20}\n",
    "\\end{equation*}\n",
    "\n",
    "donde $h=(b-a)/n$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Implementación computacional"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trapezoidal(x):\n",
    "    n = len(x)\n",
    "    h = (x[-1] - x[0]) / n\n",
    "    \n",
    "    suma = 0\n",
    "    for i in range(1, n-1):\n",
    "        suma += funcion(x[i])\n",
    "    \n",
    "    return h * (funcion(x[0]) + 2 * suma + funcion(x[-1])) / 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def funcion(x):\n",
    "    return 4 / (1 + x**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 0\n",
    "b = 1\n",
    "n = 1\n",
    "x = np.linspace(a, b, n+1)\n",
    "I = trapezoidal(x)\n",
    "I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Error en la aplicación de la regla trapezoidal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recordando que estos esquemas provienen de la serie truncada de Taylor, el error se puede obtener determinando el primer término truncado en el esquema, que para la regla trapezoidal de aplicación simple corresponde a:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "E_t=-\\frac{1}{12}f''(\\xi)(b-a)^3\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_21} \\tag{5.21}\n",
    "\\end{equation*}\n",
    "\n",
    "donde $f''(\\xi)$ es la segunda derivada en el punto $\\xi$ en el intervalo $[a,b]$, y $\\xi$ es un valor que maximiza la evaluación de esta segunda derivada. \n",
    "\n",
    "Generalizando este concepto a la aplicación múltiple de la regla trapezoidal, se pueden sumar cada uno de los errores en cada segmento para dar:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "E_t=-\\frac{(b-a)^3}{12n^3}\\sum\\limits_{i=1}^n f''(\\xi_i)\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_22} \\tag{5.22}\n",
    "\\end{equation*}\n",
    "\n",
    "el anterior resultado se puede simplificar estimando la media, o valor promedio, de la segunda derivada para todo el intervalo\n",
    "\n",
    "<a id='Ec5_23'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "\\bar{f''} \\approx \\frac{\\sum \\limits_{i=1}^n f''(\\xi_i)}{n}\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_23} \\tag{5.23}\n",
    "\\end{equation*}\n",
    "\n",
    "de esta ecuación se tiene que $\\sum f''(\\xi_i)\\approx nf''$, y reemplazando en la ecuación [(5.23)](#Ec5_23)\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "E_t \\approx \\frac{(b-a)^3}{12n^2}\\bar{f''}\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_24} \\tag{5.24}\n",
    "\\end{equation*}\n",
    "\n",
    "De este resultado se observa que si se duplica el número de segmentos, el error de truncamiento se disminuirá a una cuarta parte."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reglas de Simpson"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Las [reglas de Simpson](https://en.wikipedia.org/wiki/Simpson%27s_rule) son esquemas de integración numérica en honor al matemático [*Thomas Simpson*](https://en.wikipedia.org/wiki/Thomas_Simpson), utilizado para obtener la aproximación de la integral empleando interpolación polinomial sustituyendo a $f(x)$. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Regla de Simpson1/3 de aplicación simple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La primera regla corresponde a una interpolación polinomial de segundo orden sustituida en la ecuación [(5.8)](#Ec5_8)\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C05_Img05_SimpsonRule13.PNG?raw=true\" width=\"350\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://upload.wikimedia.org/wikipedia/commons/c/ca/Simpsons_method_illustration.svg\">wikipedia.com</a> </div>\n",
    "\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I=\\int_a^b f(x)dx \\approx \\int_a^b p_2(x)dx\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_25} \\tag{5.25}\n",
    "\\end{equation*}\n",
    "\n",
    "del esquema de interpolación de Lagrange para un polinomio de segundo grado, visto en el capitulo anterior, y remplazando en la integral arriba, se llega a \n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I\\approx\\int_{x0}^{x2} \\left[\\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}f(x_0)+\\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}f(x_1)+\\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}f(x_2)\\right]dx\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_26} \\tag{5.26}\n",
    "\\end{equation*}\n",
    "\n",
    "realizando la integración de forma analítica y un manejo algebraico, resulta\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I\\approx\\frac{h}{3} \\left[ f(x_0)+4f(x_1)+f(x_2)\\right]\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_27} \\tag{5.27}\n",
    "\\end{equation*}\n",
    "\n",
    "donde $h=(b-a)/2$ y los $x_{i+1} = x_i + h$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A continuación, vamos a comparar graficamente las funciones \"exacta\" (con muchos puntos) y una aproximada empleando alguna técnica de interpolación para $n=3$ puntos (Polinomio interpolante de orden $2$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import barycentric_interpolate\n",
    "\n",
    "# usaremos uno de los tantos métodos de interpolación dispobibles en las bibliotecas de Python\n",
    "\n",
    "n = 3                                  # puntos a interpolar para un polinomio de grado 2\n",
    "\n",
    "xp = np.linspace(a,b,n)                # generación de n puntos igualmente espaciados para la interpolación\n",
    "fp = funcion(xp)                       # evaluación de la función en los n puntos generados\n",
    "x = np.linspace(a, b, 100)             # generación de 100 puntos igualmente espaciados\n",
    "y = barycentric_interpolate(xp, fp, x) # interpolación numérica empleando el método del Baricentro\n",
    "\n",
    "fig = plt.figure(figsize=(9, 6), dpi= 80, facecolor='w', edgecolor='k')\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "l, = plt.plot(x, y)\n",
    "plt.plot(x, funcion(x), '-', c='red')\n",
    "plt.plot(xp, fp, 'o', c=l.get_color())\n",
    "plt.annotate('Función \"Real\"', xy=(.63, 1.5), xytext=(0.8, 1.25),arrowprops=dict(facecolor='black', shrink=0.05),)\n",
    "plt.annotate('Función interpolada', xy=(.72, 1.75), xytext=(0.4, 2),arrowprops=dict(facecolor='black', shrink=0.05),)\n",
    "plt.grid(True)                          # muestra la malla de fondo\n",
    "plt.show()                              # muestra la gráfica"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se observa que hay una gran diferencia entre las áreas que se estarían abarcando en la función llamada \"*real*\" (que se emplearon $100$ puntos para su generación) y la función *interpolada* (con únicamente $3$ puntos para su generación) que será la empleada en la integración numérica (aproximada) mediante la regla de *Simpson $1/3$*.\n",
    "\n",
    "Conscientes de esto, procederemos entonces a realizar el cálculo del área bajo la curva del $p_3(x)$ empleando el método de *Simpson $1/3$*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creemos un programa en *Python* para que nos sirva para cualquier función $f(x)$ que queramos integrar en cualquier intervalo $[a,b]$ empleando la regla de integración de *Simpson $1/3$*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# se ingresan los valores del intervalo [a,b]\n",
    "a = float(input('Ingrese el valor del límite inferior: '))\n",
    "b = float(input('Ingrese el valor del límite superior: '))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cuerpo del programa por la regla de Simpson 1/3\n",
    "h = (b-a)/2   # cálculo del valor de h\n",
    "\n",
    "x0 = a        # valor del primer punto para la fórmula de S1/3 \n",
    "x1 = x0 + h   # Valor del punto intermedio en la fórmula de S1/3\n",
    "x2 = b        # valor del tercer punto para la fórmula de S1/3 \n",
    "\n",
    "fx0 = funcion(x0)   # evaluación de la función en el punto x0\n",
    "fx1 = funcion(x1)   # evaluación de la función en el punto x1\n",
    "fx2 = funcion(x2)   # evaluación de la función en el punto x2\n",
    "\n",
    "int_S13 = h / 3 * (fx0 + 4*fx1 + fx2)\n",
    "\n",
    "#erel = np.abs(exacta - int_S13) / exacta * 100\n",
    "\n",
    "print('el valor aproximado de la integral por la regla de Simpson1/3 es: ', int_S13, '\\n')\n",
    "#print('el error relativo entre el valor real y el calculado es: ', erel,'%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Error en la regla de Simpson 1/3 de aplicación simple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El problema de calcular el error de esta forma es que realmente no conocemos el valor exacto. Para poder calcular el error al usar la regla de *Simpson 1/3*:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "-\\frac{h^5}{90}f^{(4)}(\\xi)\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_28} \\tag{5.28}\n",
    "\\end{equation*}\n",
    "\n",
    "será necesario derivar cuatro veces la función original: $f(x)=e^{x^2}$. Para esto, vamos a usar nuevamente el cálculo simbólico (siempre deben verificar que la respuesta obtenida es la correcta!!!):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import *\n",
    "x = symbols('x')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Derivamos cuatro veces la función $f(x)$ con respecto a $x$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "deriv4 = diff(4 / (1 + x**2),x,4)\n",
    "deriv4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "y evaluamos esta función de la cuarta derivada en un punto $0 \\leq \\xi \\leq 1$. Como la función $f{^{(4)}}(x)$ es creciente en el intervalo $[0,1]$ (compruébelo gráficamente y/o por las técnicas vistas en cálculo diferencial), entonces, el valor que hace máxima la cuarta derivada en el intervalo dado es:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = 1.0\n",
    "evald4 = deriv4.evalf(subs={x: x0})\n",
    "print('El valor de la cuarta derivada de f en x0={0:6.2f} es {1:6.4f}: '.format(x0, evald4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculamos el error en la regla de *Simpson$1/3$*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "errorS13 = abs(h**5*evald4/90)\n",
    "print('El error al usar la regla de Simpson 1/3  es: {0:6.6f}'.format(errorS13))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Entonces, podemos expresar el valor de la integral de la función $f(x)=e^{x^2}$ en el intervalo $[0,1]$ usando la *Regla de Simpson $1/3$* como:\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "$$\\color{blue}{\\int_0^1 \\frac{4}{1 + x^2}dx} = \\color{green}{3,133333} \\color{red}{+ 0.004167}$$\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si lo fuéramos a hacer \"a mano\" $\\ldots$ aplicando la fórmula directamente, con los siguientes datos:\n",
    "\n",
    "$h = \\frac{(1.0 - 0.0)}{2.0} = 0.5$\n",
    "\n",
    "$x_0 = 0.0$\n",
    "\n",
    "$x_1 = 0.5$\n",
    "\n",
    "$x_2 = 1.0$\n",
    "\n",
    "$f(x) = \\frac{4}{1 + x^2}$\n",
    "\n",
    "sustituyendo estos valores en la fórmula dada:\n",
    "\n",
    "\n",
    "$\\int_0^1\\frac{4}{1 + x^2}dx \\approx \\frac{0.5}{3} \\left[f(0)+4f(0.5)+f(1)\\right]$\n",
    "\n",
    "$\\int_0^1\\frac{4}{1 + x^2}dx \\approx \\frac{0.5}{3} \\left[ \\frac{4}{1 + 0^2} + 4\\frac{4}{1 + 0.5^2} + \\frac{4}{1 + 1^2} \\right] \\approx 3.133333$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Regla de simpson1/3 de aplicación múltiple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Al igual que en la regla Trapezoidal, las reglas de Simpson también cuentan con un esquema de aplicación múltiple (llamada también compuesta). Supongamos que se divide el intervalo $[a,b]$ se divide en $n$ sub intervalos, con $n$ par, quedando la integral\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I=\\int_{x_0}^{x_2}f(x)dx+\\int_{x_2}^{x_4}f(x)dx+\\ldots+\\int_{x_{n-2}}^{x_n}f(x)dx\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_29} \\tag{5.29}\n",
    "\\end{equation*}\n",
    "\n",
    "y sustituyendo en cada una de ellas la regla de Simpson1/3, se llega a\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I \\approx 2h\\frac{f(x_0)+4f(x_1)+f(x_2)}{6}+2h\\frac{f(x_2)+4f(x_3)+f(x_4)}{6}+\\ldots+2h\\frac{f(x_{n-2})+4f(x_{n-1})+f(x_n)}{6}\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_30} \\tag{5.30}\n",
    "\\end{equation*}\n",
    "\n",
    "\n",
    "entonces la regla de Simpson compuesta (o de aplicación múltiple) se escribe como:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I=\\int_a^bf(x)dx\\approx \\frac{h}{3}\\left[f(x_0) + 2 \\sum \\limits_{j=1}^{n/2-1} f(x_{2j}) + 4 \\sum \\limits_{j=1}^{n/2} f(x_{2j-1})+f(x_n)\\right]\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_31} \\tag{5.31}\n",
    "\\end{equation*}\n",
    "\n",
    "donde $x_j=a+jh$ para $j=0,1,2, \\ldots, n-1, n$ con $h=(b-a)/n$, $x_0=a$ y $x_n=b$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Implementación computacional regla de Simpson1/3 de aplicación múltiple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Regla de Simpson 3/8 de aplicación simple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resulta cuando se sustituye la función $f(x)$ por una interpolación de tercer orden:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I=\\int_{a}^{b}f(x)dx = \\frac{3h}{8}\\left[ f(x_0)+3f(x_1)+3f(x_2)+f(x_3) \\right]\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_32} \\tag{5.32}\n",
    "\\end{equation*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Realizando un procedimiento similar al usado para la regla de *Simpson $1/3$*, pero esta vez empleando $n=4$ puntos:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# usaremos uno de los tantos métodos de interpolación dispobibles en las bibliotecas de Python\n",
    "\n",
    "n = 4                                  # puntos a interpolar para un polinomio de grado 2\n",
    "\n",
    "xp = np.linspace(0,1,n)                # generación de n puntos igualmente espaciados para la interpolación\n",
    "fp = funcion(xp)                       # evaluación de la función en los n puntos generados\n",
    "x = np.linspace(0, 1, 100)             # generación de 100 puntos igualmente espaciados\n",
    "y = barycentric_interpolate(xp, fp, x) # interpolación numérica empleando el método del Baricentro\n",
    "\n",
    "fig = plt.figure(figsize=(9, 6), dpi= 80, facecolor='w', edgecolor='k')\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "l, = plt.plot(x, y)\n",
    "plt.plot(x, funcion(x), '-', c='red')\n",
    "plt.plot(xp, fp, 'o', c=l.get_color())\n",
    "plt.annotate('\"Real\"', xy=(.63, 1.5), xytext=(0.8, 1.25),arrowprops=dict(facecolor='black', shrink=0.05),)\n",
    "plt.annotate('Interpolación', xy=(.72, 1.75), xytext=(0.4, 2),arrowprops=dict(facecolor='black', shrink=0.05),)\n",
    "plt.grid(True)                          # muestra la malla de fondo\n",
    "plt.show()                              # muestra la gráfica"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cuerpo del programa por la regla de Simpson 3/8\n",
    "h = (b - a) / 3   # cálculo del valor de h\n",
    "\n",
    "int_S38 = 3 * h / 8 * (funcion(a) + 3*funcion(a + h) + 3*funcion(a + 2*h) + funcion(a + 3*h) + 3 * funcion(a+4*h)/8)\n",
    "\n",
    "erel = np.abs(np.pi - int_S38) / np.pi * 100\n",
    "\n",
    "print('el valor aproximado de la integral utilizando la regla de Simpson 3/8 es: ', int_S38, '\\n')\n",
    "print('el error relativo entre el valor real y el calculado es: ', erel,'%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para poder calcular el error al usar la regla de *Simpson 3/8*:\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "$$\\color{red}{-\\frac{3h^5}{80}f^{(4)}(\\xi)}$$\n",
    "</div>\n",
    "\n",
    "será necesario derivar cuatro veces la función original. Para esto, vamos a usar nuevamente el cálculo simbólico (siempre deben verificar que la respuesta obtenida es la correcta!!!):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "errorS38 = 3*h**5*evald4/80\n",
    "print('El error al usar la regla de Simpson 3/8 es: ',errorS38)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Entonces, podemos expresar el valor de la integral de la función $f(x)=e^{x^2}$ en el intervalo $[0,1]$ usando la *Regla de Simpson $3/8$* como:\n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "$$\\color{blue}{\\int_0^1\\frac{4}{1 + x^2}dx} = \\color{green}{3.138462} \\color{red}{- 0.001852}$$\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aplicando la fórmula directamente, con los siguientes datos:\n",
    "\n",
    "$h = \\frac{(1.0 - 0.0)}{3.0} = 0.33$\n",
    "\n",
    "$x_0 = 0.0$, $x_1 = 0.33$, $x_2 = 0.66$, $x_3 = 1.00$\n",
    "\n",
    "$f(x) = \\frac{4}{1 + x^2}$\n",
    "\n",
    "sustituyendo estos valores en la fórmula dada:\n",
    "\n",
    "$\\int_0^1\\frac{4}{1 + x^2}dx \\approx \\frac{3\\times0.3333}{8} \\left[ \\frac{4}{1 + 0^2} + 3\\frac{4}{1 + 0.3333^2} +3\\frac{4}{1 + 0.6666^2} + \\frac{4}{1 + 1^2} \\right] \\approx 3.138462$\n",
    "\n",
    "\n",
    "Esta sería la respuesta si solo nos conformamos con lo que podemos hacer usando word..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Regla de Simpson3/8 de aplicación múltiple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dividiendo el intervalo $[a,b]$ en $n$ sub intervalos de longitud $h=(b-a)/n$, con $n$ múltiplo de 3, quedando la integral\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I=\\int_{x_0}^{x_3}f(x)dx+\\int_{x_3}^{x_6}f(x)dx+\\ldots+\\int_{x_{n-3}}^{x_n}f(x)dx\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_33} \\tag{5.33}\n",
    "\\end{equation*}\n",
    "\n",
    "sustituyendo en cada una de ellas la regla de Simpson3/8, se llega a\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I=\\int_a^bf(x)dx\\approx \\frac{3h}{8}\\left[f(x_0) + 3 \\sum \\limits_{i=0}^{n/3-1} f(x_{3i+1}) + 3 \\sum \\limits_{i=0}^{n/3-1}f(x_{3i+2})+2 \\sum \\limits_{i=0}^{n/3-2} f(x_{3i+3})+f(x_n)\\right]\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_34} \\tag{5.34}\n",
    "\\end{equation*}\n",
    "\n",
    "donde en cada sumatoria se deben tomar los valores de $i$ cumpliendo que $i=i+3$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Implementación computacional de la regla de Simpson3/8 de aplicación múltiple"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cuadratura de Gauss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Introducción"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Retomando la idea inicial de los esquemas de [cuadratura](#Quadrature), el valor de la integral definida se estima de la siguiente manera:\n",
    "\n",
    "<a id='Ec5_35'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I=\\int_a^b f(x)dx \\approx \\sum \\limits_{i=0}^n c_if(x_i)\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_35} \\tag{5.35}\n",
    "\\end{equation*}\n",
    "\n",
    "Hasta ahora hemos visto los métodos de la regla trapezoidal y las reglas de Simpson más empleadas. En estos esquemas, la idea central es la distribución uniforme de los puntos que siguen la regla $x_i=x_0+ih$, con $i=0,1,2, \\ldots, n$ y la evaluación de la función en estos puntos.\n",
    "\n",
    "Supongamos ahora que la restricción de la uniformidad en el espaciamiento de esos puntos fijos no es más considerada y se tiene la libertad de evaluar el área bajo una recta que conecte a dos puntos cualesquiera sobre la curva. Al ubicar estos puntos en forma “inteligente”, se puede definir una línea recta que equilibre los errores negativos y positivos\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C05_Img06_GQ01.PNG?raw=true\" width=\"750\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"http://artemisa.unicauca.edu.co/~cardila/Chapra.pdf\">Chapra, S., Canale, R. Métodos Numéricos para ingenieros, 5a Ed. Mc. Graw Hill. 2007</a> </div>\n",
    "\n",
    "De la figura de la derecha, se disponen de los puntos $x_0$ y $x_1$ para evaluar la función $f(x)$. Expresando la integral bajo la curva de forma aproximada dada en la la ecuación ([5.35](#Ec5_35)), y empleando los límites de integración en el intervalo $[-1,1]$ por simplicidad (después se generalizará el concepto a un intervalo $[a,b]$), se tiene\n",
    "\n",
    "<a id='Ec5_36'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I=\\int_{-1}^1 f(x)dx \\approx c_0f(x_0)+c_1f(x_1)\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_36} \\tag{5.36}\n",
    "\\end{equation*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Determinación de los coeficientes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "se tiene una ecuación con cuatro incógnitas ($c_0, c_1, x_0$ y $x_1$) que se deben determinar. Para ello, supongamos que disponemos de un polinomio de hasta grado 3, $f_3(x)$, de donde podemos construir cuatro ecuaciones con cuatro incógnitas de la siguiente manera:\n",
    "\n",
    "- $f_3(x)=1$:\n",
    "\n",
    "<a id='Ec5_37'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "\\int_{-1}^1 1dx = c_0 \\times 1 + c_1 \\times 1 = c_0 + c_1 = 2\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_37} \\tag{5.37}\n",
    "\\end{equation*}\n",
    "\n",
    "- $f_3(x)=x$:\n",
    "\n",
    "<a id='Ec5_38'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "\\int_{-1}^1 xdx = c_0x_0 + c_1x_1 = 0\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_38} \\tag{5.38}\n",
    "\\end{equation*}\n",
    "\n",
    "- $f_3(x)=x^2$:\n",
    "\n",
    "<a id='Ec5_39'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "\\int_{-1}^1 x^2dx = c_0x^2_0 + c_1x^2_1 = \\frac{2}{3}\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_39} \\tag{5.39}\n",
    "\\end{equation*}\n",
    "\n",
    "y por último\n",
    "\n",
    "- $f_3(x)=x^3$:\n",
    "\n",
    "<a id='Ec5_40'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "\\int_{-1}^1 x^3dx = c_0x^3_0 + c_1x^3_1 = 0\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_40} \\tag{5.40}\n",
    "\\end{equation*}\n",
    "\n",
    "resolviendo simultáneamente las dos primeras ecuaciones para $c_0$ y $c_1$ en térm,inos de $x_0$ y $x_1$, se llega a\n",
    "\n",
    "<a id='Ec5_41'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "c_0=\\frac{2x_1}{x_1-x_0}, \\quad c_1=-\\frac{2x_0}{x_1-x_0}\\end{split}\n",
    "\\label{eq:Ec5_41} \\tag{5.41}\n",
    "\\end{equation*}\n",
    "\n",
    "reemplazamos estos dos valores en las siguientes dos ecuaciones\n",
    "\n",
    "<a id='Ec5_42'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "\\frac{2}{3}=\\frac{2x_0^2x_1}{x_1-x_0}-\\frac{2x_0x_1^2}{x_1-x_0}\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_42} \\tag{5.42}\n",
    "\\end{equation*}\n",
    "\n",
    "<a id='Ec5_43'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "0=\\frac{2x_0^3x_1}{x_1-x_0}-\\frac{2x_0x_1^3}{x_1-x_0}\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_43} \\tag{5.43}\n",
    "\\end{equation*}\n",
    "\n",
    "de la ecuación ([5.43](#Ec5_43)) se tiene\n",
    "\n",
    "<a id='Ec5_44'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "x_0^3x_1&=x_0x_1^3 \\\\\n",
    "x_0^2 &= x_1^2\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_44} \\tag{5.44}\n",
    "\\end{equation*}\n",
    "\n",
    "de aquí se tiene que $|x_0|=|x_1|$ (para considerar las raíces negativas recuerde que $\\sqrt{a^2}= \\pm a = |a|$), y como se asumió que $x_0<x_1$, entonces $x_0<0$ y $x_1>0$ (trabajando en el intervalo $[-1,1]$), llegándose finalmente a que $x_0=-x_1$. Reemplazando este resultado en la ecuación ([5.42](#Ec5_42))\n",
    "\n",
    "<a id='Ec5_45'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "\\frac{2}{3}=2\\frac{x_1^3+x_1^3}{2x_1}\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_45} \\tag{5.45}\n",
    "\\end{equation*}\n",
    "\n",
    "despejando, $x_1^2=1/3$, y por último se llega a que\n",
    "\n",
    "<a id='Ec5_46'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "x_0=-\\frac{\\sqrt{3}}{3}, \\quad x_1=\\frac{\\sqrt{3}}{3}\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_46} \\tag{5.46}\n",
    "\\end{equation*}\n",
    "\n",
    "reemplazando estos resultados en la ecuación ([5.41](#Ec5_41)) y de la ecuación ([5.37](#Ec5_37)), se tiene que $c_0=c_1=1$. Reescribiendo la ecuación ([5.36](#Ec5_36)) con los valores encontrados se llega por último a:\n",
    "\n",
    "<a id='Ec5_47'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I=\\int_{-1}^1 f(x)dx &\\approx c_0f(x_0)+c_1f(x_1) \\\\\n",
    "&= f \\left( \\frac{-\\sqrt{3}}{3}\\right)+f \\left( \\frac{\\sqrt{3}}{3}\\right)\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_47} \\tag{5.47}\n",
    "\\end{equation*}\n",
    "\n",
    "\n",
    "Esta aproximación realizada es \"exacta\" para polinomios de grado menor o igual a tres ($3$). La aproximación trapezoidal es exacta solo para polinomios de grado uno ($1$).\n",
    "\n",
    "***Ejemplo:*** Calcule la integral de la función $f(x)=x^3+2x^2+1$ en el intervalo $[-1,1]$ empleando tanto las técnicas analíticas como la cuadratura de Gauss vista.\n",
    "\n",
    "\n",
    "- ***Solución analítica (exacta)***\n",
    "\n",
    "$$\\int_{-1}^1 (x^3+2x^2+1)dx=\\left.\\frac{x^4}{4}+\\frac{2x^3}{3}+x \\right |_{-1}^1=\\frac{10}{3}$$\n",
    "\n",
    "\n",
    "- ***Aproximación numérica por Cuadratura de Gauss***\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "\\int_{-1}^1 (x^3+2x^2+1)dx &\\approx1f\\left(-\\frac{\\sqrt{3}}{3} \\right)+1f\\left(\\frac{\\sqrt{3}}{3} \\right) \\\\\n",
    "&=-\\frac{3\\sqrt{3}}{27}+\\frac{2\\times 3}{9}+1+\\frac{3\\sqrt{3}}{27}+\\frac{2\\times 3}{9}+1 \\\\\n",
    "&=2+\\frac{4}{3} \\\\\n",
    "&= \\frac{10}{3}\n",
    "\\end{split}\n",
    "\\end{equation*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cambios de los límites de integración"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obsérvese que los límites de integración de la ecuación ([5.47](#Ec5_47)) son de $-1$ a $1$. Esto se hizo para simplificar las matemáticas y para hacer la formulación tan general como fuera posible. Asumamos ahora que se desea determinar el valor de la integral entre dos límites cualesquiera $a$ y $b$. Supongamos también, que una nueva variable $x_d$ se relaciona con la variable original $x$ de forma lineal,\n",
    "\n",
    "<a id='Ec5_48'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "x=a_0+a_1x_d\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_48} \\tag{5.48}\n",
    "\\end{equation*}\n",
    "\n",
    "si el límite inferior, $x=a$, corresponde a $x_d=-1$, estos valores podrán sustituirse en la ecuación ([5.48](#Ec5_48)) para obtener\n",
    "\n",
    "<a id='Ec5_49'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "a=a_0+a_1(-1)\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_49} \\tag{5.49}\n",
    "\\end{equation*}\n",
    "\n",
    "de manera similar, el límite superior, $x=b$, corresponde a $x_d=1$, para dar\n",
    "\n",
    "<a id='Ec5_50'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "b=a_0+a_1(1)\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_50} \\tag{5.50}\n",
    "\\end{equation*}\n",
    "\n",
    "resolviendo estas ecuaciones simultáneamente,\n",
    "\n",
    "<a id='Ec5_51'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "a_0=(b+a)/2, \\quad a_1=(b-a)/2\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_51} \\tag{5.51}\n",
    "\\end{equation*}\n",
    "\n",
    "sustituyendo en la ecuación ([5.48](#Ec5_48))\n",
    "\n",
    "<a id='Ec5_52'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "x=\\frac{(b+a)+(b-a)x_d}{2}\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_52} \\tag{5.52}\n",
    "\\end{equation*}\n",
    "\n",
    "derivando la ecuación ([5.52](#Ec5_52)),\n",
    "\n",
    "<a id='Ec5_53'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "dx=\\frac{b-a}{2}dx_d\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_53} \\tag{5.53}\n",
    "\\end{equation*}\n",
    "\n",
    "Las ecuacio es ([5.51](#Ec5_51)) y ([5.52](#Ec5_52)) se pueden sustituir para $x$ y $dx$, respectivamente, en la evaluación de la integral. Estas sustituciones transforman el intervalo de integración sin cambiar el valor de la integral. En este caso\n",
    "\n",
    "<a id='Ec5_54'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "\\int_a^b f(x)dx = \\frac{b-a}{2} \\int_{-1}^1 f \\left( \\frac{(b+a)+(b-a)x_d}{2}\\right)dx_d\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_54} \\tag{5.54}\n",
    "\\end{equation*}\n",
    "\n",
    "Esta integral se puede aproximar como,\n",
    "\n",
    "<a id='Ec5_55'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "\\int_a^b f(x)dx \\approx \\frac{b-a}{2} \\left[f\\left( \\frac{(b+a)+(b-a)x_0}{2}\\right)+f\\left( \\frac{(b+a)+(b-a)x_1}{2}\\right) \\right]\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_55} \\tag{5.55}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Fórmulas de punto superior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fórmula anterior para la cuadratura de Gauss era de dos puntos. Se pueden desarrollar versiones de punto superior en la forma general:\n",
    "\n",
    "<a id='Ec5_56'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "I \\approx c_0f(x_0) + c_1f(x_1) + c_2f(x_2) +\\ldots+ c_{n-1}f(x_{n-1})\n",
    "\\end{split}\n",
    "\\label{eq:Ec5_56} \\tag{5.56}\n",
    "\\end{equation*}\n",
    "\n",
    "con $n$, el número de puntos.\n",
    "\n",
    "Debido a que la cuadratura de Gauss requiere evaluaciones de la función en puntos espaciados uniformemente dentro del intervalo de integración, no es apropiada para casos donde se desconoce la función. Si se conoce la función, su ventaja es decisiva.\n",
    "\n",
    "En la siguiente tabla se presentan los valores de los parámertros para $1, 2, 3, 4$ y $5$ puntos. \n",
    "\n",
    "|$$n$$ |  $$c_i$$   |     $$x_i$$   |\n",
    "|:----:|:----------:|:-------------:|\n",
    "|$$1$$ |$$2.000000$$| $$0.000000$$  |\n",
    "|$$2$$ |$$1.000000$$|$$\\pm0.577350$$|\n",
    "|$$3$$ |$$0.555556$$|$$\\pm0.774597$$|\n",
    "|      |$$0.888889$$| $$0.000000$$  |\n",
    "|$$4$$ |$$0.347855$$|$$\\pm0.861136$$|\n",
    "|      |$$0.652145$$|$$\\pm0.339981$$|\n",
    "|$$5$$ |$$0.236927$$|$$\\pm0.906180$$|\n",
    "|      |$$0.478629$$|$$\\pm0.538469$$|\n",
    "|      |$$0.568889$$| $$0.000000$$  |\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "GaussTable = [[[0], [2]], [[-1/np.sqrt(3), 1/np.sqrt(3)], [1, 1]], [[-np.sqrt(3/5), 0, np.sqrt(3/5)], [5/9, 8/9, 5/9]], [[-0.861136, -0.339981, 0.339981, 0.861136], [0.347855, 0.652145, 0.652145, 0.347855]], [[-0.90618, -0.538469, 0, 0.538469, 0.90618], [0.236927, 0.478629, 0.568889, 0.478629, 0.236927]], [[-0.93247, -0.661209, -0.238619, 0.238619, 0.661209, 0.93247], [0.171324, 0.360762, 0.467914, 0.467914, 0.360762,  0.171324]]]\n",
    "display(pd.DataFrame(GaussTable, columns=[\"Integration Points\", \"Corresponding Weights\"]))\n",
    "def IG(f, n):\n",
    "  n = int(n)\n",
    "  return sum([GaussTable[n - 1][1][i]*f(GaussTable[n - 1][0][i]) for i in range(n)])\n",
    "def f(x): return x**9 + x**8\n",
    "IG(f, 5.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplo Cuadratura de Gauss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine el valor aproximado de:\n",
    "\n",
    "$$\\int_0^1 \\frac{4}{1+x^2}dx$$\n",
    "\n",
    "empleando cuadratura gaussiana de dos puntos.\n",
    "\n",
    "Reemplazando los parámetros requeridos en la ecuación ([5.55](#Ec5_55)), donde $a=0$, $b=1$, $x_0=-\\sqrt{3}/3$ y $x_1=\\sqrt{3}/3$\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "\\int_0^1 f(x)dx &\\approx \\frac{1-0}{2} \\left[f\\left( \\frac{(1+0)+(1-0)\\left(-\\frac{\\sqrt{3}}{3}\\right)}{2}\\right)+f\\left( \\frac{(1+0)+(1-0)\\left(\\frac{\\sqrt{3}}{3}\\right)}{2}\\right) \\right]\\\\\n",
    "&= \\frac{1}{2} \\left[f\\left( \\frac{1-\\frac{\\sqrt{3}}{3}}{2}\\right)+f\\left( \\frac{1+\\frac{\\sqrt{3}}{3}}{2}\\right) \\right]\\\\\n",
    "&= \\frac{1}{2} \\left[ \\frac{4}{1 + \\left( \\frac{1-\\frac{\\sqrt{3}}{3}}{2} \\right)^2}+\\frac{4}{1 + \\left( \\frac{1+\\frac{\\sqrt{3}}{3}}{2} \\right)^2} \\right]\\\\\n",
    "&=3.147541\n",
    "\\end{split}\n",
    "\\end{equation*}\n",
    "\n",
    "Ahora veamos una breve implementación computacional"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fxG(a, b, x):\n",
    "    xG = ((b + a) + (b - a) * x) / 2\n",
    "    return funcion(xG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def GQ2(a,b):\n",
    "    c0 = 1.0\n",
    "    c1 = 1.0\n",
    "    x0 = -1.0 / np.sqrt(3)\n",
    "    x1 =  1.0 / np.sqrt(3)\n",
    "    \n",
    "    return (b - a) / 2 * (c0 * fxG(a,b,x0) + c1 * fxG(a,b,x1))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(GQ2(0,1))"
   ]
  },
  {
   "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"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": true,
   "user_envs_cfg": true
  },
  "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": "212px"
   },
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}