Cap04_InterpolacionNumerica.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\">Métodos Numéricos</h1>\n",
    "<h1 align=\"center\">Capítulo 4: Interpolació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/Cap04_InterpolacionNumerica.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": [
    "## Introducción"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La información (datos) resultante de la medición de un evento ya sea natural o social viene dada en forma discreta o tabular, es decir, se expresa como un conjunto de pares ordenados $(x_i,y_i)$. Por ejemplo, los datos obtenidos de los censos poblacionales realizados en Colombia desde 1985 según el [DANE](https://www.dane.gov.co/) son:\n",
    "\n",
    "|Año|Población*|\n",
    "|:----:|:----:|\n",
    "|1985|30802|\n",
    "|1990|34130|\n",
    "|1995|37472|\n",
    "|2000|40296|\n",
    "|2005|42889|\n",
    "|2010|45510|\n",
    "|2015|48203|\n",
    "\n",
    "(\\* en miles de habitantes)\n",
    "\n",
    "***Nota:*** Ejemplo tomado de las notas de clase del curso Simulación Computacional de la Universidad EAFIT y es autoría  del profesor Nicolás Guarín."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.interpolate\n",
    "from scipy.optimize import curve_fit\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "start = 1985.\n",
    "stop = 2015.\n",
    "step = int((stop - start) / 5) + 1\n",
    "\n",
    "x = np.linspace(start, stop, step)\n",
    "y = [30802, 34130, 37472, 40296, 42889, 45510, 48203]\n",
    "\n",
    "plt.plot(x,y, 'o')\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si quisieramos responder a la pregunta: ¿Cuál era la población de Colombia en el 2012? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(x, y, 'o')\n",
    "ax.vlines(x=2012, ymin=30000.0, ymax=48203.0, color='r')\n",
    "ax.hlines(y=46000.0, xmin=1985.0, xmax=2012.0, color='y')\n",
    "ax.hlines(y=47000.0, xmin=1985.0, xmax=2012.0, color='b')\n",
    "ax.hlines(y=48000.0, xmin=1985.0, xmax=2012.0, color='g')\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "¿cuál de los valores mostrados arriba es el más correcto? Podemos plantear varias ideas para determinar dicho valor:\n",
    "\n",
    "- Considerando que la función es constante entre los valores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(x, y, 'o')\n",
    "\n",
    "ax.hlines(y=y[0], xmin=x[0], xmax=x[0] + 2, color='b')\n",
    "ax.vlines(x=x[0] + 2, ymin=y[0], ymax=y[1], color='b')\n",
    "\n",
    "for i in range(1, len(x)-1):\n",
    "    ax.hlines(y=y[i], xmin=x[i]-3, xmax= x[i] + 2, color='b')\n",
    "    ax.vlines(x=x[i]+2, ymin=y[i], ymax=y[i+1], color='b')\n",
    "\n",
    "ax.hlines(y=y[-1], xmin=x[-1]-3, xmax=x[-1]+2, color='b')\n",
    "\n",
    "ax.vlines(x = 2012, ymin = 30000.0, ymax = y[-1], color = 'r')\n",
    "ax.hlines(y = y[-1], xmin = 1985.0, xmax = 2012.0, color = 'r')\n",
    "\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Asumiendo que la función es lineal entre valores "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(x, y, 'o--')\n",
    "ax.vlines(x = 2012, ymin = 30000.0, ymax = 46700, color = 'r')\n",
    "ax.hlines(y = 46700, xmin = 1985.0, xmax = 2012.0, color = 'r')\n",
    "\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Determinando un polinomio que pase por cada uno de los puntos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(x, y, 'o')\n",
    "\n",
    "t = np.linspace(0, 1, len(x)) # parameter t to parametrize x and y\n",
    "pxLagrange = scipy.interpolate.lagrange(t, x) # X(T)\n",
    "pyLagrange = scipy.interpolate.lagrange(t, y) # Y(T)\n",
    "n = 100\n",
    "ts = np.linspace(t[0],t[-1],n)\n",
    "xLagrange = pxLagrange(ts) # lagrange x coordinates\n",
    "yLagrange = pyLagrange(ts) # lagrange y coordinates\n",
    "ax.plot(xLagrange, yLagrange,'b-')\n",
    "\n",
    "ax.vlines(x=2012, ymin=30000.0, ymax=46700, color='r')\n",
    "ax.hlines(y=46700, xmin=1985.0, xmax=2012.0, color='r')\n",
    "\n",
    "ax.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ajustando la curva que mejor se aproxime a cada uno de los datos. En este ejemplo haremos un ajuste lineal, pero no es la única forma de hacerlo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the true objective function\n",
    "def objective(x, a, b):\n",
    "    return a * x + b\n",
    " \n",
    "# curve fit\n",
    "popt, _ = curve_fit(objective, x, y)\n",
    "\n",
    "# summarize the parameter values\n",
    "a, b = popt\n",
    "print('y = %.5f * x + %.5f' % (a, b))\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "# plot input vs output\n",
    "ax.scatter(x, y)\n",
    "\n",
    "# define a sequence of inputs between the smallest and largest known inputs\n",
    "x_line = np.arange(min(x), max(x), 1)\n",
    "\n",
    "# calculate the output for the range\n",
    "y_line = objective(x_line, a, b)\n",
    "\n",
    "# create a line plot for the mapping function\n",
    "ax.plot(x_line, y_line, '--', color='b')\n",
    "\n",
    "yfit = a * 2012 + b\n",
    "\n",
    "ax.vlines(x=2012, ymin=30000.0, ymax=yfit, color='r')\n",
    "ax.hlines(y=yfit, xmin=1985.0, xmax=2012.0, color='r')\n",
    "\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En este curso nos centraremos en los esquemas de interpolación, por lo que los esquemas de ajuste no se tratarán."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Propósitos de la interpolación"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Los problemas de interpolación surgen de muchas fuentes diferentes y pueden tener muchos propósitos diferentes. Algunos de estos incluyen:\n",
    "\n",
    "- Trazar una curva suave a través de puntos de datos discretos\n",
    "\n",
    "\n",
    "- Evaluación rápida y sencilla de una función matemática\n",
    "\n",
    "\n",
    "- Reemplazar una función difícil por una fácil\n",
    "\n",
    "\n",
    "- Diferenciar o integrar datos tabulares"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Diferencias entre Interpolación, Aproximación y Ajuste de curvas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Las técnicas para resolver el problema de determinar un valor intermedio entre dos valores conocidos se pueden enmarcar en:\n",
    "\n",
    "- [Interpolación](https://en.wikipedia.org/wiki/Interpolation)\n",
    "\n",
    "\n",
    "- [Aproximación](https://en.wikipedia.org/wiki/Approximation_theory)\n",
    "\n",
    "\n",
    "- [Ajuste de curvas](https://en.wikipedia.org/wiki/Curve_fitting)\n",
    "\n",
    "A continuación vamos a describir brevemente las diferencias entre ellas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpolación vs Aproximacion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En interpolación, se ajustan todos los puntos de datos de forma exacta mientras que la aproximación, como su nombre indica, solo se aproxima.\n",
    "\n",
    "Cuando se trata de la idoneidad, la interpolación es apropiada para suavizar esos datos ruidosos y no es apropiada cuando los puntos de datos están sujetos a errores experimentales u otras fuentes de error significativo. Tener un gran conjunto de puntos de datos también puede sobrecargar la interpolación. Por otro lado, la aproximación es principalmente apropiada para el diseño de rutinas de biblioteca para calcular funciones especiales. Esto se debe a la naturaleza de estas funciones: se considera que los valores exactos no son esenciales y, hasta cierto punto, ineficaces cuando los valores aproximados funcionan."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpolación vs Ajuste de curvas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En el ajuste de curvas, no ajustamos todos nuestros puntos de datos. Por eso tenemos el concepto de residuos. En la interpolación, se obliga a la función a ajustarse a todos los puntos de datos. Ver la referencia del [Cuarteto de ascombe](https://es.wikipedia.org/wiki/Cuarteto_de_Anscombe) como un contra ejemplo de los inconvenientes en los ajustes de curvas.\n",
    "\n",
    "Ahora que sabemos de qué categoría estamos hablando, reduzcamos a las familias de funciones utilizadas para la interpolación."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Elección de la función de interpolación"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Es importante darse cuenta de que existe cierta arbitrariedad en la mayoría de los problemas de interpolación. Hay arbitrariamente muchas funciones que interpolan un conjunto dado de datos. Simplemente requiriendo que alguna función matemática se ajuste a los puntos de datos deje exactamente abiertos tales\n",
    "preguntas como:\n",
    "\n",
    "- ¿Qué forma debe tener la función? Puede haber consideraciones matemáticas o físicas relevantes que sugieran una forma particular de interpolante.\n",
    "\n",
    "\n",
    "- ¿Cómo debería comportarse la función entre puntos de datos?\n",
    "\n",
    "\n",
    "- ¿Debería la función heredar propiedades de los datos, como monotonicidad, convexidad o periodicidad?\n",
    "\n",
    "\n",
    "- Si se grafican la función y los datos, ¿los resultados deberían ser agradables a la vista?\n",
    "\n",
    "\n",
    "- ¿Estamos interesados principalmente en los valores de los parámetros que definen la función de interpolación, o simplemente en evaluar la función en varios puntos para graficar u otros propósitos?\n",
    "\n",
    "\n",
    "La elección de la función de interpolación depende de las respuestas a estas preguntas, así como de los datos que se deben ajustar y generalmente se basa en:\n",
    "\n",
    "- Qué tan fácil es trabajar con la función (determinar sus parámetros a partir de los datos, evaluar la función en un punto dado, diferenciar o integrar la función, etc.)\n",
    "\n",
    "\n",
    "- Qué tan bien las propiedades de la función coinciden con las propiedades de los datos a ser t (suavidad, monotonicidad, convexidad, periodicidad, etc.)\n",
    "\n",
    "\n",
    "Algunas familias de funciones que se utilizan comúnmente para la interpolación incluyen:\n",
    "\n",
    "\n",
    "- [Polinomios](https://en.wikipedia.org/wiki/Polynomial_interpolation)\n",
    "\n",
    "\n",
    "- [Interpolación a trazos](https://en.wikipedia.org/wiki/Spline_interpolation)\n",
    "\n",
    "\n",
    "- Funciones trigonométricas\n",
    "\n",
    "\n",
    "- Exponenciales\n",
    "\n",
    "\n",
    "- Funciones racionales\n",
    "\n",
    "\n",
    "En este capítulo nos centraremos en la interpolación por polinomios e interpolación a trazos.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolación Polinomial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introducción"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La interpolación polinomial es el tipo de interpolación más simple y común. Una de sus características es que siempre hay un polinomio único de grado como máximo $n-1$ que pasa por $n$ puntos de datos.\n",
    "Hay muchas formas de calcular o representar un polinomio, pero se reducen a la misma función matemática. Algunos de los métodos son la base monomial, la base de [Lagrange](https://en.wikipedia.org/wiki/Lagrange_polynomial) y la base de [Newton](https://en.wikipedia.org/wiki/Newton_polynomial). Como puede observar, reciben el nombre de su base.\n",
    "\n",
    "***Inconvenientes:***\n",
    "\n",
    "- ***Polinomio de alto grado:*** una elección adecuada de funciones de base y puntos de interpolación puede mitigar algunas de las dificultades asociadas con polinomio de alto grado.\n",
    "\n",
    "\n",
    "- ***Sobreajuste ([overfitting](https://en.wikipedia.org/wiki/Overfitting)):*** ajuste de un solo polinomio a una gran cantidad de puntos de datos, lo que probablemente produciría un comportamiento de oscilación insatisfactorio en el interpolante."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La fórmula general de un polinomio de $n$-ésimo orden es:\n",
    "\n",
    "\\begin{equation*}\n",
    "f_n(x) = a_0 + a_1x + a_2x^2 +…+  a_nx^n\n",
    "\\label{eq:Ec4_1} \\tag{4.1}\n",
    "\\end{equation*}\n",
    "\n",
    "El polinomio de interpolación dado por la ecuación $\\eqref{eq:Ec4_1}$ consiste en determinar el único polinomio de n-ésimo orden que se ajusta a los $n+1$ puntos dados. Este polinomio proporciona una fórmula para calcular los valores intermedios."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpolación Lineal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El método más simple de interpolación es conectar dos puntos mediante una línea recta.\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C04_Img01_InterpolLineal.PNG?raw=true\" width=\"250\" />\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 se tiene, por semejanza de triángulos:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\frac{f_1(x)-f(x_0)}{x-x_0}=\\frac{f(x_1)-f(x_0)}{x_1-x_0}\n",
    "\\label{eq:Ec4_2} \\tag{4.2}\n",
    "\\end{equation*}\n",
    "\n",
    "reordenando,\n",
    "\n",
    "\\begin{equation*}\n",
    "f_1(x)=f(x_0)+\\frac{f(x_1)-f(x_0)}{x_1-x_0}(x-x_0)\n",
    "\\label{eq:Ec4_3} \\tag{4.3}\n",
    "\\end{equation*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplo interpolación lineal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Estime el valor de $Ln(2)$ empleando la interpolación lineal, entre $x_0=1$ y $x_1=6$\n",
    "\n",
    "Evaluando el valor del logaritmo en cada uno de los dos puntos, $Ln(1)=0$ y $Ln(6)=1.791759$\n",
    "\n",
    "$$f_1(2)=0+\\frac{1.791759-0}{6-1}(2-1)=0.3583519$$\n",
    "\n",
    "El valor exacto es $Ln(2)=0.693147$, que representa un error relativo porcentual de\n",
    "\n",
    "$$Er(\\%)=\\frac{|0.693147-0.3583519|}{0.693147}=48.3\\%$$\n",
    "\n",
    "Si se disminuye el valor del intervalo a evaluar, por ejemplo en $x_1=4$, se llega a \n",
    "\n",
    "$$f_1(2)=0+\\frac{1.386294-0}{4-1}(2-1)=0.462098$$\n",
    "\n",
    "obteniendo un error relativo porcentual del $33.3\\%$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visualización computacional"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "n = 20\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "x = np.linspace(1, 6, n)\n",
    "y = np.log(x)\n",
    "\n",
    "x1 = [x[0], x[-1]]\n",
    "y1 = [np.log(x[0]), np.log(x[-1])]\n",
    "\n",
    "x2 = [x[0], 4]\n",
    "y2 = [np.log(x[0]), np.log(4)]\n",
    "\n",
    "ax.plot(x,y, '-', x1, y1, 'o-', x2,y2, 'o-')\n",
    "ax.vlines(x = 2, ymin = 0.0, ymax = np.log(2), color = 'r', linestyles='dashed')\n",
    "\n",
    "plt.grid(True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El error en la interpolación lineal resulta de aproximar una curva con una línea recta.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Mejoras al esquema de interpolación lineal"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Disminuir el tamaño del intervalo.\n",
    "\n",
    "\n",
    "- Introducir alguna curvatura en la línea que conecta los puntos."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='DDN'></a>\n",
    "### Interpolación cuadrática"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Si tres (3) puntos de los datos están disponibles, esto puede realizarse con un polinomio de segundo grado (parábola).\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C04_Img02_InterpolCuadratica.PNG?raw=true\" width=\"350\" />\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",
    "\n",
    "La forma general de una ecuación cuadrática puede expresarse de la siguiente forma:\n",
    "\n",
    "\\begin{equation*}\n",
    "f_2(x) = b_0 + b_1(x–x_0) + b_2(x–x_0)(x–x_1)\n",
    "\\label{eq:Ec4_4} \\tag{4.4}\n",
    "\\end{equation*}\n",
    "\n",
    "Debemos determinar los valores de los coeficientes $b_i$.\n",
    "\n",
    "- Para $b_0$, en la ecuación $\\eqref{eq:Ec4_4}$, con $x = x_0$:\n",
    "\n",
    "\\begin{equation*}\n",
    "b_0 = f(x_0)\n",
    "\\label{eq:Ec4_5} \\tag{4.5}\n",
    "\\end{equation*}\n",
    "\n",
    "- Para $b_1$, sustituyendo la ecuación $\\eqref{eq:Ec4_5}$ en la la ecuación $\\eqref{eq:Ec4_4}$, y evaluando en $x = x_1$:\n",
    "\n",
    "\\begin{equation*}\n",
    "b_1 = \\frac{f(x_1)-f(x_0)}{(x_1 - x_0)}\n",
    "\\label{eq:Ec4_6} \\tag{4.6}\n",
    "\\end{equation*}\n",
    "\n",
    "- Para $b_2$, las ecuaciones $\\eqref{eq:Ec4_5}$ y $\\eqref{eq:Ec4_6}$ pueden sustituirse en la ecuación $\\eqref{eq:Ec4_4}$, evaluada en $x_2$\n",
    "\n",
    "\\begin{equation*}\n",
    "b_2=\\frac{\\frac{f(x_2)-f(x_1)}{(x_2-x_1)}-\\frac{f(x_1)-f(x_0)}{(x_1-x_0)}}{(x_2 - x_0)}\n",
    "\\label{eq:Ec4_7} \\tag{4.7}\n",
    "\\end{equation*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplo interpolación cuadrática"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Continuando con el ejemplo anterior, se van a considerar los siguientes puntos:\n",
    "\n",
    "$$x_0=1 \\hspace{1cm} f(x_0)=0.000000$$\n",
    "$$x_1=4 \\hspace{1cm} f(x_1)=1.386294$$\n",
    "$$x_2=6 \\hspace{1cm} f(x_2)=1.791759$$\n",
    "\n",
    "De las ecuaciones anteriores, \n",
    "\n",
    "$$b_0=0$$\n",
    "\n",
    "$$b_1=\\frac{1.386294-0}{4-1}=0.4620981$$\n",
    "\n",
    "$$b_2=\\frac{\\frac{1.791759-1.386294}{6-4}-0.4620981}{6-1}=-0.0518731$$\n",
    "\n",
    "Sustituyendo estos valores en la ecuación cuadrática inicial, se llega a:\n",
    "\n",
    "$$f_2(x)=0+0.4620981(x-1)-0.0518731(x-1)(x-4)$$\n",
    "\n",
    "y evaluando en $x=2$, se llega a\n",
    "\n",
    "$$f_2(2)=0.565844$$\n",
    "\n",
    "que representa un error relativo porcentual del $18.4\\%$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Implementación computacional"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def difdiv2o(x, y, xm):\n",
    "    b0 = y[0]\n",
    "    b1 = (y[1] - b0) / (x[1] - x[0])\n",
    "    b2 = ((y[2] - y[1]) / (x[2] - x[1]) - b1) / (x[2] - x[0])\n",
    "\n",
    "    return b0 + b1 * (xm - x[0]) + b2 * (xm - x[0]) * (xm - x[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "n = 20\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "x = np.linspace(1, 6, n)\n",
    "y = np.log(x)\n",
    "\n",
    "x1 = [x[0], 4, x[-1]]\n",
    "y1 = [np.log(x[0]), np.log(x1[1]), np.log(x[-1])]\n",
    "ym = difdiv2o(x1,y1,2)\n",
    "\n",
    "x1.insert(1, 2)\n",
    "y1.insert(1,ym)\n",
    "\n",
    "ax.plot(x, y, '-', x1, y1, 'o-')\n",
    "ax.vlines(x = 2, ymin = 0.0, ymax = np.log(2), color = 'r', linestyles='dashed')\n",
    "\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Polinomio de diferencias divididas de Newton"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lo anterior puede ser generalizado para ajustar un polinomio de $n$-ésimo orden a $n+1$ datos:\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C04_Img03a_DivDifNewton?raw=true\" width=\"500\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://medium.com/\">medium.com</a> </div>\n",
    "\n",
    "\\begin{equation*}\n",
    "f_n(x) = b_0+b_1(x–x_0)+\\ldots+b_n(x–x_0)(x–x_1)\\ldots(x – x_{n-1})\n",
    "\\label{eq:Ec4_8} \\tag{4.8}\n",
    "\\end{equation*}\n",
    "\n",
    "De igual manera que para las interpolaciones lineal y cuadrática, se llega a:\n",
    "\t\n",
    "\\begin{equation*}\n",
    "f_n(x)=f(x_0)+(x–x_0)f[x_1,x_0]+(x–x_0)(x–x_1)f[x_2,x-,x_0]+\\ldots+(x–x_0)(x–x_1)\\ldots(x–x_{n-1})f[x_n, x_{n-1},\\ldots,x_2,x_1,x_0]\n",
    "\\label{eq:Ec4_9} \\tag{4.9}\n",
    "\\end{equation*}\n",
    "\n",
    "Conocido como *Polinomio de interpolación por [diferencias divididas de Newton](https://en.wikipedia.org/wiki/Divided_differences)*. Las evaluaciones de las funciones puestas entre paréntesis son diferencias divididas finitas.\n",
    " \n",
    "- ***Primera diferencia dividida:***\n",
    "\n",
    "\\begin{equation*}\n",
    "f[x_i, x_j]=\\frac{f(x_i)-f(x_j)}{(x_i-x_j)}\n",
    "\\label{eq:Ec4_10} \\tag{4.10}\n",
    "\\end{equation*}\n",
    "\n",
    "- ***Segunda diferencia dividida:*** representa la diferencia de las dos primeras diferencias divididas\n",
    "\n",
    "\\begin{equation*}\n",
    "f[x_i, x_j,x_k]=\\frac{f[x_i,x_j]-f[x_j,x_k]}{(x_i-x_k)}\n",
    "\\label{eq:Ec4_11} \\tag{4.11}\n",
    "\\end{equation*}\n",
    "\n",
    "$$\\vdots$$\n",
    "\n",
    "- ***$n$-ésima diferencia dividida:*** representa la diferencia de las dos primeras diferencias divididas\n",
    "\\begin{equation*}\n",
    "f[x_n, x_{n-1},\\ldots, x_1,x_0]=\\frac{f[x_n,x_{n-1},\\ldots,x_1]-f[x_{n-1},x_{n-2},\\ldots,x_0]}{(x_n-x_0)}\n",
    "\\label{eq:Ec4_12} \\tag{4.12}\n",
    "\\end{equation*}\n",
    "\n",
    "Este proceso recursivo lo podemos visualizar de la siguiente manera:\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C04_Img03_DivDifNewton.PNG?raw=true\" width=\"500\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://upload.wikimedia.org/wikipedia/commons/0/00/Newton_method.png\">Wikimedia.org</a> </div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Implementación computacional"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En la página 513 del libro de Chapra y Canale, Figura 18.7, se tiene un algoritmo en Fortran para la implementación del código de Diferencias Divididas tipo Newton. Se invita al estudiante a que lo estudie y codifique en el lenguaje de preferencia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Escriba aquí su código"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Análisis de Error para la interpolación polinomial tipo Newton"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La ecuación $\\eqref{eq:Ec4_9}$ es similar a la *serie de expansión de Taylor*. Se agregan términos en forma secuencial para capturar el comportamiento de alto orden de la función a analizar. Estos términos son diferencias divididas finitas y, así, representan aproximaciones de derivadas de orden mayor.\n",
    "\n",
    "El error de truncamiento se expresa entonces como:\n",
    "\n",
    "\\begin{equation*}\n",
    "R_n=\\frac{f^{(n+1)}(\\xi)}{(n+1)!} \\left ( x_{i+1}-x_i\\right )^{n+1}\n",
    "\\label{eq:Ec4_13} \\tag{4.13}\n",
    "\\end{equation*}\n",
    "\n",
    "Para una interpolación de n-ésimo orden, una relación análoga para el error es\n",
    "\n",
    "<a id='Ec4_14'></a>\n",
    "\\begin{equation*}\n",
    "R_n=\\frac{f^{(n+1)}(\\xi)}{(n+1)!}(x-x_0)(x-x_1) \\ldots (x-x_n)\n",
    "\\label{eq:Ec4_14} \\tag{4.14}\n",
    "\\end{equation*}\n",
    "\n",
    "Observe que en la ecuación [(4.14)](#Ec4_14), la función debe conocerse. Para resolver esta situación, una formulación alternativa es el uso de la diferencia dividida para aproximar la derivada $(n+1)$–ésima y que no requiere el conocimiento previo de la función.\n",
    "\n",
    "<a id='Ec4_15'></a>\n",
    "\\begin{equation*}\n",
    "R_n=f_n[x_n,x_{n-1},x_{n-2},\\ldots,x_2,x_1,x_0](x-x_0)(x-x_1) \\ldots (x-x_n)\n",
    "\\label{eq:Ec4_15} \\tag{4.15}\n",
    "\\end{equation*}\n",
    "\n",
    "Debido a que la ecuación [(4.15)](#Ec4_15) contiene el término $f_n(x)$ no puede resolverse para estimar el error, pero, si se dispone de un dato adicional, la ecuación [(4.15)](#Ec4_15) puede usarse:\n",
    "\n",
    "\\begin{equation*}\n",
    "R_n=f_n[x_{n+1}, x_n,x_{n-1},x_{n-2},\\ldots,x_2,x_1,x_0](x-x_0)(x-x_1) \\ldots (x-x_n)\n",
    "\\label{eq:Ec4_16} \\tag{4.16}\n",
    "\\end{equation*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Polinomios de Interpolación de Lagrange"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Introducción"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El [polinomio de interpolación de Lagrange](https://en.wikipedia.org/wiki/Lagrange_polynomial) evita el cálculo de las diferencias divididas en el esquema de Newton. De una forma general, se representa como la [combinación lineal](https://en.wikipedia.org/wiki/Linear_combination):\n",
    "\n",
    "\\begin{equation*}\n",
    "f_n(x)=\\sum \\limits_{i=0}^n L_i(x)f(x_i)\n",
    "\\label{eq:Ec4_17} \\tag{4.17}\n",
    "\\end{equation*}\n",
    "\n",
    "donde $L_i$ son las bases polinómicas de *[Lagrange](https://es.wikipedia.org/wiki/Joseph-Louis_Lagrange)* dadas por: \n",
    "\n",
    "\\begin{equation*}\n",
    "L_i(x)=\\prod_{\\substack{j=0\\\\ j \\ne i}}^n \\frac{x-x_j}{x_i-x_j}\n",
    "\\label{eq:Ec4_18} \\tag{4.18}\n",
    "\\end{equation*}\n",
    "\n",
    "de $\\eqref{eq:Ec4_18}$ se observa que todas las funciones $L_i$ son polinomios de grado $n$ que tienen la propiedad\n",
    "\n",
    "\\begin{equation*}\n",
    "L_i(x_j)=\\delta_{ij}, \\quad \\delta_{ij} = \\left \\{\n",
    "\\begin{aligned}\n",
    "1, \\quad i=j,\\\\\n",
    "0, \\quad i \\ne j\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "\\label{eq:Ec4_19} \\tag{4.19}\n",
    "\\end{equation*}\n",
    "\n",
    "donde $\\delta_{is}$ es el [delta de Kronecker](https://en.wikipedia.org/wiki/Kronecker_delta).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Polinomio de Interpolación de Lagrange de primer grado"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tomando $n=1$ (lineal), se tiene\n",
    "\n",
    "\\begin{equation*}\n",
    "f_1(x)=\\frac{(x-x_1)}{(x_0-x_1)}f(x_0)+\\frac{(x-x_0)}{(x_1-x_0)}f(x_1)\n",
    "\\label{eq:Ec4_19a} \\tag{4.19a}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Polinomio de Interpolación de Lagrange de segundo grado"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tomando $n=2$ (cuadrático), se tiene\n",
    "\n",
    "\\begin{equation*}\n",
    "f_2(x)=\\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)\n",
    "\\label{eq:Ec4_20} \\tag{4.20}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplo polinomio de interpolación de Lagrange"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "retomando el ejemplo que se ha venido trabajando anteriormente,\n",
    "\n",
    "$$x_0=1 \\hspace{1cm} f(x_0)=0.000000$$\n",
    "$$x_1=4 \\hspace{1cm} f(x_1)=1.386294$$\n",
    "$$x_2=6 \\hspace{1cm} f(x_2)=1.791759$$\n",
    "\n",
    "- ***Polinomio de primer grado:***\n",
    "\n",
    "$$f_1(2)=\\frac{2-4}{1-4}(0)+\\frac{2-1}{4-1}(1.386294) = 0.462098$$\n",
    "\n",
    "- ***Polinomio de segundo grado:***\n",
    "\n",
    "$$f_2(2)=\\frac{(2-4)(2-6)}{(1-4)(1-6)}(0)+\\frac{(2-1)(2-6)}{(4-1)(4-6)}(1.386294)+\\frac{(2-1)(2-4)}{(6-1)(6-4)}(1.791759) = 0.565844$$\n",
    "\n",
    "- ***Nota:*** Realice una comparación con los resultados obtenidos con los correspondientes esquemas lineal y cuadrático en el esquema de [diferencias divididas de Newton](#DDN)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Implementación computacional"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lagrange(x ,i , xm ):\n",
    "    n = len(xm) - 1\n",
    "    y = 1.0\n",
    "    for j in range(n + 1):\n",
    "        if i != j:\n",
    "            y *= (x - xm[j]) / (xm[i] - xm[j])\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def interpolation(x, xm, ym):\n",
    "    n = len(xm) - 1\n",
    "    lagrpoly = np.array([lagrange(x, i, xm) for i in range(n + 1)])\n",
    "    y = np.dot(ym, lagrpoly)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "xm = np.array([1,4,6])\n",
    "ym = np.log(xm)\n",
    "#xm = np.array([1, 2, 3, 4, 5, 6])\n",
    "#ym = np.array([-3, 0, -1, 2, 1, 4])\n",
    "#ym = np.sin(xm)\n",
    "\n",
    "xplot = np.linspace(-1., 6.0, 100)\n",
    "yplot = interpolation(xplot, xm, ym)\n",
    "plt.plot(xm, ym, '--', xplot, yplot, '-')\n",
    "plt.grid('True')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Funciones de base"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para entender un poco más cómo es que trabaja la interpolación entre cada uno de los puntos, recordemos que dichos polinomios interpolantes de Lagrange deben cumplir con la propiedad del delta de Kronecker. Para visualizar esto, se implementará la descripción propuesta en el dcumento [Interpolación de Lagrange 1D](https://github.com/AppliedMechanics-EAFIT/modelacion_computacional/blob/master/notebooks/02a_interpolacion.ipynb) realizado por los profesores *Juan David Gómez Cataño* y *Nicolás Guarín Zapata* para el curso de Modelación Computacional en el programa de Ingeniería Civil de la Universidad EAFIT. Todo el crédito para ellos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# llamado a las biblitecas numéricas, de visualización y simbólicas\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import interpolate\n",
    "import sympy as sym\n",
    "sym.init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lagrange_poly(x, order, i, xi=None): \n",
    "    if xi == None:\n",
    "        xi = sym.symbols('x:%d'%(order+1))\n",
    "    index = list(range(order+1))\n",
    "    index.pop(i)\n",
    "    return sym.prod([(x - xi[j])/(xi[i] - xi[j]) for j in index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fun = lambda x: x**3 + 4.0*x**2 - 10.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npts = 200\n",
    "x_pts = np.linspace(-1, 1, npts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pts = np.array([-1, 1, 0])\n",
    "fd = fun(pts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "y_pts = fun(x_pts)\n",
    "plt.plot(x_pts , y_pts)\n",
    "plt.plot(pts, fd, 'ko')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = sym.symbols('x')                                       \n",
    "pol = []                                                \n",
    "pol.append(sym.simplify(lagrange_poly(x, 2, 0, [-1,1,0])))  \n",
    "pol.append(sym.simplify(lagrange_poly(x, 2, 1, [-1,1,0])))\n",
    "pol.append(sym.simplify(lagrange_poly(x, 2, 2, [-1,1,0])))\n",
    "pol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "for k in range(3):\n",
    "    for i in range(npts):\n",
    "        y_pts[i] = pol[k].subs([(x, x_pts[i])])\n",
    "    plt.plot(x_pts, y_pts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Dificultades de los polinomios de Lagrange"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Los polinomios de interpolación de Lagrange presentan dificultades cuando se tienen polinomios de orden muy alto, agravado cuando se tienen puntos equidistantes o se presentan saltos (discontinuidades) en la solución.\n",
    "\n",
    "A esta situación se le conoce como *[fenómeno de Runge](https://en.wikipedia.org/wiki/Runge%27s_phenomenon)*. Veamos el siguiente ejemplo:\n",
    "\n",
    "- Dada la ecuación\n",
    "\n",
    "$$f(x)=\\frac{1}{1+25x^2}$$\n",
    "\n",
    "si se interpola esta función utilizando nodos equidistantes $x_i \\in [-1, 1]$ tal que\n",
    "\n",
    "$$x_i=-1+(i-1)\\frac{2}{n} \\quad i \\in \\{1, 2, 3, \\ldots, n, n+1\\}$$"
   ]
  },
  {
   "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": [
    "n = 200\n",
    "x = [-1 + 2 * (i - 1) / n for i in range(1, n + 2)]\n",
    "y = [1 / (1+ 25 * i **2) for i in x]\n",
    "\n",
    "plt.plot(x,y)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora vamos a emplear, como ejemplo, una serie de polinomios interpolantes para evaluar su comportamiento, desde orden 1 (2 puntos) hasta de orden 9 (10 puntos)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Polinomios interpolantes de orden 1-9 (2-10 puntos)\n",
    "\n",
    "#data = [2, 3, 4, 6, 10]\n",
    "data = [12]\n",
    "#fig = plt.figure(figsize=(16, 6), dpi= 50, facecolor='w', edgecolor='k')\n",
    "#ax = fig.add_subplot(111)\n",
    "\n",
    "plt.plot(x, y, '-', label = 'Exacta')\n",
    "plt.title('Fenómeno de Runge')\n",
    "xplot = np.linspace(-1., 1.0, 100)\n",
    "\n",
    "\n",
    "for i in data:\n",
    "    xRi = np.linspace(-1, 1, i)\n",
    "    yRi = [1 / (1+ 25 * j **2) for j in xRi]\n",
    "    \n",
    "    yploti = interpolation(xplot, xRi, yRi)\n",
    "    string = \"P_orden \" + str(i-1)\n",
    "    plt.plot(xplot, yploti, '-', label = string)\n",
    "\n",
    "#    idx = np.argwhere(np.diff(np.sign(f - g))).flatten()\n",
    "#    plt.plot(x[idx], f[idx], 'ro')\n",
    "\n",
    "    plt.legend()\n",
    "\n",
    "plt.plot(x, y, '-', label = 'Exacta')\n",
    "plt.title('Fenómeno de Runge')\n",
    "plt.grid('True')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se observa que a medida que se aumenta el orden del polinomio, intentando obtener un mejor ajuste, se presentan oscilaciones en los puntos extremos.\n",
    "\n",
    "***Como actividad complementaria, se invita al estudiante reproducir el fenómeno de Runge que se representa en la gráfica de la Figura 18.14 (p. 526) que se encuentra en el libro de Chapra y Canale, 5a Ed.***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolación mediante trazadores"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introducción"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En las secciones anteriores, se usaron polinomios de $n$-ésimo grado para interpolar entre $n+1$ puntos que se tenían como datos, por ejemplo, para $10$ puntos se puede obtener un polinomio exacto de noveno grado. Esta curva podría agrupar todas las curvas al menos hasta, e incluso, la novena derivada. No obstante, hay casos, como el que se acaba de observar, donde estas funciones llevan a resultados erróneos a causa de los errores de redondeo y los puntos lejanos (*fenómeno de Runge*). \n",
    "\n",
    "Como alternativa para intentar mitigar esta situación se pueden implementar polinomios de menor grado en subconjuntos de los datos. Tales polinomios conectores se denominan trazadores o [splines](https://en.wikipedia.org/wiki/Spline_(mathematics)).\n",
    "\n",
    "Supongamos que empleamos curvas de tercer grado para unir dos conjuntos de datos, cada una de esas funciones se pueden construir de tal forma que las conexiones entre ecuaciones cúbicas adyacentes resulten visualmente suaves. Podría parecer que la aproximación de tercer grado de los trazadores sería inferior a la expresión de noveno grado, entonces, por qué un trazador resulta preferible?\n",
    "\n",
    "El concepto de trazador se originó en la técnica de dibujo que usa una cinta delgada y flexible (llamada spline, en inglés), para dibujar curvas suaves a través de un conjunto de puntos. Se coloca un papel sobre una mesa y alfileres en el papel en la ubicación de los datos. Una curva cúbica suave resulta al entrelazar la cinta entre los alfileres. De aquí que se haya adoptado el nombre de Trazador Cúbico para los polinomios de este tipo.\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C04_Img04_Spline.PNG?raw=true\" width=\"350\" />\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trazadores lineales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La unión más simple entre dos puntos es una línea recta. Los trazadores de primer grado para un grupo de datos ordenados pueden definirse como un conjunto de funciones lineales:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "f(x) & = f(x_0) + m_0(x - x_0), \\quad x_0 \\le x \\le x_1 \\\\\n",
    "f(x) & = f(x_1) + m_1(x - x_1), \\quad x_1 \\le x \\le x_2 \\\\\n",
    "f(x) & = f(x_2) + m_2(x - x_2), \\quad x_2 \\le x \\le x_3 \\\\\n",
    "&\\vdots \\\\\n",
    "f(x) & = f(x_{n-1}) + m_{n-1}(x - x_{n-1}), \\quad x_{n-1} \\le x \\le x_n\n",
    "\\end{split}\n",
    "\\label{eq:Ec4_21} \\tag{4.21}\n",
    "\\end{equation*}\n",
    "\n",
    "donde $m_i=\\frac{f(x_{i+1})-f(x_i)}{(x_{i+1}-x_i)}$ es la pendiente de la línea recta que une los puntos. La principal desventaja de los trazadores de primer grado es que no son suaves. En los puntos donde se encuentran dos trazadores, la pendiente cambia de forma abrupta. La primer derivada de la función es discontinua en esos puntos."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trazadores cuadráticos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para asegurar que las derivadas $m$-ésimas sean continuas en los nodos, se debe emplear un trazador de un grado de, al menos, $m+1$. El objetivo de los trazadores cuadráticos es obtener un polinomio de segundo grado para cada intervalo entre los datos. De manera general, el polinomio en cada intervalo se representa como:\n",
    "\n",
    "<a id='Ec4_22'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "f(x_i) = a_ix^2+b_ix+c_i\n",
    "\\end{split}\n",
    "\\label{eq:Ec4_22} \\tag{4.22}\n",
    "\\end{equation*}\n",
    "\n",
    "Para $n+1$ datos ($i=0, 1, 2,\\ldots, n$) existen $n$ intervalos y, en consecuencia, $3n$ constantes desconocidas ($a$, $b$ y $c$) por evaluar. Por lo tanto, se requieren $3n$ ecuaciones o condiciones para evaluar las incógnitas. Éstas son:\n",
    "\n",
    "1. Los valores de la función de polinomios adyacentes deben ser iguales en los nodos interiores. Esta condición se representa como: \n",
    "\n",
    "<a id='Ec4_23'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "a_{i−1}x_{i−1}^2+b_{i−1}x_{i−1}+c_{i−1}&=f(x_{i−1}) \\\\\n",
    "a_i x_{i−1}^2+b_i x_{i−1}+c_i&=f(x_{i−1})\n",
    "\\end{split}\n",
    "\\label{eq:Ec4_23} \\tag{4.23}\n",
    "\\end{equation*}\n",
    "\n",
    "&emsp; &emsp;para $i=2$ a $n$. Como sólo se emplean nodos interiores, las ecuaciones anteriores proporcionan, cada una, $n–1$ condiciones; en total, $2n–2$ condiciones. \n",
    "\n",
    "2. La primera y la última función deben pasar a través de los puntos extremos. Esto agrega dos ecuaciones más:\n",
    "\n",
    "<a id='Ec4_24'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "a_{1}x_{0}^2+b_{1}x_{0}+c_{1}&=f(x_{0}) \\\\\n",
    "a_n x_{n}^2+b_n x_{n}+c_n&=f(x_{n})\n",
    "\\end{split}\n",
    "\\label{eq:Ec4_24} \\tag{4.24}\n",
    "\\end{equation*}\n",
    "\n",
    "&emsp; &emsp;En total se tienen $2n−2+2=2n$ condiciones\n",
    "\n",
    "3. La primera derivada de la ecuación: $f_i(x)=a_ix^2+b_i x+c$ es: $f'(x_𝑖)=2a_ix+b$. De manera general, esta condición se representa como:\n",
    "\n",
    "<a id='Ec4_25'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "2a_{i−1} x_{i−1}+b_{i−1}=2a_i x_{i−1}+b_i\n",
    "\\end{split}\n",
    "\\label{eq:Ec4_25} \\tag{4.25}\n",
    "\\end{equation*}\n",
    "\n",
    "&emsp; &emsp;para $i = 2an$. Esto proporciona otras $n–1$ condiciones, llegando a un total de $2n+n−1=3n−1$. Como se tienen $3n$ incógnitas, falta una condición más.\n",
    "\n",
    "4. Suponga que en el primer punto la segunda derivada es cero. Esto es: $f''_𝑖(x)=2a_i$, que se puede expresar matemáticamente como: \n",
    "\n",
    "<a id='Ec4_26'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "a_1=0\n",
    "\\end{split}\n",
    "\\label{eq:Ec4_26} \\tag{4.26}\n",
    "\\end{equation*}\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C04_Img05_Spline2.PNG?raw=true\" width=\"350\" />\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>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplo Trazadores Cuadráticos "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dado el siguiente conjunto de datos:\n",
    "\n",
    "|x|f(x)|\n",
    "|:--:|:--:|\n",
    "|3.0|2.5|\n",
    "|4.5|1.0|\n",
    "|7.0|2.5|\n",
    "|9.0|0.5|\n",
    "\n",
    "Evalúe el valor en $x=5.0$ empleando trazadores cuadráticos y cúbicos\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "De la tabla se tienen cuatro puntos y $n=3$ intervalos, por lo que se deben determinar $3n=3\\times3=9$ incógnitas.\n",
    "\n",
    "1. De las ecuaciones [(4.23)](#Ec4_23) se determinan $2\\times3-2=4$ condiciones así,\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "4.5^{2}a_{1}+4.5b_{1}+c_{1}&=1.0 \\\\\n",
    "4.5^{2}a_{2}+4.5b_{2}+c_{2}&=1.0 \\\\\n",
    "7.0^{2}a_{2}+7.0b_{2}+c_{2}&=2.5 \\\\\n",
    "7.0^{2}a_{3}+7.0b_{3}+c_{3}&=2.5\n",
    "\\end{split}\n",
    "\\end{equation*}\n",
    "\n",
    "2. La primera y última función pasan por los puntos extremos, agregando 2 ecuaciones más, ecuación [(4.24)](#Ec4_24):\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "3.0^2a_{1}+3.0b_{1}+c_{1}&=2.5 \\\\\n",
    "9.0^2a_{3}+9.0b_{3}+c_{3}&=0.5 \n",
    "\\end{split}\n",
    "\\end{equation*}\n",
    "\n",
    "3. La continuidad de las derivadas crean adicionalmente $3-1=2$ condiciones, ecuación [(4.25)](#Ec4_25):\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "9.0a_{1}+b_{1}=9.0a_{2}+b_{2} \\\\\n",
    "14.0a_{2}+b_{2}=14.0a_{3}+b_{3} \n",
    "\\end{split}\n",
    "\\end{equation*}\n",
    "\n",
    "4. De la ecuación [(4.26)](#Ec4_26) se obtiene la ecuación faltante, es decir,\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "a_{1}=0\n",
    "\\end{split}\n",
    "\\end{equation*}\n",
    "\n",
    "Esta última ecuación especifica de forma exacta el valor de una de las $9$ condiciones requridas, por lo que el problema se reduce a encontrar las restantes $8$ ecuaciones. Estas condiciones se expresan de forma matricial de la siguiente manera:\n",
    "\n",
    "\\begin{align*}\n",
    "\\left[\\begin{array}{cccc}\n",
    "  4.5 & 1 & 0     & 0   & 0 & 0 & 0 & 0 \\\\\n",
    "  0   & 0 & 20.25 & 4.5 & 1 & 0 & 0 & 0 \\\\\n",
    "  0   & 0 & 49    & 7   & 1 & 0 & 0 & 0 \\\\\n",
    "  0   & 0 & 0     & 0   & 0 & 49 & 7 & 1 \\\\\n",
    "  3   & 1 & 0     & 0   & 0 & 0  & 0 & 0 \\\\\n",
    "  0   & 0 & 0     & 0   & 0 & 81 & 9 & 1 \\\\\n",
    "  1   & 0 & -9    & -1  & 0 & 0 & 0 & 0 \\\\\n",
    "  0   & 0 & 14    & 1   & 0 & -14 & -1 & 0 \\\\\n",
    "\\end{array}\\right]\n",
    "\\begin{Bmatrix}\n",
    "  b_{1}  \\\\\n",
    "  c_{1}  \\\\\n",
    "  a_{2}  \\\\\n",
    "  b_{2}  \\\\\n",
    "  c_{2}  \\\\\n",
    "  a_{3}  \\\\\n",
    "  b_{3}  \\\\\n",
    "  c_{3}  \\\\\n",
    "\\end{Bmatrix}\n",
    "= \\begin{Bmatrix}\n",
    "  1.0  \\\\\n",
    "  1.0  \\\\\n",
    "  2.5  \\\\\n",
    "  2.5  \\\\\n",
    "  2.5  \\\\\n",
    "  0.5  \\\\\n",
    "  0.0  \\\\\n",
    "  0.0\n",
    "\\end{Bmatrix}\n",
    "\\end{align*}\n",
    "\n",
    "Empleando una de las técnicas de resolución de sistemas de ecuaciones lineales vistas en el capítulo anterior, se llega a las respuetas:\n",
    "\n",
    "\\begin{array}{crl}\n",
    "a_1&=0.0, &b_1&=-1, &c_1&=5.5 \\\\\n",
    "a_2&=0.64, &b_2&=-6.76, &c_2&=18.46 \\\\\n",
    "a_3&=-1.6, &b_3&=24.6, &c_3&=-91.3\n",
    "\\end{array}\n",
    "\n",
    "Sustituyendo estos valores en las ecuaciones cuadráticas originales,\n",
    "\n",
    "\\begin{array}{crl}\n",
    "f_1(x)=-x+5.5, &3.0\\le x \\le 4.5 \\\\\n",
    "f_2(x)=0.64x^2-6.76x+18.46, &4.5 \\le x \\le 7.0 \\\\\n",
    "f_3(x)=-1.6x^2+24.6x-91.3, &7.0\\le x \\le 9.0 \\\\\n",
    "\\end{array}\n",
    "\n",
    "Por último, se emplea la ecuación $f_2$, para predecir el valor en $x=5$, que es $f_2(5)=0.66$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trazadores cúbicos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El objetivo es obtener polinomios de tercer grado para cada intervalo entre los nodos:\n",
    "\n",
    "<a id='Ec4_27'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "f_i(x)=a_i x^3+b_i x^2+c_i x+d_i\n",
    "\\end{split}\n",
    "\\label{eq:Ec4_27} \\tag{4.27}\n",
    "\\end{equation*}\n",
    "\n",
    "válida para $x_i \\leq x \\leq x_{i+1}$ para $i = 1, 2, \\ldots, n$. Para encontrar la función de interpolación, primero debemos determinar los coeficientes $a_i, b_i, c_i, d_i$ para cada función cúbica. Así para los $n$ puntos ($i=1,2,\\ldots,n$) existen $n-1$ intervalos y, por lo tanto, se requerirán de $4(n-1)$ condiciones para evaluar las incógnitas. Estas son:\n",
    "\n",
    "1. Los valores de la función deben ser iguales en los nodos interiores: $2n−2$ condiciones.\n",
    "\n",
    "2. La primera y la última función deben pasar a través de los puntos extremos: $2$ condiciones.\n",
    "\n",
    "3. Las primeras derivadas en los nodos interiores deben ser iguales: $n−1$ condiciones.\n",
    "\n",
    "4. Las segundas derivadas en los nodos interiores deben ser iguales: $n−1$.\n",
    "\n",
    "5. Las segundas derivadas en los nodos extremos son cero: $2$ condiciones.\n",
    "\n",
    "Las ecuaciones son:\n",
    "\n",
    "- Primero sabemos que las funciones cúbicas deben cortar los datos de los puntos de la izquierda y de la derecha:\n",
    "\n",
    "<a id='Ec4_28'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "f_i(x_i)&=y_i, i = 1, 2, \\ldots, n-1 \\\\\n",
    "f_i(x_{i+1})&=y_{i+1}, i = 1, 2, \\ldots, n-1\n",
    "\\end{split}\n",
    "\\label{eq:Ec4_28} \\tag{4.28}\n",
    "\\end{equation*}\n",
    "\n",
    "Lo que nos da $2(n−1)$ ecuaciones. A continuación, queremos que cada función cúbica se una lo más suavemente posible con sus vecinas, por lo que restringimos los splines para que tengan derivadas primera y segunda continuas en los puntos de datos $i=2,…,n−1$.\n",
    "\n",
    "<a id='Ec4_29'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "f'_i(x_{i+1})&=f'_{i+1}(x_{i+1}), i = 1, 2, \\ldots, n-2 \\\\\n",
    "f''_i(x_{i+1})&=f''_{i+1}(x_{i+1}), i = 1, 2, \\ldots, n-2\n",
    "\\end{split}\n",
    "\\label{eq:Ec4_29} \\tag{4.29}\n",
    "\\end{equation*}\n",
    "\n",
    "que nos da otras $2(n−2)$ ecuaciones.\n",
    "\n",
    "Se requieren dos ecuaciones más para calcular los coeficientes de $S_i(x)$. Estas dos últimas restricciones son arbitrarias y pueden elegirse para que se ajusten a las circunstancias de la interpolación que se realiza. Un conjunto común de restricciones finales es suponer que las segundas derivadas son cero en los puntos extremos. Esto significa que la curva es una \"línea recta\" en los puntos finales. Explícitamente,\n",
    "\n",
    "<a id='Ec4_30'></a>\n",
    "\\begin{equation*}\n",
    "\\begin{split}\n",
    "f''_1(x_1)&=0 \\\\\n",
    "f''_{n-1}(x_n)&=0\n",
    "\\end{split}\n",
    "\\label{eq:Ec4_30} \\tag{4.30}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplo Trazadores Cúbicos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Utilice CubicSpline para trazar la interpolación spline cúbica del conjunto de datos $x = [0, 1, 2]$ e $y = [1, 3, 2]$ para $0 \\leq x \\leq 2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.interpolate import CubicSpline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#plt.style.use('seaborn-poster')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [0, 1, 2]\n",
    "y = [1, 3, 2]\n",
    "\n",
    "# use bc_type = 'natural' adds the constraints as we described above\n",
    "f = CubicSpline(x, y, bc_type='natural')\n",
    "x_new = np.linspace(0, 2, 100)\n",
    "y_new = f(x_new)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize = (10,8))\n",
    "plt.plot(x_new, y_new, 'b')\n",
    "plt.plot(x, y, 'ro')\n",
    "plt.title('Cubic Spline Interpolation')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para determinar los coeficientes de cada función cúbica, escribimos las restricciones explícitamente como un sistema de ecuaciones lineales con $4(n−1)$ incógnitas. Para $n$ puntos de datos, las incógnitas son los coeficientes $a_i,b_i,c_i,d_i$ del spline cúbico, $f_i$ que une los puntos $x_i$ y $x_{i+1}$.\n",
    "\n",
    "- Para las restricciones $f_i(x_i)=y_i$ tenemos:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "\\begin{split}\n",
    "    a_1 x_1^3     &+         b_1 x_1^2 &+         c_1 x_1 &+     d_1 &= y_1 \\\\\n",
    "    a_2 x_2^3     &+         b_2 x_2^2 &+         c_2 x_2 &+     d_2 &= y_2 \\\\\n",
    "    \\vdots \\\\\n",
    "a_{n-1} x_{n-1}^3 &+ b_{n-1} x_{n-1}^2 &+ c_{n-1} x_{n-1} &+ d_{n-1} &=y_{n-1}\n",
    "\\end{split}\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "- Para las restricciones $f_i(x_{i+1})=y_{i+1}$ tenemos:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "\\begin{split}\n",
    "    a_1 x_2^3     &+         b_1 x_2^2 &+         c_1 x_2 &+     d_1 &= y_2 \\\\\n",
    "    a_2 x_3^3     &+         b_2 x_3^2 &+         c_2 x_3 &+     d_2 &= y_3 \\\\\n",
    "    \\vdots \\\\\n",
    "a_{n-1} x_{n}^3 &+ b_{n-1} x_{n}^2 &+ c_{n-1} x_{n} &+ d_{n-1} &=y_{n}\n",
    "\\end{split}\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "- Para las restricciones $f'_i(x_{i+1})=f'_{i+1}(x_{i+1})$ tenemos:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "\\begin{split}\n",
    "    3a_1 x_2^2 &+ 2b_1 x_2 &+ c_1 & -3a_2 x_2^2   &- 2b_2 x_2 &- c_2 &= 0 \\\\\n",
    "    3a_2 x_3^2 &+ 2b_2 x_3 &+ c_2 & -3a_3 x_3^2   &- 2b_3 x_3 &- c_3 &= 0 \\\\\n",
    "    \\vdots \\\\\n",
    "3a_{n-2} x_{n-1}^2 &+ 2b_{n-2} x_{n-1} &+ c_{n-2} & -3a_{n-1} x_{n-1}^2 &- 2b_{n-1} x_{n-1} &- c_{n-1} &= 0\n",
    "\\end{split}\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "- Para las restricciones $f''_i(x_{i+1})=f''_{i+1}(x_{i+1})$ tenemos:\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "\\begin{split}\n",
    "        6a_1 x_2 &+     2b_1 &-        6a_2 x_2  &-     2b_2 &= 0 \\\\\n",
    "        6a_2 x_3 &+     2b_2 &-        6a_3 x_3  &-     2b_3 &= 0 \\\\\n",
    "    \\vdots \\\\\n",
    "6a_{n-2} x_{n-1} &+ 2b_{n-2} &-6a_{n-1} x_{n-1}  &- 2b_{n-1} &= 0\n",
    "\\end{split}\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "- finalmente para las restricciones en los puntos finales $f'''_1(x_1)=0$ y $f'''_{n-1}(x_n)=0$, tenemos:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "\\begin{split}\n",
    "    6a_1 x_1 &+ 2b_1 &= 0 \\\\\n",
    "    6a_{n-1} x_{n} &+ 2b_{n-1}&= 0\n",
    "\\end{split}\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "Estas ecuaciones son lineales en los coeficientes desconocidos $a_i,b_i,c_i y d_i$. Podemos ponerlos en forma matricial y resolver los coeficientes de cada spline mediante división por la izquierda. Recuerda que siempre que resolvamos la ecuación matricial $Ax=b$ para $x$, debemos asegurarnos de que $A$ sea cuadrado e invertible. En el caso de encontrar ecuaciones spline cúbicas, la matriz $A$ siempre es cuadrada e invertible siempre que los valores $x_i$ en el conjunto de datos sean únicos.\n",
    "\n",
    "Encuentre la interpolación spline cúbica en $x = 1.5$ según los datos $x = [0, 1, 2]$, $y = [1, 3, 2]$. Primero creamos el sistema de ecuaciones apropiado y encontramos los coeficientes de los splines cúbicos resolviendo el sistema en forma matricial.\n",
    "\n",
    "La forma matricial del sistema de ecuaciones es:\n",
    "\n",
    "\\begin{align*}\n",
    "\\left[\\begin{array}{cccc}\n",
    "  0 & 0 & 0 & 1 &  0 &  0 &  0 & 0 \\\\\n",
    "  0 & 0 & 0 & 0 &  1 &  1 &  1 & 1 \\\\\n",
    "  1 & 1 & 1 & 1 &  0 &  0 &  0 & 0 \\\\\n",
    "  0 & 0 & 0 & 0 &  8 &  4 &  2 & 1 \\\\\n",
    "  3 & 2 & 1 & 0 & -3 & -2 & -1 & 0 \\\\\n",
    "  6 & 2 & 0 & 0 & -6 & -2 &  0 & 0 \\\\\n",
    "  0 & 2 & 0 & 0 &  0 &  0 &  0 & 0 \\\\\n",
    "  0 & 0 & 0 & 0 & 12 &  2 &  0 & 0\n",
    "\\end{array}\\right]\n",
    "\\begin{Bmatrix}\n",
    "  a_{1}  \\\\\n",
    "  b_{1}  \\\\\n",
    "  c_{1}  \\\\\n",
    "  d_{1}  \\\\\n",
    "  a_{2}  \\\\\n",
    "  b_{2}  \\\\\n",
    "  c_{2}  \\\\\n",
    "  d_{2}  \n",
    "\\end{Bmatrix}\n",
    "= \\begin{Bmatrix}\n",
    "  1 \\\\\n",
    "  3 \\\\\n",
    "  3 \\\\\n",
    "  2 \\\\\n",
    "  0 \\\\\n",
    "  0 \\\\\n",
    "  0 \\\\\n",
    "  0 \n",
    "\\end{Bmatrix}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.array([1, 3, 3, 2, 0, 0, 0, 0])\n",
    "b = b[:, np.newaxis]\n",
    "A = np.array([[0, 0, 0, 1, 0, 0, 0, 0], \n",
    "              [0, 0, 0, 0, 1, 1, 1, 1], \n",
    "              [1, 1, 1, 1, 0, 0, 0, 0], \n",
    "              [0, 0, 0, 0, 8, 4, 2, 1], \n",
    "              [3, 2, 1, 0, -3, -2, -1, 0], \n",
    "              [6, 2, 0, 0, -6, -2, 0, 0],\n",
    "              [0, 2, 0, 0, 0, 0, 0, 0], \n",
    "              [0, 0, 0, 0, 12, 2, 0, 0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "np.dot(np.linalg.inv(A), b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "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": 4
}