Cap03_SistemaEcuacionesLineales.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\">Métodos Numéricos</h1>\n",
    "<h1 align=\"center\">Capítulo 3: Sistemas de Ecuaciones Lineales</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/Cap03_SistemaEcuacionesLineales.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": {
    "toc": true
   },
   "source": [
    "<h1>Tabla de Contenidos<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Sistema-de-Ecuaciones-Lineales\" data-toc-modified-id=\"Sistema-de-Ecuaciones-Lineales-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Sistema de Ecuaciones Lineales</a></span><ul class=\"toc-item\"><li><span><a href=\"#Introducción\" data-toc-modified-id=\"Introducción-1.1\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>Introducción</a></span></li><li><span><a href=\"#Singularidad-y-No-Singularidad\" data-toc-modified-id=\"Singularidad-y-No-Singularidad-1.2\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>Singularidad y No Singularidad</a></span></li><li><span><a href=\"#Existencia-y-Unicidad\" data-toc-modified-id=\"Existencia-y-Unicidad-1.3\"><span class=\"toc-item-num\">1.3&nbsp;&nbsp;</span>Existencia y Unicidad</a></span></li><li><span><a href=\"#Interpretación-geométrica\" data-toc-modified-id=\"Interpretación-geométrica-1.4\"><span class=\"toc-item-num\">1.4&nbsp;&nbsp;</span>Interpretación geométrica</a></span></li></ul></li><li><span><a href=\"#Métodos-directos-o-de-eliminación\" data-toc-modified-id=\"Métodos-directos-o-de-eliminación-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Métodos directos o de eliminación</a></span><ul class=\"toc-item\"><li><span><a href=\"#Eliminación-de-Gauss-Simple\" data-toc-modified-id=\"Eliminación-de-Gauss-Simple-2.1\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Eliminación de Gauss Simple</a></span><ul class=\"toc-item\"><li><span><a href=\"#Introducción\" data-toc-modified-id=\"Introducción-2.1.1\"><span class=\"toc-item-num\">2.1.1&nbsp;&nbsp;</span>Introducción</a></span></li><li><span><a href=\"#Algoritmo-general\" data-toc-modified-id=\"Algoritmo-general-2.1.2\"><span class=\"toc-item-num\">2.1.2&nbsp;&nbsp;</span>Algoritmo general</a></span></li><li><span><a href=\"#Etapas-en-el-procesos-de-eliminación\" data-toc-modified-id=\"Etapas-en-el-procesos-de-eliminación-2.1.3\"><span class=\"toc-item-num\">2.1.3&nbsp;&nbsp;</span>Etapas en el procesos de eliminación</a></span></li><li><span><a href=\"#Algoritmo-proceso-de-Eliminación\" data-toc-modified-id=\"Algoritmo-proceso-de-Eliminación-2.1.4\"><span class=\"toc-item-num\">2.1.4&nbsp;&nbsp;</span>Algoritmo proceso de Eliminación</a></span></li><li><span><a href=\"#Resultado-del-proceso-de-Eliminación\" data-toc-modified-id=\"Resultado-del-proceso-de-Eliminación-2.1.5\"><span class=\"toc-item-num\">2.1.5&nbsp;&nbsp;</span>Resultado del proceso de Eliminación</a></span></li><li><span><a href=\"#Implementación-computacional-proceso-de-Eliminación\" data-toc-modified-id=\"Implementación-computacional-proceso-de-Eliminación-2.1.6\"><span class=\"toc-item-num\">2.1.6&nbsp;&nbsp;</span>Implementación computacional proceso de Eliminación</a></span></li><li><span><a href=\"#Sustitución-regresiva-(o-despeje)\" data-toc-modified-id=\"Sustitución-regresiva-(o-despeje)-2.1.7\"><span class=\"toc-item-num\">2.1.7&nbsp;&nbsp;</span>Sustitución regresiva (o despeje)</a></span></li><li><span><a href=\"#Algoritmo-proceso-de-sustitución-regresiva\" data-toc-modified-id=\"Algoritmo-proceso-de-sustitución-regresiva-2.1.8\"><span class=\"toc-item-num\">2.1.8&nbsp;&nbsp;</span>Algoritmo proceso de sustitución regresiva</a></span></li><li><span><a href=\"#Implementación-computacional\" data-toc-modified-id=\"Implementación-computacional-2.1.9\"><span class=\"toc-item-num\">2.1.9&nbsp;&nbsp;</span>Implementación computacional</a></span></li><li><span><a href=\"#Ejemplo-teórico-algoritmo-Eliminación-de-Gauss-simple\" data-toc-modified-id=\"Ejemplo-teórico-algoritmo-Eliminación-de-Gauss-simple-2.1.10\"><span class=\"toc-item-num\">2.1.10&nbsp;&nbsp;</span>Ejemplo teórico algoritmo Eliminación de Gauss simple</a></span></li><li><span><a href=\"#Ejemplo-numérico-algoritmo-Eliminación-de-Gauss-Simple\" data-toc-modified-id=\"Ejemplo-numérico-algoritmo-Eliminación-de-Gauss-Simple-2.1.11\"><span class=\"toc-item-num\">2.1.11&nbsp;&nbsp;</span>Ejemplo numérico algoritmo Eliminación de Gauss Simple</a></span></li><li><span><a href=\"#Ejemplo-computacional-algoritmo-de-Eliminación-de-Gauss-Simple\" data-toc-modified-id=\"Ejemplo-computacional-algoritmo-de-Eliminación-de-Gauss-Simple-2.1.12\"><span class=\"toc-item-num\">2.1.12&nbsp;&nbsp;</span>Ejemplo computacional algoritmo de Eliminación de Gauss Simple</a></span></li></ul></li><li><span><a href=\"#Método-de-Gauss---Jordan\" data-toc-modified-id=\"Método-de-Gauss---Jordan-2.2\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span>Método de Gauss - Jordan</a></span></li><li><span><a href=\"#Aspectos-computacionales,-inconvenitentes-de-los-métodos-de-eliminación\" data-toc-modified-id=\"Aspectos-computacionales,-inconvenitentes-de-los-métodos-de-eliminación-2.3\"><span class=\"toc-item-num\">2.3&nbsp;&nbsp;</span>Aspectos computacionales, inconvenitentes de los métodos de eliminación</a></span><ul class=\"toc-item\"><li><span><a href=\"#Análisis-de-desempeño-algoritmo-de-Eliminación-de-Gauss-Simple\" data-toc-modified-id=\"Análisis-de-desempeño-algoritmo-de-Eliminación-de-Gauss-Simple-2.3.1\"><span class=\"toc-item-num\">2.3.1&nbsp;&nbsp;</span>Análisis de desempeño algoritmo de Eliminación de Gauss Simple</a></span></li><li><span><a href=\"#División-entre-cero\" data-toc-modified-id=\"División-entre-cero-2.3.2\"><span class=\"toc-item-num\">2.3.2&nbsp;&nbsp;</span>División entre cero</a></span></li><li><span><a href=\"#Errores-de-redondeo\" data-toc-modified-id=\"Errores-de-redondeo-2.3.3\"><span class=\"toc-item-num\">2.3.3&nbsp;&nbsp;</span>Errores de redondeo</a></span></li><li><span><a href=\"#Condicionamiento-del-sistema\" data-toc-modified-id=\"Condicionamiento-del-sistema-2.3.4\"><span class=\"toc-item-num\">2.3.4&nbsp;&nbsp;</span>Condicionamiento del sistema</a></span></li><li><span><a href=\"#Escalamiento\" data-toc-modified-id=\"Escalamiento-2.3.5\"><span class=\"toc-item-num\">2.3.5&nbsp;&nbsp;</span>Escalamiento</a></span></li></ul></li><li><span><a href=\"#Técnicas-para-mejorar-la-solución-en-los-métodos-de-Eliminación\" data-toc-modified-id=\"Técnicas-para-mejorar-la-solución-en-los-métodos-de-Eliminación-2.4\"><span class=\"toc-item-num\">2.4&nbsp;&nbsp;</span>Técnicas para mejorar la solución en los métodos de Eliminación</a></span><ul class=\"toc-item\"><li><span><a href=\"#Escalamiento\" data-toc-modified-id=\"Escalamiento-2.4.1\"><span class=\"toc-item-num\">2.4.1&nbsp;&nbsp;</span>Escalamiento</a></span></li><li><span><a href=\"#Pivoteo\" data-toc-modified-id=\"Pivoteo-2.4.2\"><span class=\"toc-item-num\">2.4.2&nbsp;&nbsp;</span>Pivoteo</a></span></li></ul></li><li><span><a href=\"#Factorización-LU\" data-toc-modified-id=\"Factorización-LU-2.5\"><span class=\"toc-item-num\">2.5&nbsp;&nbsp;</span>Factorización LU</a></span><ul class=\"toc-item\"><li><span><a href=\"#Introducción\" data-toc-modified-id=\"Introducción-2.5.1\"><span class=\"toc-item-num\">2.5.1&nbsp;&nbsp;</span>Introducción</a></span></li><li><span><a href=\"#Determinación-de-la-descomposición-LU\" data-toc-modified-id=\"Determinación-de-la-descomposición-LU-2.5.2\"><span class=\"toc-item-num\">2.5.2&nbsp;&nbsp;</span>Determinación de la descomposición LU</a></span></li><li><span><a href=\"#Estrategia-de-solución\" data-toc-modified-id=\"Estrategia-de-solución-2.5.3\"><span class=\"toc-item-num\">2.5.3&nbsp;&nbsp;</span>Estrategia de solución</a></span></li><li><span><a href=\"#Obtención-de-las-matrices-L-y-U\" data-toc-modified-id=\"Obtención-de-las-matrices-L-y-U-2.5.4\"><span class=\"toc-item-num\">2.5.4&nbsp;&nbsp;</span>Obtención de las matrices <em>L</em> y <em>U</em></a></span></li><li><span><a href=\"#Ejemplo-LU\" data-toc-modified-id=\"Ejemplo-LU-2.5.5\"><span class=\"toc-item-num\">2.5.5&nbsp;&nbsp;</span>Ejemplo <em>LU</em></a></span></li><li><span><a href=\"#Implementación-computacional-del-algoritmo-LU\" data-toc-modified-id=\"Implementación-computacional-del-algoritmo-LU-2.5.6\"><span class=\"toc-item-num\">2.5.6&nbsp;&nbsp;</span>Implementación computacional del algoritmo LU</a></span></li><li><span><a href=\"#Ejemplo-computacional-algoritmo-LU:-Cálculo-de-la-Inversa\" data-toc-modified-id=\"Ejemplo-computacional-algoritmo-LU:-Cálculo-de-la-Inversa-2.5.7\"><span class=\"toc-item-num\">2.5.7&nbsp;&nbsp;</span>Ejemplo computacional algoritmo LU: Cálculo de la Inversa</a></span></li></ul></li></ul></li><li><span><a href=\"#Análisis-de-Error-y-condición-del-Sistema\" data-toc-modified-id=\"Análisis-de-Error-y-condición-del-Sistema-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Análisis de Error y condición del Sistema</a></span><ul class=\"toc-item\"><li><span><a href=\"#Condicionamiento-del-sistema\" data-toc-modified-id=\"Condicionamiento-del-sistema-3.1\"><span class=\"toc-item-num\">3.1&nbsp;&nbsp;</span>Condicionamiento del sistema</a></span></li><li><span><a href=\"#Normas-vectoriales\" data-toc-modified-id=\"Normas-vectoriales-3.2\"><span class=\"toc-item-num\">3.2&nbsp;&nbsp;</span>Normas vectoriales</a></span><ul class=\"toc-item\"><li><span><a href=\"#Norma-1\" data-toc-modified-id=\"Norma-1-3.2.1\"><span class=\"toc-item-num\">3.2.1&nbsp;&nbsp;</span>Norma-1</a></span></li><li><span><a href=\"#Norma-2\" data-toc-modified-id=\"Norma-2-3.2.2\"><span class=\"toc-item-num\">3.2.2&nbsp;&nbsp;</span>Norma-2</a></span></li><li><span><a href=\"#Norma-$\\infty$\" data-toc-modified-id=\"Norma-$\\infty$-3.2.3\"><span class=\"toc-item-num\">3.2.3&nbsp;&nbsp;</span>Norma-$\\infty$</a></span></li><li><span><a href=\"#Conclusión\" data-toc-modified-id=\"Conclusión-3.2.4\"><span class=\"toc-item-num\">3.2.4&nbsp;&nbsp;</span>Conclusión</a></span></li></ul></li><li><span><a href=\"#Normas-matriciales\" data-toc-modified-id=\"Normas-matriciales-3.3\"><span class=\"toc-item-num\">3.3&nbsp;&nbsp;</span>Normas matriciales</a></span><ul class=\"toc-item\"><li><span><a href=\"#Norma-1\" data-toc-modified-id=\"Norma-1-3.3.1\"><span class=\"toc-item-num\">3.3.1&nbsp;&nbsp;</span>Norma-1</a></span></li><li><span><a href=\"#Norma-2\" data-toc-modified-id=\"Norma-2-3.3.2\"><span class=\"toc-item-num\">3.3.2&nbsp;&nbsp;</span>Norma-2</a></span></li><li><span><a href=\"#Norma-$\\infty$\" data-toc-modified-id=\"Norma-$\\infty$-3.3.3\"><span class=\"toc-item-num\">3.3.3&nbsp;&nbsp;</span>Norma-$\\infty$</a></span></li><li><span><a href=\"#Conclusión\" data-toc-modified-id=\"Conclusión-3.3.4\"><span class=\"toc-item-num\">3.3.4&nbsp;&nbsp;</span>Conclusión</a></span></li></ul></li><li><span><a href=\"#Número-de-condición-de-una-matriz\" data-toc-modified-id=\"Número-de-condición-de-una-matriz-3.4\"><span class=\"toc-item-num\">3.4&nbsp;&nbsp;</span>Número de condición de una matriz</a></span><ul class=\"toc-item\"><li><span><a href=\"#Definición\" data-toc-modified-id=\"Definición-3.4.1\"><span class=\"toc-item-num\">3.4.1&nbsp;&nbsp;</span>Definición</a></span></li><li><span><a href=\"#Propiedades\" data-toc-modified-id=\"Propiedades-3.4.2\"><span class=\"toc-item-num\">3.4.2&nbsp;&nbsp;</span>Propiedades</a></span></li><li><span><a href=\"#Cálculo-del-número-de-condición\" data-toc-modified-id=\"Cálculo-del-número-de-condición-3.4.3\"><span class=\"toc-item-num\">3.4.3&nbsp;&nbsp;</span>Cálculo del número de condición</a></span></li></ul></li><li><span><a href=\"#Precisión-de-las-soluciones\" data-toc-modified-id=\"Precisión-de-las-soluciones-3.5\"><span class=\"toc-item-num\">3.5&nbsp;&nbsp;</span>Precisión de las soluciones</a></span><ul class=\"toc-item\"><li><span><a href=\"#Residual-de-una-solución\" data-toc-modified-id=\"Residual-de-una-solución-3.5.1\"><span class=\"toc-item-num\">3.5.1&nbsp;&nbsp;</span>Residual de una solución</a></span></li><li><span><a href=\"#Ejemplo-numérico\" data-toc-modified-id=\"Ejemplo-numérico-3.5.2\"><span class=\"toc-item-num\">3.5.2&nbsp;&nbsp;</span>Ejemplo numérico</a></span></li><li><span><a href=\"#Estimación-de-la-precisión\" data-toc-modified-id=\"Estimación-de-la-precisión-3.5.3\"><span class=\"toc-item-num\">3.5.3&nbsp;&nbsp;</span>Estimación de la precisión</a></span></li><li><span><a href=\"#Interpretación-geométrica\" data-toc-modified-id=\"Interpretación-geométrica-3.5.4\"><span class=\"toc-item-num\">3.5.4&nbsp;&nbsp;</span>Interpretación geométrica</a></span></li><li><span><a href=\"#Conclusión\" data-toc-modified-id=\"Conclusión-3.5.5\"><span class=\"toc-item-num\">3.5.5&nbsp;&nbsp;</span>Conclusión</a></span></li><li><span><a href=\"#Ejemplo-numérico---Matriz-de-Hilbert\" data-toc-modified-id=\"Ejemplo-numérico---Matriz-de-Hilbert-3.5.6\"><span class=\"toc-item-num\">3.5.6&nbsp;&nbsp;</span>Ejemplo numérico - Matriz de Hilbert</a></span></li></ul></li><li><span><a href=\"#Tipos-especiales-de-sistemas-lineales\" data-toc-modified-id=\"Tipos-especiales-de-sistemas-lineales-3.6\"><span class=\"toc-item-num\">3.6&nbsp;&nbsp;</span>Tipos especiales de sistemas lineales</a></span><ul class=\"toc-item\"><li><span><a href=\"#Simétrica\" data-toc-modified-id=\"Simétrica-3.6.1\"><span class=\"toc-item-num\">3.6.1&nbsp;&nbsp;</span>Simétrica</a></span></li><li><span><a href=\"#Dispersa\" data-toc-modified-id=\"Dispersa-3.6.2\"><span class=\"toc-item-num\">3.6.2&nbsp;&nbsp;</span>Dispersa</a></span></li><li><span><a href=\"#Definida-positiva\" data-toc-modified-id=\"Definida-positiva-3.6.3\"><span class=\"toc-item-num\">3.6.3&nbsp;&nbsp;</span>Definida positiva</a></span><ul class=\"toc-item\"><li><span><a href=\"#Ejemplo:-Factorización-de-Cholesky\" data-toc-modified-id=\"Ejemplo:-Factorización-de-Cholesky-3.6.3.1\"><span class=\"toc-item-num\">3.6.3.1&nbsp;&nbsp;</span>Ejemplo: Factorización de Cholesky</a></span></li><li><span><a href=\"#Implementación-computacional-algoritmo-de-cholesky\" data-toc-modified-id=\"Implementación-computacional-algoritmo-de-cholesky-3.6.3.2\"><span class=\"toc-item-num\">3.6.3.2&nbsp;&nbsp;</span>Implementación computacional algoritmo de cholesky</a></span></li></ul></li><li><span><a href=\"#Matriz-de-banda\" data-toc-modified-id=\"Matriz-de-banda-3.6.4\"><span class=\"toc-item-num\">3.6.4&nbsp;&nbsp;</span>Matriz de banda</a></span></li></ul></li></ul></li><li><span><a href=\"#Métodos-Iterativos\" data-toc-modified-id=\"Métodos-Iterativos-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Métodos Iterativos</a></span><ul class=\"toc-item\"><li><span><a href=\"#Introducción\" data-toc-modified-id=\"Introducción-4.1\"><span class=\"toc-item-num\">4.1&nbsp;&nbsp;</span>Introducción</a></span></li><li><span><a href=\"#Método-de-Jacobi\" data-toc-modified-id=\"Método-de-Jacobi-4.2\"><span class=\"toc-item-num\">4.2&nbsp;&nbsp;</span>Método de <em>Jacobi</em></a></span><ul class=\"toc-item\"><li><span><a href=\"#Algoritmo-del-método-de-Jacobi\" data-toc-modified-id=\"Algoritmo-del-método-de-Jacobi-4.2.1\"><span class=\"toc-item-num\">4.2.1&nbsp;&nbsp;</span>Algoritmo del método de Jacobi</a></span></li></ul></li><li><span><a href=\"#Método-de-Gauss---Seidel\" data-toc-modified-id=\"Método-de-Gauss---Seidel-4.3\"><span class=\"toc-item-num\">4.3&nbsp;&nbsp;</span>Método de Gauss - Seidel</a></span></li><li><span><a href=\"#Forma-matricial-de-los-métodos-iterativos\" data-toc-modified-id=\"Forma-matricial-de-los-métodos-iterativos-4.4\"><span class=\"toc-item-num\">4.4&nbsp;&nbsp;</span>Forma matricial de los métodos iterativos</a></span></li><li><span><a href=\"#Convergencia-de-los-métodos-iterativos\" data-toc-modified-id=\"Convergencia-de-los-métodos-iterativos-4.5\"><span class=\"toc-item-num\">4.5&nbsp;&nbsp;</span>Convergencia de los métodos iterativos</a></span></li><li><span><a href=\"#Método-SOR\" data-toc-modified-id=\"Método-SOR-4.6\"><span class=\"toc-item-num\">4.6&nbsp;&nbsp;</span>Método SOR</a></span><ul class=\"toc-item\"><li><span><a href=\"#Algoritmo-método-SOR\" data-toc-modified-id=\"Algoritmo-método-SOR-4.6.1\"><span class=\"toc-item-num\">4.6.1&nbsp;&nbsp;</span>Algoritmo método SOR</a></span></li></ul></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sistema de Ecuaciones Lineales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introducción"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un [Sistema de Ecuaciones Lineales](https://en.wikipedia.org/wiki/System_of_linear_equations), *SEL*, es un conjunto de $m$ ecuaciones lineales con $n$ incógnitas, cuya solución es un conjunto de valores para las incógnitas con el que se satisfacen todas las ecuaciones.\n",
    "\n",
    "Asumamos que los elementos que componen las ecuaciones están definidos en los $\\mathbb{R}$ y que siempre tendremos la misma cantidad de ecuaciones que de incógnitas, es decir $m=n$.\n",
    "\n",
    "$$a_{11}x_1 + a_{12}x_2 + a_{13}x_3 + \\ldots + a_{1n}x_n = b_1 \\hspace{1 cm} (ec. 1)$$\n",
    "\n",
    "$$a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + \\ldots + a_{2n}x_n = b_2 \\hspace{1 cm} (ec. 2)$$\n",
    "\n",
    "$$a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + \\ldots + a_{3n}x_n = b_3 \\hspace{1 cm} (ec. 3)$$\n",
    "\n",
    "$$\\vdots$$\n",
    "\n",
    "$$a_{n1}x_1 + a_{n2}x_2 + a_{n3}x_3 + \\ldots + a_{nn}x_n = b_n \\hspace{1 cm} (ec. n)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En notación matricial, un sistema de ecuaciones algebraicas lineal tiene la forma:\n",
    "\n",
    "$$\\mathbf{Ax=b}$$\n",
    "\n",
    "donde\n",
    "\n",
    "- $\\mathbf{A}$: Matriz de coeficientes $[n \\times n]$\n",
    "\n",
    "\n",
    "- $\\mathbf{b}$: vector de términos independientes $[n \\times 1]$\n",
    "\n",
    "\n",
    "- $\\mathbf{x}$: Vector de incógnitas $[n \\times 1]$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "\\left[\\begin{array}{cccc}\n",
    "  a_{11} & a_{12} & \\ldots &a_{1n} \\\\\n",
    "  a_{21} & a_{22} & \\ldots &a_{2n} \\\\\n",
    "  a_{31} & a_{32} & \\ldots &a_{3n} \\\\\n",
    "  \\vdots & \\vdots & \\ddots &\\vdots \\\\\n",
    "  a_{n1} & a_{n2} & \\ldots &a_{nn} \\\\\n",
    "\\end{array}\\right]\n",
    "\\begin{Bmatrix}\n",
    "  x_{1}  \\\\\n",
    "  x_{2}  \\\\\n",
    "  x_{3}  \\\\\n",
    "  \\vdots \\\\\n",
    "  x_{n}\n",
    "\\end{Bmatrix}\n",
    "= \\begin{Bmatrix}\n",
    "  b_{1}  \\\\\n",
    "  b_{2}  \\\\\n",
    "  b_{3}  \\\\\n",
    "  \\vdots \\\\\n",
    "  b_{n}\n",
    "\\end{Bmatrix}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lo siguiente conduce a la pregunta *¿Puede el vector $\\mathbf{b}$ ser expresado como una combinación lineal de las columnas de la matriz $\\mathbf{A}$?*\n",
    "\n",
    "- Los coeficientes de esta combinación lineal están dados por las componentes del vector solución $\\mathbf{x}$.\n",
    "\n",
    "\n",
    "- Puede o no tener solución.\n",
    "\n",
    "\n",
    "- Puede no ser única.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Singularidad y No Singularidad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una matriz $\\mathbf{A}_{[n,n]}$ se dice que es singular si presenta una de las siguientes propiedades:\n",
    "\n",
    "\n",
    "- $\\mathbf{A}$ No tiene inversa: No existe una matriz $\\mathbf{M}$ tal que $\\mathbf{AM = MA = I}$.\n",
    "\n",
    "\n",
    "- $\\text{det}(\\mathbf{A}) = 0$\n",
    "\n",
    "\n",
    "- $\\text{Rango}(\\mathbf{A}) < n$: El rango de una matriz es el máximo número de filas o columnas linealmente independientes.\n",
    "\n",
    "\n",
    "- $\\mathbf{Az = 0}$, para cualquier vector $\\mathbf{z} \\neq \\mathbf{0}$.\n",
    "\n",
    "\n",
    "La solubilidad de un sistema de ecuaciones lineales está determinado si la matriz $\\mathbf{A}$ es, o no, singular.\n",
    "\n",
    "- ***No Singular:*** Tiene solución y es única\n",
    "\n",
    "\n",
    "- ***Singular:*** No tiene solución o tiene infinitas soluciones"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Existencia y Unicidad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- La Existencia y la Unicidad de la solución de $\\mathbf{Ax=b}$ depende de si $\\mathbf{A}$ es singular o no singular.\n",
    "\n",
    "\n",
    "- Puede también depender de $\\mathbf{b}$, pero únicamente en el caso singular\n",
    "\n",
    "\n",
    "- Si $\\mathbf{b} \\in \\text{span}(\\mathbf{A})$, el sistema se dice que es consistente\n",
    "\n",
    "\n",
    "|$\\mathbf{A}$ | $\\mathbf{b}$           | #soluciones|\n",
    "|-----------|:------------------------:|-----------|\n",
    "|No singular| puede ser arbitrario     |única      |\n",
    "|singular   | $\\mathbf{b} \\in \\text{span}(\\mathbf{A})$   |infinitas  |\n",
    "|singular   | $\\mathbf{b} \\notin \\text{span}(\\mathbf{A})$|ninguna    |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpretación geométrica"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En $2D$, cada ecuación determina una línea recta en el plano. Si dos líneas rectas son no paralelas (no singular), entonces el punto de intersección es único. La solución es el punto de intersección de las dos líneas\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img01_Sln2D_01.PNG?raw=true\" width=\"350\" />\n",
    "</p>\n",
    "\n",
    "\n",
    "Si las dos líneas son paralelas (singular), entonces o las líneas no se intersectan (no hay solución) \n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img02_Sln2D_02.PNG?raw=true\" width=\"350\" />\n",
    "</p>\n",
    "\n",
    "o coinciden en todos los puntos a lo largo de las líneas (infinitas soluciones)\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img03_Sln2D_03.PNG?raw=true\" width=\"350\" />\n",
    "</p>\n",
    "\n",
    "En mayores dimensiones, cada ecuación determina un hiperplano; si la matriz es no singular, la intersección de los hiperplanos es la única solución.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ***Ejemplo: No singularidad***\n",
    "\n",
    "Dado el Sistema $2 \\times 2$\n",
    "\n",
    "$$2x_1 + 3x_2 = b_1 $$\n",
    "$$5x_1 + 4x_2 = b_2 $$\n",
    "\n",
    "en notación matricial se tiene:\n",
    "\n",
    "\\begin{align*}\n",
    "\\left[\\begin{array}{cc}\n",
    "  2 & 3 & \\\\\n",
    "  5 & 4 & \\\\\n",
    " \\end{array}\\right]\n",
    "\\begin{Bmatrix}\n",
    "  x_{1}  \\\\\n",
    "  x_{2}  \\\\\n",
    "\\end{Bmatrix}\n",
    "= \\begin{Bmatrix}\n",
    "  b_{1}  \\\\\n",
    "  b_{2}  \\\\\n",
    "\\end{Bmatrix}\n",
    "\\end{align*}\n",
    "\n",
    "- Es no singular, independientemente del valor de $\\mathbf{b}$\n",
    "\n",
    "\n",
    "- Por ejemplo, si $\\mathbf{b}=\\{ 8 \\hspace{.25cm} 13\\}^T$, entonces $\\mathbf{x}=\\{ 1 \\hspace{0.25cm} 2 \\}^T$ es la solución única del sistema.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ***Ejemplo: Singularidad***\n",
    "\n",
    "Dado el Sistema $2 \\times 2$\n",
    "\n",
    "$$2x_1 + 3x_2 = b_1 $$\n",
    "$$4x_1 + 6x_2 = b_2 $$\n",
    "\n",
    "en notación matricial se tiene:\n",
    "\n",
    "\\begin{align*}\n",
    "\\left[\\begin{array}{cc}\n",
    "  2 & 3 & \\\\\n",
    "  4 & 6 & \\\\\n",
    " \\end{array}\\right]\n",
    "\\begin{Bmatrix}\n",
    "  x_{1}  \\\\\n",
    "  x_{2}  \\\\\n",
    "\\end{Bmatrix}\n",
    "= \\begin{Bmatrix}\n",
    "  b_{1}  \\\\\n",
    "  b_{2}  \\\\\n",
    "\\end{Bmatrix}\n",
    "\\end{align*}\n",
    "\n",
    "- Es singular, independientemente del valor de $\\mathbf{b}$\n",
    "\n",
    "\n",
    "- Por ejemplo, si $\\mathbf{b}=\\{ 4 \\hspace{.25cm} 7\\}^T$, entonces no existe solución.\n",
    "\n",
    "\n",
    "- Si $\\mathbf{b}=\\{ 4 \\hspace{0.25cm} 8 \\}^T$ , $\\mathbf{x}=\\{ \\gamma \\hspace{0.25cm} (4-2\\gamma)/3 \\}^T$ es solución para cualquier número real $\\gamma$, entonces tiene infinitas soluciones"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Ejemplo: mal condicionamiento***\n",
    "\n",
    "Adicionalmente a esto, los sistemas muy próximos a ser singulares, pueden ocasionar muchos problemas para su solución.\n",
    "\n",
    "Dado el Sistema $2 \\times 2$\n",
    "\n",
    "$$-\\frac{2.3}{5}x_1 + x_2 = 1.1 $$\n",
    "$$-\\frac{1}{2}x_1 + x_2 = 1 $$\n",
    "\n",
    "en notación matricial se tiene:\n",
    "\n",
    "\\begin{align*}\n",
    "\\left[\\begin{array}{cc}\n",
    "  -\\frac{2.3}{5} & 1 & \\\\\n",
    "  -\\frac{1}{2} & 1 & \\\\\n",
    " \\end{array}\\right]\n",
    "\\begin{Bmatrix}\n",
    "  x_{1}  \\\\\n",
    "  x_{2}  \\\\\n",
    "\\end{Bmatrix}\n",
    "= \\begin{Bmatrix}\n",
    "  1.1  \\\\\n",
    "  1  \\\\\n",
    "\\end{Bmatrix}\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img03a_Sln2D_04.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",
    "Se observa que, visualmente, no es posible determinar cuál es la solución del sistema. Igual inconveniente podría tener el computador, al no poder determinar (aproximar) cuál es el valor de la solución."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Métodos directos o de eliminación"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Un método para la solución de *SEL* es directo si se define como un conjunto finito de operaciones y de pasos para hallar la solución. El método genera una solución que sería exacta si no fuese por los errores de redondeo al efectuar las operaciones aritméticas (sumas, restas, multiplicaciones y divisiones).\n",
    "\n",
    "Dado un SEL de la forma $\\mathbf{Ax=b}$, se aplican *“operaciones elementales de fila”* para transformar el sistema de ecuaciones en uno equivalente y más simple de resolver.\n",
    "\n",
    "***Operaciones Elementales de Fila:*** Dada una matriz $\\mathbf{A}$, se dice que se realiza una operación elemental de fila si se hace una de las siguientes operaciones:\n",
    "\n",
    "- Intercambio de dos filas;\n",
    "\n",
    "\n",
    "- Sumar a una fila un múltiplo de otra fila; \n",
    "\n",
    "\n",
    "- Multiplicar una fila por un escalar no nulo."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='EGS'></a>\n",
    "### Eliminación de Gauss Simple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Introducción"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El método de [Eliminación de Gauss Simple](https://en.wikipedia.org/wiki/Gaussian_elimination#:~:text=In%20mathematics%2C%20Gaussian%20elimination%2C%20also,the%20corresponding%20matrix%20of%20coefficients.) se descompone en:\n",
    "\n",
    "- ***Eliminación:*** Dado el sistema $\\mathbf{Ax=b}$, se realizan operaciones fundamentales de fila hasta obtener un sistema triangular superior equivalente $\\mathbf{Ux=B}$\n",
    "\n",
    "\n",
    "- ***Despeje:*** Del sistema triangular superior, se halla primero el valor de $x_n$, con este valor se halla $x_{n-1}$, con estos dos valores se halla $x_{n-2}$ y así sucesivamente hasta obtener el valor de $x_{1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Algoritmo general"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color:rgba(0, 0, 0, 0.0470588); padding:10px 0;font-family:monospace;\">\n",
    "Lea <font color = \"red\">A, b</font><br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp; <font color = \"orange\">U, B $\\leftarrow$ Eliminación(A,b) </font><br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp; <font color = \"orange\">x $\\leftarrow$ Despeje(U,B) </font><br>\n",
    "Imprima <font color = \"red\">x</font><br>\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='Etapaseliminacion'></a>\n",
    "#### Etapas en el procesos de eliminación"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ***Eliminación Etapa 1:*** Convertir todos los elementos $a_{i1}=0$ seleccionando un multiplicador de forma adecuada, en este caso $M_{i1}=a_{i1}/a_{11}$, y calcular la nueva $i–\\text{ésima}$ fila como: $F_i^{\\text{new}}=F_i - M_{i1}F_1$, con $2 \\leq i \\leq n$. Para todo $1 \\leq j \\leq n+1$ se cumple: $a_{ij}^{(new)}=a_{ij}-M_{i1}a_{1j}$.\n",
    "\n",
    "\n",
    "- ***Eliminación Etapa 2:*** Convertir todos los elementos $a_{i2}=0$ seleccionando un multiplicador de forma adecuada, en este caso $M_{i2}=a_{i2}/a_{22}$, y calcular la nueva $i–\\text{ésima}$ fila como: $F_i^{\\text{new}}=F_i - M_{i2}F_2$, con $3 \\leq i \\leq n$. Para todo $2 \\leq j \\leq n+1$ se cumple: $a_{ij}^{(new)}=a_{ij}-M_{i2}a_{2j}$.\n",
    "\n",
    "$$\\vdots$$\n",
    "\n",
    "\n",
    "- ***Eliminación Etapa k:*** Convertir todos los elementos $a_{ik}=0$ seleccionando un multiplicador de forma adecuada, en este caso $M_{ik}=a_{ik}/a_{kk}$, y calcular la nueva $i–\\text{ésima}$ fila como: $F_i^{\\text{new}}=F_i - M_{ik}F_k$, con $1 \\leq k \\leq n-1$ y $k-1 \\leq i \\leq n$. Para todo $k \\leq j \\leq n+1$ se cumple: $a_{ij}^{(new)}=a_{ij}-M_{ik}a_{kj}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Algoritmo proceso de Eliminación"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El presente seudocódigo muestra el algoritmo en el proceso de eliminación, empleando la matriz amentada $[\\mathbf{A}|\\mathbf{b}]$. \n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img04_Eliminacion.PNG?raw=true\" width=\"500\" />\n",
    "</p>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Resultado del proceso de Eliminación"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El resultado del proceso de eliminación sobre la matriz aumentada es el siguiente\n",
    "\n",
    "\n",
    "\\begin{align*}\n",
    "\\left[\\begin{array}{ccccccc}\n",
    "  a_{11} & a_{12} & a_{13} & \\ldots & a_{1n} & \\vdots & a_{1,n+1} \\\\\n",
    "  0      & a_{22} & a_{23} & \\ldots & a_{2n} & \\vdots & a_{2,n+1} \\\\\n",
    "  0      & 0      & a_{33} & \\ldots & a_{3n} & \\vdots & a_{3,n+1} \\\\\n",
    "  \\vdots & \\vdots & \\vdots & \\ddots & \\vdots & \\vdots & \\vdots     \\\\\n",
    "  0      & 0      & 0      & \\vdots & a_{nn} & \\vdots & a_{n,n+1} \\\\\n",
    "\\end{array}\\right]\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='Algoritmoeliminacion'></a>\n",
    "#### Implementación computacional proceso de Eliminación"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A continuación se presenta una implementación en `Python` del algoritmo de eliminación"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def eliminacion(Ab):\n",
    "    \"\"\"\n",
    "    función en Python 3.8 que realiza el proceso de eliminación (encontrar la triangular superior)\n",
    "    en un Sistema de Ecuaciones Lineales (SEL)\n",
    "    \n",
    "    input: Se ingresa la matriz ampliada Ab[n x n+1]\n",
    "    \n",
    "    output: devuelve la matriz ampliada Ab[n x n+1] triangular superior\n",
    "    \"\"\"\n",
    "    n = Ab.shape[0]\n",
    "    for k in range(0,n-1):\n",
    "        print(\"\\n\\nEtapa\", k+1, \"Eliminación columna: \", k+1, \"\\n\")\n",
    "        for i in range(k+1, n):\n",
    "            multiplicador = Ab[i][k]/Ab[k][k]\n",
    "            for j in range(k, n+1):\n",
    "                Ab[i][j] = Ab[i][j] - multiplicador * Ab[k][j]\n",
    "            print(\"\\n\",Ab)\n",
    "    return Ab\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sustitución regresiva (o despeje)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una vez se obtiene la matriz triangular superior, se procede a realizar el procedimiento de sustitución regresiva, que consiste en despejar la última incógnita de la última ecuación (observe que es una ecuación con una sola incógnita) quedando\n",
    "\n",
    "$$x_n = \\frac{a_{n,n+1}}{a_{n,n}}$$\n",
    "\n",
    "después se pasa a la penúltima ecuación y se despeja la penúltima incógnita, que también es una ecuación con una sola incógnita, ya que se conoce el valor de $x_n$ dado en el paso anterior, quedando:    \n",
    "\n",
    "$$x_{n-1} = \\frac{a_{n-1,n+1}-a_{n,n+1} \\times x_n}{a_{n,n}}$$\n",
    "\n",
    "se sigue este procedimiento hasta que se llega a la primera ecuación, resulta la siguiente fórmula recursiva:\n",
    "\n",
    "$$x_i=\\frac{a_{i,n+1}-\\sum\\limits_{j=i+1}^{n}a_{ij}x_j}{a_{ii}}, i=n, n-1, n-2, \\ldots, 1.$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Algoritmo proceso de sustitución regresiva"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img05_SustRegresiva.PNG?raw=true\" width=\"500\" />\n",
    "</p>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Implementación computacional"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A continuación se presenta una implementación en `Python` del algoritmo de sustitución regresiva"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sustitucionReg(Ab):\n",
    "    \"\"\"\n",
    "    función en Python 3.8 que realiza el proceso de sustitución regresiva\n",
    "    \n",
    "    input: Se ingresa la matriz triangular superior ampliada Ab[n * n+1]\n",
    "    \n",
    "    output: devuelve el vector solución x[n * 1]\n",
    "    \"\"\"    \n",
    "    \n",
    "    n = Ab.shape[0]\n",
    "    x = np.zeros(n)\n",
    "    x[n-1] = Ab[n-1][n] / Ab[n-1][n-1]\n",
    "    for i in range(n-2,-1,-1):\n",
    "        suma = 0\n",
    "        for j in range(i+1, n):\n",
    "            suma = suma + Ab[i][j] * x[j]\n",
    "        x[i] = (Ab[i][n] - suma) / Ab[i][i]\n",
    "    [print(\"x[{0}] = {1:6.4f}\".format(i,x[i])) for i in range(n)]\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplo teórico algoritmo Eliminación de Gauss simple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A manera de ejemplo paso-a-paso, vamos a construir un sistema de $5$ ecuaciones con $5$ incógnitas:\n",
    "\n",
    "\n",
    "$$a_{11}x_1 + a_{12}x_2 + a_{13}x_3 + a_{14}x_4 + a_{15}x_5 = b_1 \\hspace{1 cm} (ec. 1)$$\n",
    "\n",
    "$$a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + a_{24}x_4 + a_{25}x_5 = b_2 \\hspace{1 cm} (ec. 2)$$\n",
    "\n",
    "$$a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + a_{34}x_4 + a_{35}x_5 = b_3 \\hspace{1 cm} (ec. 3)$$\n",
    "\n",
    "$$a_{41}x_1 + a_{42}x_2 + a_{43}x_3 + a_{44}x_4 + a_{45}x_5 = b_4 \\hspace{1 cm} (ec. 4)$$\n",
    "\n",
    "$$a_{51}x_1 + a_{52}x_2 + a_{53}x_3 + a_{54}x_4 + a_{55}x_5 = b_5 \\hspace{1 cm} (ec. 5)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ***Eliminación:*** Dado el sistema $\\mathbf{Ax=b}$ se realizan operaciones fundamentales de fila hasta obtener un sistema triangular superior equivalente $\\mathbf{Ux=B}$.\n",
    "\n",
    "> ***Etapa $1$:*** Convertir todos los elementos $a_{i1}=0$ seleccionando un multiplicador de forma adecuada, en este caso $M_{i1}=a_{i1}/a_{11}$ , y calcular la nueva $i$–ésima fila como: $F_i^{(nueva)}=F_{i}−M_{i1} \\times F_1$, con $2\\leq i \\leq n$. \n",
    "\n",
    "\n",
    "> ***Etapa $2$:*** Convertir todos los elementos $a_{i2}=0$ seleccionando un multiplicador de forma adecuada, en este caso $M_{i2}=a_{i2}/a_{22}$ , y calcular la nueva $i$–ésima fila como: $F_i^{(nueva)}=F_{i}−M_{i2} \\times F_2$, con $3\\leq i \\leq n$.\n",
    "\n",
    "\n",
    "> ***Etapa $3$:*** Convertir todos los elementos $a_{i3}=0$ seleccionando un multiplicador de forma adecuada, en este caso $M_{i3}=a_{i3}/a_{33}$ , y calcular la nueva $i$–ésima fila como: $F_i^{(nueva)}=F_{i}−M_{i3} \\times F_3$, con $4\\leq i \\leq n$.\n",
    "\n",
    "\n",
    "> $\\vdots$\n",
    "\n",
    "resultando"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{matrix} \n",
    "a_{11}x_1 &+& a_{12}x_2 &+& a_{13}x_3 &+& a_{14}x_4 &+& a_{15}x_5 &=& b_1 \\\\ \n",
    "          & & a_{22}^{(1)}x_2 &+& a_{23}^{(1)}x_3 &+& a_{24}^{(1)}x_4 &+& a_{25}^{(1)}x_5 &=& b_2^{(1)} \\\\\n",
    "          & &           & & a_{33}^{(2)}x_3 &+& a_{34}^{(2)}x_4 &+& a_{35}^{(2)}x_5 &=& b_3^{(2)} \\\\\n",
    "          & &           & &           & & a_{44}^{(3)}x_4 &+& a_{45}^{(3)}x_5 &=& b_4^{(3)} \\\\\n",
    "          & &           & &           & &           & & a_{55}^{(4)}x_5 &=& b_5^{(4)} \\\\\n",
    "  \\end{matrix}$$\n",
    "\n",
    "Los superíndices entre paréntesis indican la cantidad de veecs que fue transformado ese término, debido a las operaciones de fila."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ***Sustitución Regresiva*** del sistema triangular superior: se halla primero el valor de $𝑥_5$, con este valor se halla $𝑥_4$, con estos dos valores se halla $𝑥_3$ y así sucesivamente hasta obtener el valor de $𝑥_1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{matrix} \n",
    "x_5 &=& \\frac{b_5^{(4)}}{a_{55}^{(4)}}\\\\\n",
    "x_4 &=& \\frac{b_4^{(3)}-a_{45}^{(3)}x_5}{a_{44}^{(3)}} \\\\\n",
    "x_3 &=& \\frac{b_3^{(2)}-a_{34}^{(2)}x_4 - a_{35}^{(2)}x_5}{a_{33}^{(2)}} \\\\\n",
    "x_2 &=& \\frac{b_2^{(1)}-a_{23}^{(1)}x_3 - a_{24}^{(1)}x_4 - a_{25}^{(1)}x_5 }{a_{22}^{(1)}} \\\\\n",
    "x_1 &=& \\frac{b_1 - a_{12}x_2 - a_{13}x_3 - a_{14}x_4 - a_{15}x_5 }{a_{11}} \\\\ \n",
    "  \\end{matrix}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='EjemploNumerico'></a>\n",
    "#### Ejemplo numérico algoritmo Eliminación de Gauss Simple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resolver el sistema $4 \\times 4$:\n",
    "\n",
    "$$6x_1  - 2x_2  + 2x_3 +  4x_4 = 12$$\n",
    "\n",
    "$$12x_1 - 8x_2  + 6x_3 + 10x_4 = 34$$\n",
    "\n",
    "$$3x_1  - 13x_2 + 9x_3 + 3x_4  = 27$$\n",
    "\n",
    "$$-6x_1 + 4x_2  +  x_3 - 18x_4 = -38$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='Multiplicadores'></a>\n",
    "- ***Eliminación:***\n",
    "> ***Etapa $1$:*** Convertir en ceros todos los elementos de la columna $1$ debajo del elemento pivote $(a_{11}=6)$, seleccionando un multiplicador para cada fila, de forma adecuada. En este caso:\n",
    "\n",
    "$$M_{21}=a_{21}/a_{11}=12/6 = 2$$\n",
    "\n",
    "$$M_{31}=a_{31}/a_{11}=3/6 = 1/2$$\n",
    "\n",
    "$$M_{41}=a_{41}/a_{11}=-6/6 = -1$$\n",
    "\n",
    "\n",
    "y calcular:\n",
    "\n",
    "> la nueva $2$a fila como: $F_2^{(nueva)}=F_{2}− M_{21} F_1 = F_{2}− 2 F_1$. \n",
    "\n",
    "> la nueva $3$a fila como: $F_3^{(nueva)}=F_{3}− M_{31} F_1 = F_{3}− F_1/2$. \n",
    "\n",
    "> la nueva $4$a fila como: $F_4^{(nueva)}=F_{4}- M_{41} F_1 = F_{4}- (-1)F_1 = F_{4} + F_1$. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{matrix} \n",
    "6x_1 &-&  2x_2 &+& 2x_3 &+&  4x_4 &=& 12 \\\\ \n",
    "  0  &-&  4x_2 &+& 2x_3 &+&  2x_4 &=& 10 \\\\\n",
    "  0  &-& 12x_2 &+& 8x_3 &+&   x_4 &=& 21 \\\\\n",
    "  0  & &  2x_2 &+& 3x_3 &-& 14x_4 &=& -26\\\\\n",
    "  \\end{matrix}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ***Etapa $2$:*** Convertir en ceros todos los elementos de la columna $2$ debajo del elemento pivote $(a_{22}=-4)$, seleccionando un multiplicador para cada fila, de forma adecuada. En este caso:\n",
    "\n",
    "$$M_{32}=a_{32}/a_{22}=-12/-4 = 3$$\n",
    "\n",
    "$$M_{42}=a_{42}/a_{22}= 2/-4 = -1/2$$\n",
    "\n",
    "y calcular \n",
    "\n",
    "> la nueva $3$a fila como: $F_3^{(nueva)}=F_{3}−3 F_2$. \n",
    "\n",
    "> la nueva $4$a fila como: $F_4^{(nueva)}=F_{4}+ F_2/2$. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{matrix} \n",
    "6x_1 &-&  2x_2 &+& 2x_3 &+&  4x_4 &=& 12 \\\\ \n",
    "  0  &-&  4x_2 &+& 2x_3 &+&  2x_4 &=& 10 \\\\\n",
    "  0  & &    0  & & 2x_3 &-&  5x_4 &=& -9 \\\\\n",
    "  0  & &    0  & & 4x_3 &-& 13x_4 &=& -21\\\\\n",
    "  \\end{matrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> ***Etapa $3$:*** Convertir en ceros todos los elementos de la columna $3$ debajo del elemento pivote $(a_{33}=2)$, seleccionando un multiplicador de forma adecuada. En este caso:\n",
    "\n",
    "$$M_{43}=a_{43}/a_{33}=4/2 = 2$$\n",
    "\n",
    "y calcular la nueva $4$a fila como: $F_4^{(nueva)}=F_{4}−2 F_3$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='MEliminada'></a>\n",
    "$$\\begin{matrix} \n",
    "6x_1 &-&  2x_2 &+& 2x_3 &+&  4x_4 &=& 12 \\\\ \n",
    "  0  &-&  4x_2 &+& 2x_3 &+&  2x_4 &=& 10 \\\\\n",
    "  0  & &    0  & & 2x_3 &-&  5x_4 &=& -9 \\\\\n",
    "  0  & &    0  & &  0   &-&  3x_4 &=& -3 \\\\\n",
    "  \\end{matrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se observa que el sistema resultante es triangular superior.\n",
    "\n",
    "Por último, aplicando el algoritmo de la *sustitución regresiva*:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{matrix} \n",
    "x_4 &=& \\frac{-3}{-3}= 1\\\\\n",
    "x_3 &=& \\frac{-9+5 \\times 1}{2} = -2 \\\\\n",
    "x_2 &=& \\frac{10-2(2)-2(1)}{-4} = -3\\\\\n",
    "x_1 &=& \\frac{12+2(-3)-2(-2)-4(1)}{6}=1\\\\\n",
    "\\end{matrix}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplo computacional algoritmo de Eliminación de Gauss Simple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vamos a resolver el problema usando programación en el lenguaje `Python`. Se importarán algunas bibliotecas matemáticas para trabajar con arreglos en vez de listas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "usando los algoritmos vistos:\n",
    "\n",
    "- Lea la matriz de coeficientes $\\mathbf{A}$ y el vector de términos independientes $\\mathbf{b}$\n",
    "\n",
    "\n",
    "- Realice el proceso de eliminación para obtener la matríz triangular superior\n",
    "\n",
    "\n",
    "- Realice el proceso de sustitución regresiva para obtener el vector solución\n",
    "\n",
    "Generando la matriz de coeficientes $\\mathbf{A}$ y el vector de términos independientes $\\mathbf{b}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  6.  -2.   2.   4.]\n",
      " [ 12.  -8.   6.  10.]\n",
      " [  3. -13.   9.   3.]\n",
      " [ -6.   4.   1. -18.]]\n"
     ]
    }
   ],
   "source": [
    "#Ingreso de la matriz de coeficientes A, de forma manual\n",
    "\n",
    "A = np.array([[ 6., - 2, 2,   4],\n",
    "              [12, - 8, 6,  10],\n",
    "              [ 3, -13, 9,   3],\n",
    "              [-6,   4, 1, -18]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 12.  34.  27. -38.]\n"
     ]
    }
   ],
   "source": [
    "# Ingreso del vector de términos independientes b, de forma manual\n",
    "b = np.array([12., 34, 27, -38])\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Creando la matriza aumentada $[\\mathbf{A}|\\mathbf{b}]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  6.  -2.   2.   4.  12.]\n",
      " [ 12.  -8.   6.  10.  34.]\n",
      " [  3. -13.   9.   3.  27.]\n",
      " [ -6.   4.   1. -18. -38.]]\n",
      "(4, 5)\n"
     ]
    }
   ],
   "source": [
    "Ab = np.c_[A,b]\n",
    "print(Ab)\n",
    "print(np.shape(Ab))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ***Proceso de Eliminación***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Etapa 1 Eliminación columna:  1 \n",
      "\n",
      "\n",
      " [[  6.  -2.   2.   4.  12.]\n",
      " [  0.  -4.   2.   2.  10.]\n",
      " [  3. -13.   9.   3.  27.]\n",
      " [ -6.   4.   1. -18. -38.]]\n",
      "\n",
      " [[  6.  -2.   2.   4.  12.]\n",
      " [  0.  -4.   2.   2.  10.]\n",
      " [  0. -12.   8.   1.  21.]\n",
      " [ -6.   4.   1. -18. -38.]]\n",
      "\n",
      " [[  6.  -2.   2.   4.  12.]\n",
      " [  0.  -4.   2.   2.  10.]\n",
      " [  0. -12.   8.   1.  21.]\n",
      " [  0.   2.   3. -14. -26.]]\n",
      "\n",
      "\n",
      "Etapa 2 Eliminación columna:  2 \n",
      "\n",
      "\n",
      " [[  6.  -2.   2.   4.  12.]\n",
      " [  0.  -4.   2.   2.  10.]\n",
      " [  0.   0.   2.  -5.  -9.]\n",
      " [  0.   2.   3. -14. -26.]]\n",
      "\n",
      " [[  6.  -2.   2.   4.  12.]\n",
      " [  0.  -4.   2.   2.  10.]\n",
      " [  0.   0.   2.  -5.  -9.]\n",
      " [  0.   0.   4. -13. -21.]]\n",
      "\n",
      "\n",
      "Etapa 3 Eliminación columna:  3 \n",
      "\n",
      "\n",
      " [[ 6. -2.  2.  4. 12.]\n",
      " [ 0. -4.  2.  2. 10.]\n",
      " [ 0.  0.  2. -5. -9.]\n",
      " [ 0.  0.  0. -3. -3.]]\n"
     ]
    }
   ],
   "source": [
    "#timein = np.time.time()\n",
    "Ab = eliminacion(Ab)\n",
    "#timeout = np.time.time()\n",
    "#elepsedtimeE = timeout - timein"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ***Proceso de Sustitución regresiva***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x[0] = 1.0000\n",
      "x[1] = -3.0000\n",
      "x[2] = -2.0000\n",
      "x[3] = 1.0000\n"
     ]
    }
   ],
   "source": [
    "#timein = np.time.time()\n",
    "x = sustitucionReg(Ab)\n",
    "#timeout = np.time.time()\n",
    "#elepsedtimeS = timeout - timein\n",
    "\n",
    "#totalT = elepsedtimeS + elepsedtimeE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resolviendo el *SEL* por el método que incorpora `Python` por defecto:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1. -3. -2.  1.]\n"
     ]
    }
   ],
   "source": [
    "x = np.linalg.solve(A,b)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comprobando que la solución es correcta, multiplicando el vector solución por la matriz de coeficientes para obtener el vector de términos independientes, es decir:\n",
    "\n",
    "$$\\mathbf{Ax=b}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([False,  True, False,  True])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.dot(x) == b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 12.  34.  27. -38.]\n"
     ]
    }
   ],
   "source": [
    "print(A.dot(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 12.  34.  27. -38.]\n"
     ]
    }
   ],
   "source": [
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Método de Gauss - Jordan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "El método de eliminación de [*Gauss - Jordan*](https://es.wikipedia.org/wiki/Eliminaci%C3%B3n_de_Gauss-Jordan) consiste básicamente en continuar con el proceso de eliminación dado en el algoritmo de eliminación simple de Gauss, hasta obtener una matriz diagonal (únicamente tiene elementos diferentes de cero en la diagonal principal) o identidad (únicamente contiene los elementos de la diagonal principal diferentes a 0 e iguales a 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{matrix} \n",
    "a_{11}^{(7)}x_1 & &    0       & &   0       & &  0         &=& b_{1}^{(7)}  \\\\ \n",
    "  0       & &  a_{22}^{(6)}x_2 & &   0       & &  0         &=& b_{2}^{(6)} \\\\\n",
    "  0       & &    0       & & a_{33}^{(5)}x_3 & &  0         &=& b_{3}^{(5)} \\\\\n",
    "  0       & &    0       & &   0       & &  a_{44}^{(4)}x_4 &=& b_{4}^{(4)} \\\\\n",
    "  \\end{matrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La solución del sistema es \"inmediata\":"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{matrix} \n",
    "x_1 &=& \\frac{b_{1}^{(7)}}{a_{11}^{(7)}}  \\\\ \n",
    "x_2 &=& \\frac{b_{2}^{(7)}}{a_{22}^{(6)}}  \\\\ \n",
    "x_3 &=& \\frac{b_{3}^{(7)}}{a_{33}^{(5)}}  \\\\ \n",
    "x_4 &=& \\frac{b_{4}^{(7)}}{a_{44}^{(4)}}  \\\\ \n",
    "\\end{matrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Continuando con el ejemplo anterior:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{matrix} \n",
    "6x_1 &-&  2x_2 &+& 2x_3 &+&  4x_4 &=& 12 \\\\ \n",
    "  0  &-&  4x_2 &+& 2x_3 &+&  2x_4 &=& 10 \\\\\n",
    "  0  & &    0  & & 2x_3 &-&  5x_4 &=& -9 \\\\\n",
    "  0  & &    0  & &  0   &-&  3x_4 &=& -3 \\\\\n",
    "  \\end{matrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se emplea el mismo algoritmo que se venía utilizando, pero ahora se empezará desde la columna $4$, fila $1$, y se eliminan todos los elementos hasta el elemento que está encima del elemento pivote. Una vez eliminados todos los elementos encima del pivote, se continua con la columna $3$ y así sucesivamente hasta la columna $2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- <strong>*Etapa $4$*:</strong> Convertir en ceros todos los elementos de la columna $4$ encima del elemento pivote $(a_{44}=-3)$, seleccionando un multiplicador para cada fila, de forma adecuada. En este caso:\n",
    "\n",
    "> $M_{14} = a_{14}/a_{44} =  4/(-3)$\n",
    "\n",
    "> $M_{24} = a_{24}/a_{44} =  2/(-3)$ \n",
    "\n",
    "> $M_{34} = a_{34}/a_{44} = -5/(-3)$ \n",
    "\n",
    "\n",
    "y calcular:\n",
    "\n",
    "> la nueva $1$a fila como: $F_1^{(nueva)}=F_{1}− M_{14} F_4 = F_{1} + 4/3 F_4$. \n",
    "\n",
    "> la nueva $2$a fila como: $F_2^{(nueva)}=F_{2}− M_{24} F_4 = F_{2} + 2/3 F_4$. \n",
    "\n",
    "> la nueva $3$a fila como: $F_3^{(nueva)}=F_{3}- M_{34} F_4 = F_{3} - 5/3 F_4 $. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{matrix} \n",
    "6x_1 &-&  2x_2 &+& 2x_3 & &     0 &=&  8 \\\\ \n",
    "  0  &-&  4x_2 &+& 2x_3 & &     0 &=&  8 \\\\\n",
    "  0  & &    0  & & 2x_3 & &     0 &=& -4 \\\\\n",
    "  0  & &    0  & &  0   &-&  3x_4 &=& -3 \\\\\n",
    "  \\end{matrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- <strong>*Etapa $5$*:</strong> Convertir en ceros todos los elementos de la columna $3$ encima del elemento pivote $(a_{33}=2)$, seleccionando un multiplicador para cada fila, de forma adecuada. En este caso:\n",
    "\n",
    "> $M_{13} = a_{13}/a_{33} = 2/2 = 1$ \n",
    "\n",
    "> $M_{23} = a_{23}/a_{33} = 2/2 = 1$ \n",
    "\n",
    "y calcular \n",
    "\n",
    "> la nueva $1$a fila como: $F_1^{(nueva)} = F_{1} − F_3$. \n",
    "\n",
    "> la nueva $2$a fila como: $F_2^{(nueva)} = F_{2} - F_3$. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{matrix} \n",
    "6x_1 &-&  2x_2 & & 0    & &     0 &=&  12 \\\\ \n",
    "  0  &-&  4x_2 & & 0    & &     0 &=&  12 \\\\\n",
    "  0  & &    0  & & 2x_3 & &     0 &=& -4 \\\\\n",
    "  0  & &    0  & &  0   &-&  3x_4 &=& -3 \\\\\n",
    "  \\end{matrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- <strong>*Etapa $6$*:</strong> Convertir en ceros todos los elementos de la columna $2$ encima del elemento pivote $(a_{22}=-4)$, seleccionando un multiplicador de forma adecuada. En este caso:\n",
    "\n",
    "> $M_{12} = a_{12}/a_{22} = -2/-4 = 1/2$ \n",
    "\n",
    "y calcular la nueva $1$a fila como: $F_1^{(nueva)}=F_{1}− F_2/2$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{matrix} \n",
    "6x_1 & &     0 & & 0    & &     0 &=&  6 \\\\ \n",
    "  0  &-&  4x_2 & & 0    & &     0 &=&  12 \\\\\n",
    "  0  & &    0  & & 2x_3 & &     0 &=& -4 \\\\\n",
    "  0  & &    0  & &  0   &-&  3x_4 &=& -3 \\\\\n",
    "  \\end{matrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La solución es inmediata:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{matrix} \n",
    "x_1 &=& \\frac{6}{6} &=& 1 \\\\ \n",
    "x_2 &=& \\frac{12}{-4} &=& -3  \\\\ \n",
    "x_3 &=& \\frac{-4}{-2} &=& -2  \\\\ \n",
    "x_4 &=& \\frac{-3}{-3} &=& 1  \\\\ \n",
    "\\end{matrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert alert-success\">\n",
    "$\\color{red}{\\textbf{Laboratorio computacional:}}$ \n",
    "    \n",
    "- Implemente el algoritmo del método de eliminación de Gauss - Jordan (Emplee el algoritmo de eliminación dado como punto de partida). \n",
    "    \n",
    "- Realice un análisis comparativo de la cantidad de operaciones necesarias en cada uno de ellos. ¿Cuál considera qué es mejor? ¿Por qué? ¿Qué criterios tuvo en cuenta para dicho análisis.\n",
    "    \n",
    "- ¿En qué casos emplearía un método u otro? ¿Por qué?\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Aspectos computacionales, inconvenitentes de los métodos de eliminación"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Análisis de desempeño algoritmo de Eliminación de Gauss Simple"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La cantidad de operaciones (sumas, restas, multiplicaciones y/o divisiones) que hay qué realizar para alcanzar el sistema triangular superior (eliminación) de $n$ ecuaciones, es de\n",
    "\n",
    "$$\\frac{3(n-1)n}{2}+\\frac{(n-1)n(2n-1)}{3}$$\n",
    "\n",
    "y para efectuar el proceso de sustitución regresiva (solución del sistema) se requieren de $n^2$ operaciones aritméticas. \n",
    "\n",
    "En total, para resolver un *SEL* empleando el método de Eliminación de Gauss simple (sin contabilizar las operaciones de pivoteo) es de:\n",
    "\n",
    "$$\\frac{4n^3+9n^2-7n}{6}$$\n",
    "\n",
    "Esta expresión se puede escribir de forma mas compacta como:\n",
    "\n",
    "$$ \\frac{2n^3}{3}+O(n^2)$$\n",
    "\n",
    "La anterior ecuación indica que el número de [operaciones de punto flotante por segundo](https://en.wikipedia.org/wiki/FLOPS \"flops\") (flops, en inglés) es del orden de $n^3$ más una componente de orden $n^2$ y menores. Es decir, para resolver un sistema de 10 ecuaciones con 10 incógnitas por el método de Eliminacion de Gauss se requieren del orden de 800 operaciones aritméticas. En la siguiente tabla se observa un resumen de las magnitudes de operaciones dependiendo del orden de $n$.\n",
    "\n",
    "|  n   | Eliminación | Sust. atrás | flops | 2$n^3$/3 | % elim. |\n",
    "|-----:|------------:|------------:|-----:|----------------:|-------:|\n",
    "| 10   |      375    |     55      |   430 |    333           |  87.21% |\n",
    "| 100  |   338250    |   5050      | 343300 |   333333           |  98.53% |\n",
    "| 1000 | 3.34$\\times 10^8$ | 500500| 3.34 $\\times 10^8$ | 3.33$\\times 10^8$|  99.85% |\n",
    "\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": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### División entre cero"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se presenta cuando un elemento de la diagonal principal es cero, $a_{ii}=0$ ($i=1,2,3,\\ldots, n$). Se evita realizando un procedimiento de pivoteo, parcial o total, de filas y/o columnas. Esta técnica se verá en otro apartado más adelante."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Errores de redondeo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observando la respuesta del ejemplo anterior y haciendo la comparación $\\mathbf{Ax=b}$, se observa que uno puede llegar a pensar que efectivamente la solución encontrada es la solución del sistema en custión. Pero realizando una comprobación computacional, evaluando la igualdad dada por el sistema, se tiene:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.dot(x) == b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(A.dot(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se observa que hay dos componentes que no cumplen con la igualdad."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bGS = A.dot(x)\n",
    "[print(\"bGS[{0}] = {1:6.20f}\".format(i,bGS[i])) for i in range(len(x))];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Error relativo %\n",
    "[print(\"error[{0}] = {1:6.20f}\".format(i,abs(bGS[i]-b[i])/abs(b[i])*100)) for i in range(len(x))];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La discrepancia, aunque pequeña en este simple ejemplo, es debida al redondeo realizado en las diferentes operaciones aplicadas. Si se hubiera trabajado con fraccionarios en lugar de decimales, evitando los errores de redondeo, la respuesta hubiera sido exacta."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Condicionamiento del sistema"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Los sistemas bien condicionados son aquellos en los que un pequeño cambio en alguno de los coeficientes de la matriz $\\mathbf{A}$ provoca un cambio pequeño en la solución. Lo contrario serían los sistemas mal condicionados. Una interpretación alternativa puede ser que un amplio rango de resultados puede satisfacer las ecuaciones en forma aproximada, como es el caso cuando se emplean computadores para la solución.\n",
    "\n",
    "- ***Ejemplo:*** Resolver el siguiente SEL\n",
    "\n",
    "$$x_1 + 2x_2 = 10 $$\n",
    "$$1.1 x_1 + 2x_2 = 10.4 $$\n",
    "\n",
    "La solución es inmediata:\n",
    "\n",
    "$$x_1=\\frac{2(10)-2(10.4)}{1(2)-2(1.1)}=4$$\n",
    "$$x_2=\\frac{1(10.4)-1.1(10)}{1(2)-2(1.1)}=3$$\n",
    "\n",
    "Ahora realizando un pequeño cambio en el coeficiente $a_{21}$ de $1.1$ a $1.05$, el resultado es:\n",
    "\n",
    "$$x_1=\\frac{2(10)-2(10.4)}{1(2)-2(1.0.5)}=8$$\n",
    "$$x_2=\\frac{1(10.4)-1.1(10)}{1(2)-2(1.0.5)}=1$$\n",
    "\n",
    "La diferencia en los resultados se da porque en el denominador se tiene la resta de dos números casi iguales, este tipo de diferencia es muy sensible a pequeños cambios en las cantidades empleadas.\n",
    "\n",
    "En un sistema mal condicionado, las pendientes de las rectas son casi iguales, y visualmente es difícil percibir el punto exacto donde se cruzan. Revisemos matemáticamente esta situación. Escribiendo un sistema $2\\times2$ de forma general, se tiene:\n",
    "\n",
    "$$a_{11}x_1 + a_{12}x_2 = b_1 $$\n",
    "$$a_{21}x_1 + a_{22}x_2 = b_2 $$\n",
    "\n",
    "Reescribiendo estas ecuaciones en la forma alternativa:\n",
    "\n",
    "$$x_2=-\\frac{a_{11}}{a_{12}}x_1+\\frac{b_1}{a_{12}}$$\n",
    "$$x_2=-\\frac{a_{21}}{a_{22}}x_1+\\frac{b_2}{a_{22}}$$\n",
    "\n",
    "y si las pendientes son casi iguales,\n",
    "\n",
    "$$\\frac{a_{11}}{a_{12}}\\approx \\frac{a_{21}}{a_{22}}$$\n",
    "\n",
    "llegando finalmente a\n",
    "\n",
    "$$a_{11}a_{22}-a_{12}a_{21}\\approx 0$$\n",
    "\n",
    "Recordando del curso de Algebra Lineal, el [determinante](https://en.wikipedia.org/wiki/Determinant) de un sistema $2\\times2$ es justamente $a_{11}a_{22}-a_{12}a_{21}$, por lo que se concluye que ***un sistema mal condicionado es aquel cuyo determinante es cercano a cero***."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='Escalamiento'></a>\n",
    "#### Escalamiento "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La gran pregunta a ser resuelta del ítem anterior es: \"***qué es estar cerca al cero***\". Esto es más complicado aún porque el valor del determinante puede variar al multiplicar las ecuaciones involucradas por algún factor de escalamiento sin alterar la solución del sistema. \n",
    "\n",
    "- ***Ejemplo:***\n",
    "\n",
    "Tomando como ejemplo el mismo sistema $2\\times2$ del ejemplo anterior, se encuentra que su determinante es:\n",
    "\n",
    "$$D=|1(2)-2(1.1)|=0.2$$\n",
    "\n",
    "Se corrobora que su determinante al ser relativamente cercano a cero se tiene un sistema mal condicionado. Ahora si multiplicamos todo el sistema por un factor de $10$, se tiene:\n",
    "\n",
    "$$10x_1 + 20x_2 = 100 $$\n",
    "$$11x_1 + 20x_2 = 104 $$\n",
    "\n",
    "y su determinante es:\n",
    "\n",
    "$$D=|10(20)-20(11)|=20 \\gg 0$$\n",
    "\n",
    "que efectivamente podemos afirmar es mucho mayor que cero, pero aún persiste el problema del mal condicionamiento, pues la solución es la misma."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Técnicas para mejorar la solución en los métodos de Eliminación"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Escalamiento"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como se observó en el numeral anterior, el efecto del escalamiento esconde lo que realmente puede suceder en un sistema mal condicionado. Para evitar dicho inconveniente, se puede escalar cada una de las ecuaciones de tal manera que el máximo valor de cualquiera de los coeficientes de dicha ecuación sea igual a $1$. Esto se hace simplemente dividiendo por el valor del mayor coeficiente presente en cada ecuación.\n",
    "\n",
    "- ***Ejemplo:***\n",
    "\n",
    "Realizando el escalamiento en el sistema de ecuaciones\n",
    "\n",
    "$$x_1 + 2x_2 = 10 $$\n",
    "$$1.1 x_1 + 2x_2 = 10.4 $$\n",
    "\n",
    "dividiendo por $2$ ambas ecuaciones, para que el mayor coeficiente sea $1.0$ en cada una de las dos ecuaciones, resultado en\n",
    "\n",
    "$$0.5x_1 + x_2 = 5 $$\n",
    "$$0.55 x_1 + x_2 = 5.2 $$\n",
    "\n",
    "cuyo determinante es \n",
    "\n",
    "$$D=|0.5(1)-1(0.55)|=0.05$$\n",
    "\n",
    "y ahora realizando el cálculo del determinante sobre las ecuaciones afectadas por un factor de $10$, se llega al mismo resultado: $0.05$ y se evidencia claramente la dificultad del mal condicionamiento del sistema.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pivoteo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El [pivote](https://en.wikipedia.org/wiki/Pivot_element \"pivote\") es el elemento de una matriz que se selecciona primero mediante un algoritmo (por ejemplo, eliminación gaussiana, algoritmo simplex, etc.), para realizar ciertos cálculos. En este tema, el elemento pivote será el elemento que se encuentra en la posición $a_{ii}$, correspondiente a la diagonal principal.\n",
    "\n",
    "Antes de empezar con el proceso de eliminación, es conveniente determinar si existe algún elemento igual, o cercano, a cero en los elementos de la diagonal principal ($a_{ii}\\approx 0$), ya que esto ocasionará grandes errores de redondeo, afectando por lo tanto la solución. En caso afirmativo, la solución a este problema es casi obvia: si el elemento de la diagonal principal, $a_{ii}=0$ en la etapa $k$, entonces intercambie la fila $k$ del sistema con alguna fila posterior cuya entrada $a_{ii} \\neq0$. Tal intercambio no altera la solución del sistema. Con una entrada diagonal distinta de cero como pivote, el proceso puede continuar como de costumbre.\n",
    "\n",
    "En principio, cualquier valor distinto de cero servirá como pivote para calcular los multiplicadores, pero en la práctica la elección debe hacerse con cierto cuidado para minimizar el error. Cuando la porción restante de la matriz se multiplica por la matriz de eliminación elemental resultante, debemos intentar limitar el crecimiento de las entradas de la matriz transformada para no amplificar los errores de redondeo. Por esta razón, es deseable que los multiplicadores no excedan de $1$ en magnitud. Este requisito se puede cumplir eligiendo la entrada de mayor magnitud en o debajo de la diagonal como pivote, tal como se explicó en el apartado de [escalamiento](#Escalamiento). Esta política se denomina pivoteo parcial y es esencial en la práctica para una implementación numéricamente estable de la eliminación gaussiana para sistemas lineales generales.\n",
    "\n",
    "- ***Ejemplo:*** Revise el ejemplo 9.9 propuesto en el libro de Chapra y Canale.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert alert-success\">\n",
    "$\\color{red}{\\textbf{Laboratorio computacional:}}$ \n",
    "    \n",
    "- Complemente los algoritmos computacionales vistos en la primera parte del capítulo con la rutina de Pivoteo. (***Ayuda:*** Ver la figura 9.6 del libro de Chapra y Canale, p. 272)\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Factorización LU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Introducción"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como se indicó al comienzo de este capítulo, el método de [Eliminación de Gauss simple](#EGS) se emplea para resolver sistemas de ecuaciones lineales de la forma $\\mathbf{Ax=b}$. A pesar de ser un método bastante simple para ser aplicado, no está exento de inconvenientes. A la ineficiencia en el desempeño computacional y posibilidades de errores numéricos debido al redondeo, se le suma lo poco práctico al momento de resolver sistemas en los que el vector de términos independientes $\\mathbf{b}$ es diferente en cada problema a desarrollar, manteniendo la matriz de coeficientes $\\mathbf{A}$ la misma. Este es un típico problema que se presenta cuando se quiere evaluar el comportamiento de un sistema físico al cambiar únicamente las condiciones de frontera y/o iniciales. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Determinación de la descomposición LU "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Este proceso de eliminación también es conocido como factorización, o descomposición, $\\mathbf{LU}$ porque descompone la matriz $\\mathbf{A}$ en un producto de una matriz triangular inferior unitaria, $\\mathbf{L}$, y una matriz triangular superior, $\\mathbf{U}$. \n",
    "\n",
    "En forma matricial, al sistema $\\mathbf{Ax=b}$ lo podemos reordenar como\n",
    "\n",
    "$$[\\mathbf{A}]\\{\\mathbf{x}\\}-\\{\\mathbf{b}\\}=\\{\\mathbf{0}\\}$$\n",
    "\n",
    "Ahora supongamos que podemos expresar el proceso de eliminación de la forma\n",
    "\n",
    "$$[\\mathbf{U}]\\{\\mathbf{x}\\}-\\{\\mathbf{d}\\}=\\{\\mathbf{0}\\}$$\n",
    "\n",
    "donde $\\mathbf{U}$ es una matriz triangular superior, quedando el sistema\n",
    "\n",
    "\\begin{align*}\n",
    "\\left[\\begin{array}{ccc}\n",
    "  u_{11} & u_{12} & u_{13} \\\\\n",
    "  0      & u_{22} & u_{23} \\\\\n",
    "  0      & 0      & u_{33} \\\\  \n",
    " \\end{array}\\right]\n",
    "\\begin{Bmatrix}\n",
    "  x_{1}  \\\\\n",
    "  x_{2}  \\\\\n",
    "  x_{3}  \\\\  \n",
    "\\end{Bmatrix}\n",
    "= \\begin{Bmatrix}\n",
    "  d_{1}  \\\\\n",
    "  d_{2}  \\\\\n",
    "  d_{3}  \\\\\n",
    "\\end{Bmatrix}\n",
    "\\end{align*}\n",
    "\n",
    "Igualmente supongamos que existe una matriz diagonal inferior con diagonal unitaria, $\\mathbf{L}$\n",
    "\n",
    "\\begin{align*}L=\n",
    "\\left[\\begin{array}{ccc}\n",
    "  1      & 0      & 0 \\\\\n",
    "  l_{21} & 1      & 0 \\\\\n",
    "  l_{31} & l_{32} & 1 \\\\  \n",
    " \\end{array}\\right]\n",
    "\\end{align*}\n",
    "\n",
    "cuya propiedad es que al posmultiplicarla por $[\\mathbf{U}]\\{\\mathbf{x}\\}-\\{\\mathbf{d}\\}$, obtenemos el sistema matricial inicia, es decir:\n",
    "\n",
    "$$[\\mathbf{L}]\\{[\\mathbf{U}]\\{\\mathbf{x}\\}-\\{\\mathbf{d}\\}\\}=[\\mathbf{A}]\\{\\mathbf{x}\\}-\\{\\mathbf{b}\\}$$\n",
    "\n",
    "de aquí se obtiene\n",
    "\n",
    "$$[\\mathbf{L}][\\mathbf{U}]=[\\mathbf{A}]$$\n",
    "\n",
    "y\n",
    "\n",
    "$$[\\mathbf{L}][\\mathbf{d}]=[\\mathbf{b}]$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Estrategia de solución"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El anterior procedimiento determina una estrategia de dos pasos para llegar a la solución:\n",
    "\n",
    "1. ***Descomposición $\\mathbf{LU}$:*** $\\mathbf{A}$ se factoriza en las matrices triangulares inferior $\\mathbf{L}$ y superior $\\mathbf{U}$.\n",
    "\n",
    "2. ***Sustitución:*** $\\mathbf{L}$ y $\\mathbf{U}$ se emplean para obtener $\\mathbf{x}$ para un valor del vector de términos independientes $\\mathbf{b}$.\n",
    "\n",
    "Resumiendo, veámos el siguiente gráfico:\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img06_LU.PNG?raw=true\" width=\"500\" />\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": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Obtención de las matrices *L* y *U*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La obtención de las matrices $\\mathbf{L}$ y $\\mathbf{U}$ es inmediata, resultantes de la descomposición realizada en el método de la [Eliminación de Gauss simple](#EGS). \n",
    "\n",
    "- $\\mathbf{U}$, es la misma resultante del proceso de eliminación, donde se obtiene la matriz triangular superior.\n",
    "\n",
    "onde $\\mathbf{U}$ es una matriz triangular superior, quedando el sistema\n",
    "\n",
    "\\begin{align*}U=\n",
    "\\left[\\begin{array}{ccc}\n",
    "  a_{11} & a_{12} & a_{13} \\\\\n",
    "  0      & a'_{22} & a'_{23} \\\\\n",
    "  0      & 0      & a''_{33} \\\\  \n",
    " \\end{array}\\right]\n",
    "\\end{align*}\n",
    "\n",
    "donde cada apóstrofe indica la cantidad de variaciones que ha sufrido cada coeficiente $a_{ij}$ durante el proceso de eliminación. \n",
    "\n",
    "- $\\mathbf{L}$, resulta de almacenar los multiplicadores $M_{ik}=a_{ik}/a_{kk}$ que se emplearon en las [etapas de eliminación](#Etapaseliminacion), con la diagonal principal unitaria (el subindice $k$ corresponde a la $k$-ésima etapa de eliminación, que coincide con el número de la columna en la matriz de coeficientes).\n",
    "\n",
    "\\begin{align*}L=\n",
    "\\left[\\begin{array}{ccc}\n",
    "  1      & 0      & 0 \\\\\n",
    "  M_{21} & 1      & 0 \\\\\n",
    "  M_{31} & M_{32} & 1 \\\\  \n",
    " \\end{array}\\right]\n",
    "\\end{align*}\n",
    "\n",
    "entonces,\n",
    "\n",
    "$$M_{21}=\\frac{a_{21}}{a_{11}}$$\n",
    "\n",
    "$$M_{31}=\\frac{a_{31}}{a_{11}}$$\n",
    "\n",
    "$$M_{32}=\\frac{a'_{32}}{a'_{22}}$$\n",
    "\n",
    "$$\\vdots$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplo *LU*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Del [Ejemplo numérico](#EjemploNumerico) de Eliminación de Gauss simple, se tiene el siguiente sistema de ecuaciones lineales\n",
    "\n",
    "$$6x_1  - 2x_2  + 2x_3 +  4x_4 = 12$$\n",
    "\n",
    "$$12x_1 - 8x_2  + 6x_3 + 10x_4 = 34$$\n",
    "\n",
    "$$3x_1  - 13x_2 + 9x_3 + 3x_4  = 27$$\n",
    "\n",
    "$$-6x_1 + 4x_2  +  x_3 - 18x_4 = -38$$\n",
    "\n",
    "\n",
    "la matriz $\\mathbf{U}$ es la [matriz de eliminación](#MEliminada), dada por:\n",
    "\n",
    "\\begin{align*}U=\n",
    "\\left[\\begin{array}{cccc}\n",
    "  6 & -2 & 2 &  4 \\\\\n",
    "  0 & -4 & 2 &  2 \\\\\n",
    "  0 &  0 & 2 & -5 \\\\\n",
    "  0 &  0 & 0 & -3 \\\\\n",
    "  \\end{array}\\right]\n",
    "\\end{align*}\n",
    "\n",
    "y la matriz $\\mathbf{L}$ está dada por cada uno de los elementos de los [multiplicadores](#Multiplicadores) $M_{ik}$:\n",
    "\n",
    "$$M_{21}=a_{21}/a_{11}=12/6 = 2$$\n",
    "\n",
    "$$M_{31}=a_{31}/a_{11}=3/6 = 1/2$$\n",
    "\n",
    "$$M_{41}=a_{41}/a_{11}=-6/6 = -1$$\n",
    "\n",
    "$$M_{32}=a_{32}/a_{22}=-12/-4 = 3$$\n",
    "\n",
    "$$M_{42}=a_{42}/a_{22}= 2/-4 = -1/2$$\n",
    "\n",
    "$$M_{43}=a_{43}/a_{33}=4/2 = 2$$\n",
    "\n",
    "\\begin{align*}L=\n",
    "\\left[\\begin{array}{cccc}\n",
    "    1  &  0 & 0 & 0 \\\\\n",
    "    2  & 1 & 0 & 0 \\\\\n",
    "   1/2 & 3 & 1 & 0 \\\\\n",
    "   -1  & -1/2 & 2 & 1 \\\\\n",
    "  \\end{array}\\right]\n",
    "\\end{align*}\n",
    "\n",
    "Verificando que $\\mathbf{LU=A}$\n",
    "\n",
    "\\begin{align*}\n",
    "\\begin{bmatrix}\n",
    "    1  &  0 & 0 & 0 \\\\\n",
    "    2  & 1 & 0 & 0 \\\\\n",
    "   1/2 & 3 & 1 & 0 \\\\\n",
    "   -1  & -1/2 & 2 & 1 \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "  6 & -2 & 2 &  4 \\\\\n",
    "  0 & -4 & 2 &  2 \\\\\n",
    "  0 &  0 & 2 & -5 \\\\\n",
    "  0 &  0 & 0 & -3 \\\\\n",
    "\\end{bmatrix}\n",
    "= \\begin{bmatrix}\n",
    "   6 & - 2 & 2 &   4 \\\\\n",
    "  12 & - 8 & 6 &  10 \\\\\n",
    "   3 & -13 & 9 &   3 \\\\\n",
    "  -6 &   4 & 1 & -18 \\\\\n",
    "\\end{bmatrix}\n",
    "\\end{align*}\n",
    "\n",
    "Comprobando computacionalmente la factorización $\\mathbf{LU=A}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.   0.   0.   0. ]\n",
      " [ 2.   1.   0.   0. ]\n",
      " [ 0.5  3.   1.   0. ]\n",
      " [-1.  -0.5  2.   1. ]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "#Ingreso de la matriz de coeficientes A, de forma manual\n",
    "\n",
    "L = np.array([[ 1,  0, 0, 0],\n",
    "              [ 2,  1, 0, 0],\n",
    "              [.5,  3, 1, 0],\n",
    "              [-1,-.5, 2, 1]])\n",
    "print(L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 6 -2  2  4]\n",
      " [ 0 -4  2  2]\n",
      " [ 0  0  2 -5]\n",
      " [ 0  0  0 -3]]\n"
     ]
    }
   ],
   "source": [
    "#Ingreso de la matriz de coeficientes A, de forma manual\n",
    "\n",
    "U = np.array([[ 6, -2, 2,  4],\n",
    "              [ 0, -4, 2,  2],\n",
    "              [ 0,  0, 2, -5],\n",
    "              [ 0,  0, 0, -3]])\n",
    "print(U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  6.  -2.   2.   4.]\n",
      " [ 12.  -8.   6.  10.]\n",
      " [  3. -13.   9.   3.]\n",
      " [ -6.   4.   1. -18.]]\n"
     ]
    }
   ],
   "source": [
    "print(L.dot(U))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'A' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-4-71d5d644394d>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mL\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mU\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m==\u001b[0m \u001b[0mA\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'A' is not defined"
     ]
    }
   ],
   "source": [
    "L.dot(U) == A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Implementación computacional del algoritmo LU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La implementación computacional se realizará modificando la función que se creó para el proceso de [eliminación](#Algoritmoeliminacion) para determinar la matriz $\\mathbf{L}$, pero esta teniendo en cuenta únicamente la matriz de coeficientes $\\mathbf{A}$ (no la ampliada), ya que la idea es poder emplear diferentes valores del vector de términos independientes $\\mathbf{b}$.\n",
    "\n",
    "Una forma de hacerlo es empleando el algoritmo de eliminación del método de eliminación de Gauss simple visto y almacenar los valores de los multiplicadores en las respectivas componentes de la nueva matriz $\\mathbf{L}$. Para ello, se usa una matriz identidad y se reemplazan los ceros por los valores de los multiplicadores en cada una de las posiciones indicadas. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def eliminacionLU(U):\n",
    "    \"\"\"\n",
    "    función en Python 3.8 que realiza el proceso de eliminación LU\n",
    "    en un Sistema de Ecuaciones Lineales (SEL) \n",
    "    \n",
    "    input: ingresa la matriz A[n x n]\n",
    "    \n",
    "    output: devuelve las matrices triangular inferior, L, y triangular superior,\n",
    "            U, de [n x n].\n",
    "    \"\"\"\n",
    "    n = U.shape[0]\n",
    "    L = np.identity(n)\n",
    "    \n",
    "    for k in range(0,n-1):\n",
    "        for i in range(k+1, n):\n",
    "            m = U[i][k]/U[k][k]\n",
    "            L[i,k] = m\n",
    "            for j in range(k, n):\n",
    "                U[i][j] = U[i][j] - m * U[k][j]\n",
    "    return (L, U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Esta aproximación es simple, pero computacionalmente costosa, ya que se crea una matriz con casi el 50% de los elementos iguales a cero. Una alternativa podría ser emplear la estructura de datos tipo [diccionario](https://docs.python.org/3/tutorial/datastructures.html#dictionaries) de Python, que facilita almacenar tanto el valor como las \"coordenadas\" matriciales (fila/columna) diferentes de cero y así no almacenar ceros innecesariamente. Se deja como actividad a ser realizada por el estudiante.\n",
    "\n",
    "Veamos ahora los resultados obtenidos tras ejecutar este algoritmo y la comprobación $\\mathbf{LU=A}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "U = A.copy()\n",
    "L,U = eliminacionLU(U)\n",
    "print(L)\n",
    "print(U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(L.dot(U))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L.dot(U) == A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Actividad para el estudiante:*** A qué considera que no todos los valores en esta comprobación son `True`?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplo computacional algoritmo LU: Cálculo de la Inversa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determinar la inversa de la matriz $\\mathbf{A}$ del ejemplo que hemos venido trabajando.\n",
    "\n",
    "La inversa de una matriz empleando la factorización $\\mathbf{LU}$ se puede obtener mediante el cálculo sucesivo del procedimiento descrito, cambiando cada vez el vector del lado derecho por un vector columna que tiene todos sus elementos iguales a cero,excepto por el valor de uno en la posición que corresponde a la columna de la matriz inversa a ser calculada. Este proceso generará cada vez como respuesta la columna correspondiente a la matriz inversa. Veamos.\n",
    "\n",
    "En el ejemplo anterior ya se obtuvo la factorización $\\mathbf{LU}$, ahora realicemos el procedimiento para cada vector $\\mathbf{b}$, y cada una de las respuestas $\\mathbf{x}$ serán las columnas de la matriz inversa.\n",
    "\n",
    "Primero, obtengamos la primera columna de la matriz inversa. Para ello, creamos un vector:\n",
    "\n",
    "\\begin{align}\n",
    "  b_1 &= \\begin{bmatrix}\n",
    "         1 \\\\\n",
    "         0 \\\\\n",
    "         0 \\\\\n",
    "         0\n",
    "       \\end{bmatrix}\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "b1 = np.array([1, 0, 0, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora vamos a resolver el sistema $\\mathbf{Ld=b_1}$ por sustitución hacia adelante. El algoritmo que tenemos de sustitución es hacia atrás. Toca modificarlo para que resuelva hacia adelante (progresiva)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sustitucionProg(A,b):\n",
    "    \"\"\"\n",
    "    función en Python 3.8 que realiza el proceso de sustitución progresiva\n",
    "    \n",
    "    input: Se ingresa la matriz triangular inferior L[n * n]\n",
    "    \n",
    "    output: devuelve el vector solución b[n * 1]\n",
    "    \"\"\"\n",
    "    \n",
    "    n = A.shape[0]\n",
    "    d = np.zeros(n)\n",
    "    d[0] = b[0]\n",
    "    for i in range(1,n):\n",
    "        suma = 0\n",
    "        for j in range(n):\n",
    "            suma = suma + A[i][j] * d[j]\n",
    "        d[i] = b[i] - suma\n",
    "    return d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  1.   -2.    5.5 -11. ]\n"
     ]
    }
   ],
   "source": [
    "d1 = sustitucionProg(L,b1)\n",
    "print(d1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Este resultado preliminar se emplea al lado derecho para resolver el sistema  $\\mathbf{Ux=d}$ mediante sustitución regresiva, donde el vector $\\mathbf{x}$ será la solución que corresponderá a la primera columna de la matriz inversa.\n",
    "\n",
    "Observe que el algoritmo de sustitución regresiva que se tienen implementado es con la matriz aumentada, por lo que toca antes aumentar la matriz $\\mathbf{U}$ con el vector de términos independientes, $\\mathbf{d}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  6.   -2.    2.    4.    1. ]\n",
      " [  0.   -4.    2.    2.   -2. ]\n",
      " [  0.    0.    2.   -5.    5.5]\n",
      " [  0.    0.    0.   -3.  -11. ]]\n",
      "(4, 5)\n"
     ]
    }
   ],
   "source": [
    "Ud1 = np.c_[U,d1]\n",
    "print(Ud1)\n",
    "print(np.shape(Ud1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x[0] = -3.4861\n",
      "x[1] = 8.2917\n",
      "x[2] = 11.9167\n",
      "x[3] = 3.6667\n",
      "[-3.48611111  8.29166667 11.91666667  3.66666667]\n"
     ]
    }
   ],
   "source": [
    "x1 = sustitucionReg(Ud1)\n",
    "print(x1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El vector x1 calculado corresponde a la primera columna de la matriz inversa. Repitiendo este procedimiento, cambiando cada vez el vector de términos independientes, para el cálculo de la segunda columna, sería:\n",
    "\n",
    "\\begin{align}\n",
    "  b_2 &= \\begin{bmatrix}\n",
    "         0 \\\\\n",
    "         1 \\\\\n",
    "         0 \\\\\n",
    "         0\n",
    "       \\end{bmatrix}\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x[0] = 2.1528\n",
      "x[1] = -4.7917\n",
      "x[2] = -6.9167\n",
      "x[3] = -2.1667\n",
      "[ 2.15277778 -4.79166667 -6.91666667 -2.16666667]\n"
     ]
    }
   ],
   "source": [
    "b2 = np.array([0, 1, 0, 0])\n",
    "d2 = sustitucionProg(L,b2)\n",
    "Ud2 = np.c_[U,d2]\n",
    "x2 = sustitucionReg(Ud2)\n",
    "print(x2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para el siguiente término se emplea el vector de términos independientes\n",
    "\n",
    "\\begin{align}\n",
    "  b_3 &= \\begin{bmatrix}\n",
    "         0 \\\\\n",
    "         0 \\\\\n",
    "         1 \\\\\n",
    "         0\n",
    "       \\end{bmatrix}\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x[0] = -0.6944\n",
      "x[1] = 1.4167\n",
      "x[2] = 2.1667\n",
      "x[3] = 0.6667\n",
      "[-0.69444444  1.41666667  2.16666667  0.66666667]\n"
     ]
    }
   ],
   "source": [
    "b3 = np.array([0, 0, 1, 0])\n",
    "d3 = sustitucionProg(L,b3)\n",
    "Ud3 = np.c_[U,d3]\n",
    "x3 = sustitucionReg(Ud3)\n",
    "print(x3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para el último término, se emplea el vector de términos independientes\n",
    "\n",
    "\\begin{align}\n",
    "  b_4 &= \\begin{bmatrix}\n",
    "         0 \\\\\n",
    "         0 \\\\\n",
    "         0 \\\\\n",
    "         1\n",
    "       \\end{bmatrix}\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x[0] = 0.3056\n",
      "x[1] = -0.5833\n",
      "x[2] = -0.8333\n",
      "x[3] = -0.3333\n",
      "[ 0.30555556 -0.58333333 -0.83333333 -0.33333333]\n"
     ]
    }
   ],
   "source": [
    "b4 = np.array([0, 0, 0, 1])\n",
    "d4 = sustitucionProg(L,b4)\n",
    "Ud4 = np.c_[U,d4]\n",
    "x4 = sustitucionReg(Ud4)\n",
    "print(x4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Por último, montando la matriz inversa con cada uno de los vectores obtenidos anteriormente. Debe tener presente que los vectores obtenidos son en forma de vector fila, por lo que hay que realizar la transposición de la misma."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-3.48611111  2.15277778 -0.69444444  0.30555556]\n",
      " [ 8.29166667 -4.79166667  1.41666667 -0.58333333]\n",
      " [11.91666667 -6.91666667  2.16666667 -0.83333333]\n",
      " [ 3.66666667 -2.16666667  0.66666667 -0.33333333]]\n"
     ]
    }
   ],
   "source": [
    "Ainv = np.matrix([x1, x2, x3, x4]).T\n",
    "print(Ainv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comprobando que efectivamente se obtuvo la matriz inversa, realizamos el cálculo: \n",
    "\n",
    "$$\\mathbf{A^{-1}A=I}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.00000000e+00 -1.77635684e-15  1.77635684e-15  3.55271368e-15]\n",
      " [ 7.10542736e-15  1.00000000e+00  3.55271368e-15  7.10542736e-15]\n",
      " [ 0.00000000e+00 -7.10542736e-15  1.00000000e+00  7.10542736e-15]\n",
      " [ 0.00000000e+00 -1.77635684e-15 -1.77635684e-15  1.00000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "Identity = Ainv.dot(A)\n",
    "print(Identity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Análisis de Error y condición del Sistema"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Condicionamiento del sistema"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La inversa de una matriz proporciona una forma para determinar si el sistema está mal condicionado. Se tienen tres métodos para esto:\n",
    "\n",
    "- Escalar la matriz de coeficientes.\n",
    "\n",
    "- Multplicar la inversa por la matriz de coeficientes original y estimar si el resultado es lo suficientemente cercano a la matriz identidad, es decir: $\\mathbf{A^{-1}A \\approx I}$\n",
    "\n",
    "- Invertir la matriz inversa y estimar si el resultado está lo suficientemente cerca a la matriz de coeficientes original, es decir, $\\mathbf{[A^{-1}]^{-1}\\approx A}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normas vectoriales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img07a_norms.gif?raw=true\" width=\"350\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://upload.wikimedia.org/wikipedia/commons/1/1f/Lp_space_animation.gif\">Wikipedia</a> </div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para medir errores y sensibilidad en la resolución de sistemas lineales, necesitamos alguna noción del \"*tamaño*\" de vectores y matrices. El concepto escalar de magnitud, módulo o valor absoluto se puede generalizar al concepto de normas para vectores y matrices.\n",
    "\n",
    "- ***[Norma](https://en.wikipedia.org/wiki/Norm_(mathematics)):*** *es una función que toma valores reales y que proporciona una medida del tamaño o \"longitud\" de entidades matemáticas multicomponentes, como los vectores y las matrices*.\n",
    "\n",
    "Aunque es posible una definición más general, todas las normas vectoriales que usaremos son instancias de $p-normas$, que para un entero $p>0$ y un vector $\\mathbf{x}$ de dimensión $n$ se definen por:\n",
    "\n",
    "$$\\|\\mathbf{x}\\|_p=\\left(\\sum_{i=1}^n |x_i|^p \\right)^{1/p}$$\n",
    "\n",
    "Algunos casos especiales de normas son:\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Norma-1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "También conocida como *Distancia Manhattan* o *Geometría del taxista* ([Taxicab geometry](https://en.wikipedia.org/wiki/Taxicab_geometry)). Es la distancia entre dos puntos como la suma de las diferencias (absolutas) de sus coordenadas. El nombre *Manhattan* alude al diseño en cuadrícula de la mayoría de las calles de la isla de *[Manhattan](http://www.mappery.com/map-of/Manhattan-Tourist-Map)*, lo que causa que el camino más corto que un auto puede tomar entre dos puntos de la ciudad tengan la misma distancia que dos puntos en la geometría del taxista. Fue presentado por [Hermann Minkowski](https://en.wikipedia.org/wiki/Hermann_Minkowski).\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img08_Manhattan.png?raw=true\" width=\"250\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://es.wikipedia.org/wiki/Geometr%C3%ADa_del_taxista#/media/Archivo:Manhattan_distance.svg\">Wikipedia</a> </div>\n",
    "\n",
    "\n",
    "***Definición:***\n",
    "\n",
    "$$\\|\\mathbf{x}\\|_1=\\sum_{i=1}^n |x_i|$$\n",
    "\n",
    "En un espacio de dimensión 2 (2D), un punto es representado por $(x,y)$. Considere dos puntos $P1: (x_1,y_1)$ y $P2: (x_2,y_2)$, la $Norma-1$, o distancia Manhattan, entre $P1$ y $P2$ está dada por:\n",
    "\n",
    "\n",
    "$$\\|\\mathbf{x}\\|_1=\\{|x_1-x_2|\\ +\\ |y_1-y_2|\\}$$\n",
    "\n",
    "\n",
    "***Propiedades:***\n",
    "\n",
    "Algunas propiedades de esta norma son:\n",
    "\n",
    "- Hay varios caminos (finitos) entre dos puntos cuya longitud es igual a la distancia de *Manhattan*\n",
    "\n",
    "\n",
    "- Un camino recto con una longitud igual a la distancia de *Manhattan* tiene dos movimientos permitidos:\n",
    "\n",
    "  - Vertical (una dirección)\n",
    "  - Horizontal (una dirección)\n",
    "\n",
    "\n",
    "- Para un punto dado, el otro punto a una distancia de *Manhattan* determinada se encuentra en un cuadrado:\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img09_L1Norm.png?raw=true\" width=\"150\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://miro.medium.com/max/219/1*cd3uQPINUGlpPaEganGG9Q.png\">Medium</a> </div>\n",
    "\n",
    "\n",
    "***Aplicaciones:***\n",
    "\n",
    "- ***Análisis de regresión:*** se utiliza en regresión lineal para encontrar una línea recta que se ajuste a un conjunto dado de puntos.\n",
    "\n",
    "\n",
    "- ***Detección comprimida:*** al resolver un sistema indeterminado de ecuaciones lineales, el término de regularización para el vector de parámetros se expresa en términos de la distancia de Manhattan. Este enfoque aparece en el marco de recuperación de señales llamado detección comprimida\n",
    "\n",
    "\n",
    "- ***Distribución de frecuencia:*** se utiliza para evaluar las diferencias en distribuciones de frecuencia discretas.\n",
    "\n",
    "\n",
    "- ***[Regularización](https://towardsdatascience.com/regularization-in-machine-learning-76441ddcf99a):*** en modelos de Machine Learning para evitar el overfitting."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Norma-2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "También llamada *[norma Euclideana](https://en.wikipedia.org/wiki/Euclidean_distance)* (o *distancia euclideana*). La *distancia euclidiana* es la distancia más corta entre dos puntos en un espacio de dimensión $N$, también conocido como *[espacio euclidiano](https://en.wikipedia.org/wiki/Euclidean_space)*. Se usa como una métrica común para medir la similitud entre dos puntos de datos y se usa en varios campos como geometría, minería de datos, aprendizaje profundo y otros.\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img10_Euclidean_distance_2d.svg?raw=true\" width=\"350\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://en.wikipedia.org/wiki/Euclidean_distance#/media/File:Euclidean_distance_2d.svg\">Wikipedia</a> </div>\n",
    "\n",
    "\n",
    "***Definición***\n",
    "\n",
    "$$\\|\\mathbf{x}\\|_2=\\sqrt{\\sum_{i=1}^n |x_i|^2}$$\n",
    "\n",
    "\n",
    "En un espacio de dimensión 2 (2D), un punto es representado por $(x,y)$. Considere dos puntos $p: (x_1,y_1)$ y $q: (x_2,y_2)$, la $Norma-2$, o *distancia euclideana*, entre $p$ y $q$ está dada por:\n",
    "\n",
    "\n",
    "$$\\|\\mathbf{x}\\|_2=\\sqrt{(x_1-x_2)^2 + (y_1-y_2)^2}$$\n",
    "\n",
    "\n",
    "***Propiedades:***\n",
    "\n",
    "- Hay un camino único entre dos puntos cuya longitud es igual a la distancia euclidiana.\n",
    "\n",
    "\n",
    "- Para un punto dado, el otro punto se encuentra en un círculo tal que la distancia euclidiana es fija. El radio del círculo es la distancia euclidiana fija.\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img11_L2Norm.png?raw=true\" width=\"150\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://miro.medium.com/max/219/1*y1BAWpMRMeHcO7O1XEbk-Q.png\">Medium</a> </div>\n",
    "\n",
    "\n",
    "\n",
    "***Aplicaciones:***\n",
    "\n",
    "- ***Geometría euclidiana:*** para encontrar la distancia más corta entre dos puntos en un espacio euclidiano y la longitud de una línea recta entre dos puntos.\n",
    "\n",
    "\n",
    "- ***Análisis de clústeres:*** esta métrica se usa comúnmente en algoritmos de clústeres como [K-means](https://en.wikipedia.org/wiki/K-means_clustering)\n",
    "\n",
    "\n",
    "- ***Ciencia de datos:*** se utiliza como una métrica simple para medir la [similitud](https://en.wikipedia.org/wiki/Similarity_measure) entre dos puntos de datos\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Norma-$\\infty$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "También conocida como *[distancia de Chebyshev](https://en.wikipedia.org/wiki/Chebyshev_distance)*. Es una métrica que es la distancia absoluta máxima entre dos pntos $N-dimensionales$. Tiene aplicaciones del mundo real en ajedrez, logística de almacén y muchos otros campos. También se conoce como distancia de Tchebychev, métrica máxima, distancia del tablero de ajedrez y métrica $L_{\\infty}$. Lleva el nombre de *[Pafnuty Chebyshev](https://en.wikipedia.org/wiki/Pafnuty_Chebyshev)*, matemático ruso que recibió conocido por su trabajo sobre probabilidad, estadística, mecánica, geometría analítica y teoría de números. \n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img12_ChebyshevDistance.png?raw=true\" width=\"250\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://en.wikipedia.org/wiki/Chebyshev_distance\">Wikipedia</a> </div>\n",
    "\n",
    "\n",
    "***Definición***\n",
    "\n",
    "\n",
    "$$\\|\\mathbf{x}\\|_{\\infty}=\\underset{1 \\leq i \\leq n}{max}|x_i|$$\n",
    "\n",
    "\n",
    "En un espacio de dimensión 2 (2D), un punto es representado por $(x,y)$. Considere dos puntos $P1: (x_1,y_1)$ y $P2: (x_2,y_2)$, la $Norma-\\infty$, entre $P1$ y $P2$ está dada por:\n",
    "\n",
    "\n",
    "$$\\|\\mathbf{x}\\|_{\\infty}=max\\{|x_1-x_2| ,|y_1-y_2|\\}$$\n",
    "\n",
    "\n",
    "\n",
    "***Propiedades***\n",
    "\n",
    "Las propiedades de la distancia de Chebyshev son:\n",
    "\n",
    "- La distancia bidimensional de Manhattan también tiene círculos en forma de cuadrados, con lados de longitud $\\sqrt{2}r$, orientados en un ángulo de $\\pi/4$ $(45°)$ con respecto a los ejes de coordenadas, por lo que la distancia plana de Chebyshev puede verse como equivalente por rotación y escalado a la distancia plana de Manhattan.\n",
    "\n",
    "- Una esfera formada usando la distancia de Chebyshev como métrica es un cubo con cada cara perpendicular a uno de los ejes de coordenadas, pero una esfera formada usando la distancia de Manhattan es un octaedro: estos son poliedros duales, pero entre cubos, solo el cuadrado (y un segmento de línea 1-dimensional) son politopos auto-duales.\n",
    "\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img13_L_infty_Norm.png?raw=true\" width=\"150\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://miro.medium.com/max/219/1*b7Zm7DSvo-u8D7-4rVs53Q.png\">Medium</a> </div>\n",
    "\n",
    "\n",
    "***Aplicaciones***\n",
    "\n",
    "Algunas aplicaciones de la distancia de Chebyshev son:\n",
    "\n",
    "- ***Ajedrez:*** el número mínimo de movimientos que necesita un rey para ir de un cuadrado a otro en un tablero de ajedrez es igual a la distancia de Chebyshev entre los centros de los cuadrados.\n",
    "\n",
    "\n",
    "- ***Logística de almacenaje:*** la distancia de Chebyshev se utiliza a veces en la logística de almacenaje (o bodegaje), ya que mide de forma eficaz el tiempo que tarda una grúa puente en mover un objeto.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Conclusión"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Todas estas normas dan los mismos resultados cualitativamente, pero en ciertas circunstancias puede ser más fácil trabajar con una norma particular analítica o computacionalmente. La $norma-1$ o la $norma-\\infty$ se usan generalmente para analizar la sensibilidad de las soluciones a los sistemas lineales. También usaremos la $norma-2$ más adelante en otros contextos. Las diferencias entre estas normas se ilustran en la siguiente figura, que muestra la esfera unitaria, $\\{\\mathbf{x}: \\|\\mathbf{x}\\| = 1\\}$, en dos dimensiones para cada norma.\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img07_norms.PNG?raw=true\" width=\"350\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://heath.cs.illinois.edu/scicomp/\">Heath, Michael. Scientific Computing. An introductory survey</a> </div>\n",
    "\n",
    "La norma de un vector es simplemente el factor por el cual la esfera unitaria correspondiente debe expandirse o encogerse para abarcar el vector. Por ejemplo, las normas tienen los siguientes valores para el vector:\n",
    "\n",
    "$$\\|\\mathbf{x}\\|_1=|(-1.6,1.2)|=2.8$$\n",
    "$$\\|\\mathbf{x}\\|_2=|(-1.6,1.2)|=2.0$$\n",
    "$$\\|\\mathbf{x}\\|_{\\infty}=|(-1.6,1.2)|=1.6$$\n",
    "\n",
    "en general, para cualquier vector $\\mathbf{x} \\in \\mathbb{R}^n$, tenemos que:\n",
    "\n",
    "$$\\|\\mathbf{x}\\|_1 \\ge \\|\\mathbf{x}\\|_2 \\ge \\|\\mathbf{x}\\|_{\\infty}$$\n",
    "\n",
    "Para cualquier norma vectorial, se cumplen las siguientes propiedades importantes, donde $\\mathbf{x}$ e $\\mathbf{y}$ son cualquier vector:\n",
    "\n",
    "- $\\|\\mathbf{x}\\|>0 \\text{ si } \\mathbf{x} \\ne 0$\n",
    "\n",
    "\n",
    "- $\\|\\gamma \\mathbf{x}\\|=|\\gamma| \\cdot \\|\\mathbf{x}\\| \\text{ para cualquier escalar } \\gamma$\n",
    "\n",
    "\n",
    "- $\\|\\mathbf{x+y}\\| \\le \\|\\mathbf{x}\\|+\\|\\mathbf{y}\\| \\text{ (desigualdad triangular)}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normas matriciales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "También se necesita alguna forma de medir el tamaño o la magnitud de las matrices. Todas las normas matriciales que se usarran en este capítulo se definen en términos de una norma vectorial subyacente. Específicamente, dada una norma vectorial, definimos la [norma matricial](https://en.wikipedia.org/wiki/Matrix_norm) correspondiente de una matriz $\\mathbf{A}$ de la siguiente manera:\n",
    "\n",
    "\n",
    "$$\\|\\mathbf{A}\\|=\\underset{x \\ne 0}{max} \\frac{\\|\\mathbf{Ax}\\|}{\\|\\mathbf{x}\\|}$$\n",
    "\n",
    "Se dice que dicha norma matricial está subordinada a la norma vectorial. Intuitivamente, la norma de una matriz mide el estiramiento máximo que la matriz hace a cualquier vector, medido en la norma de vector dada.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Norma-1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se realiza una sumatoria de los valores absolutos de los coeficientes para cada columna, y la mayor de estas sumatorias se toma como la norma. Esto se conoce como la norma columna-suma.\n",
    "  \n",
    "$$\\|\\mathbf{A}\\|_1=\\underset{1 \\leq j \\leq n}{max}\\sum_{i=1}^n |a_{ij}|$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Norma-2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una forma útil de recordar es que las normas de la matriz concuerdan con las normas vectoriales correspondientes para una matriz $n \\times 1$. la $norma-2$ también es llamada *norma de Frobenius* o de *Hilbert-Shmidt*. Desafortunadamente, la norma matricial correspondiente a la $norma-2$ vectorial no es tan simple de calcular; resulta ser igual a la raíz cuadrada del valor propio más grande de la matriz $\\mathbf{A}$ ([traza](https://en.wikipedia.org/wiki/Trace_(linear_algebra))) o como el [valor singular](https://en.wikipedia.org/wiki/Singular_value) más grande de $\\mathbf{A}$.\n",
    "\n",
    "\n",
    "$$\\|\\mathbf{A}\\|_2=\\sqrt{\\sum_{i=1}^n \\sum_{j=1}^n |a_{ij}|^2} = \\sqrt{traza(\\mathbf{A}^T\\mathbf{A})} = \\sqrt{\\sum_{i=1}^n \\sigma_i^2(\\mathbf{A})}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Norma-$\\infty$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar a la $norma-1$ para matrices, se puede hacer para los renglones y resulta una matriz-uniforme o norma renglón-suma\n",
    "\n",
    "$$\\|\\mathbf{A}\\|_{\\infty}=\\underset{1 \\leq i \\leq n}{max}\\sum_{j=1}^n |a_{ij}|$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Conclusión"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Las normas matriciales que se han definido satisfacen las siguientes propiedades importantes, donde $\\mathbf{A}$ y $\\mathbf{B}$ son matrices cualesquiera:\n",
    "\n",
    "- $\\|\\mathbf{A}\\|>0 \\text{ si } \\mathbf{A} \\ne 0$\n",
    "\n",
    "\n",
    "- $\\|\\gamma \\mathbf{A}\\|=|\\gamma| \\cdot \\|\\mathbf{A}\\| \\text{ para cualquier escalar } \\gamma$\n",
    "\n",
    "\n",
    "- $\\|\\mathbf{A+B}\\| \\le \\|\\mathbf{A}\\|+\\|\\mathbf{B}\\|$\n",
    "\n",
    "\n",
    "- $\\|\\mathbf{AB}\\| \\leq \\|\\mathbf{A}\\| \\cdot \\|\\mathbf{B}\\|$\n",
    "\n",
    "\n",
    "- $\\|\\mathbf{Ax}\\| \\leq \\|\\mathbf{A}\\| \\cdot \\|\\mathbf{x}\\| \\text{ para cualquier vector } \\mathbf{x}$\n",
    "\n",
    "En un tratamiento más general, las tres primeras propiedades pueden tomarse como la definición de una norma matricial. Las dos propiedades restantes, conocidas como condiciones submultiplicativas o de consistencia, pueden o no ser válidas para estas normas matriciales más generales, pero siempre son válidas para las normas matriciales subordinadas a las $p-normas$ vectoriales."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Número de condición de una matriz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Definición"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El <a id='NCOND'></a>[número de condición](https://en.wikipedia.org/wiki/Condition_number) de una matriz cuadrada no singular $\\mathbf{A}$ con respecto a una norma dada se define como:\n",
    "\n",
    "$$cond(\\mathbf{A})=\\|\\mathbf{A}\\| \\cdot \\|\\mathbf{A^{-1}}\\|$$\n",
    "\n",
    "por convención, $cond(\\mathbf{A})=\\infty$, si $\\mathbf{A}$ es singular. Ya que\n",
    "\n",
    "$$\\|\\mathbf{A}\\| \\cdot \\|\\mathbf{A^{-1}}\\|=\\left(\\underset{x \\ne 0}{max} \\frac{\\|\\mathbf{Ax}\\|}{\\|\\mathbf{x}\\|}\\right)\\cdot \\left(\\underset{x \\ne 0}{min} \\frac{\\|\\mathbf{Ax}\\|}{\\|\\mathbf{x}\\|}\\right)^{-1}$$\n",
    "\n",
    "El número de condición de una matriz mide la relación entre el estiramiento máximo que hace la matriz a cualquier vector distinto de cero y el encogimiento máximo. Este concepto es consistente con la noción general de número de condición definido en el capítulo 1 - Teoría de errores. el número de condición de la matriz limita la razón del cambio relativo en la solución de un sistema lineal a un cambio relativo en los datos de entrada.\n",
    "\n",
    "El número de condición es una medida de qué tan cerca está una matriz de ser singular: una matriz con un número de condición grande es casi singular, mientras que una matriz con un número de condición cercano a $1$ está lejos de ser singular. Tenga en cuenta que el determinante de una matriz no es un buen indicador de la singularidad cercana: aunque una matriz $\\mathbf{A}$ es singular si $det(\\mathbf{A})=0$, la magnitud de un determinante distinto de cero, grande o pequeño, no da información sobre qué tan cerca del singular es la matriz puede ser. Por ejemplo, $det(\\alpha\\mathbf{I}n)= \\alpha^n$, que puede ser arbitrariamente pequeño para $|\\alpha|<1$, pero la matriz está perfectamente acondicionada para cualquier valor de $\\alpha$ distinto de cero, con un número de condición cercano a $1$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Propiedades"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Algunas propiedades importantes del número de condición son:\n",
    "\n",
    "\n",
    "- Para cualquier matriz $\\mathbf{A}$, $cond(\\mathbf{A}) \\ge 1$.\n",
    "\n",
    "\n",
    "- Para la matriz identidad, $cond(\\mathbf{I}) = 1$.\n",
    "\n",
    "\n",
    "- Para cualquier matriz de permutacion $\\mathbf{P}$, $cond(\\mathbf{P}) = 1$.\n",
    "\n",
    "\n",
    "- Para cualquier matriz $\\mathbf{A}$ y un escalar $\\gamma$ diferente de cero, $cond(\\gamma \\mathbf{A}) = cond(\\mathbf{A})$.\n",
    "\n",
    "\n",
    "- Para cualquier matriz diagonal $\\mathbf{D}=diag(d_i)$, $cond(\\mathbf{D}) = (max|d_i|) / (min|d_i|)$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cálculo del número de condición"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La utilidad del número de condición radica en evaluar la precisión de las soluciones de los sistemas lineales. Dado que la definición del número de condición implica el cálculo de la inversa de la matriz, que es obviamente una tarea no trivial. De hecho, calcular el número de condición literalmente requeriría mucho más trabajo que resolver el sistema lineal cuya precisión debe evaluarse utilizando el número de condición. En la práctica, por lo tanto, el número de condición se estima simplemente, quizás dentro de un orden de magnitud, como un subproducto relativamente económico del proceso de solución.\n",
    "\n",
    "La norma matricial $\\|\\mathbf{A}\\|$ se calcula fácilmente como la máxima suma absoluta de columnas (o suma de filas, según la norma utilizada). Es la estimación de $\\|\\mathbf{A}^{-1}\\|$ a bajo costo lo que presenta un desafío. De las propiedades de las normas, sabemos que si $\\mathbf{z}$ es la solución de $\\mathbf{Az=y}$, entonces\n",
    "\n",
    "$$\\frac{\\|\\mathbf{z}\\|}{\\|\\mathbf{y}\\|} \\leq \\|\\mathbf{A}^{-1}\\|$$\n",
    "\n",
    "y el límite se logra para algún vector $\\mathbf{y}$ elegido de manera \"óptima\". Por lo tanto, deseamos elegir un vector $\\mathbf{y}$ de modo que la relación $\\|\\mathbf{z}\\|/\\|\\mathbf{y}\\|$ sea lo más grande posible y, por lo tanto, sea una estimación razonable de $\\|\\mathbf{A}^{-1}\\|$. Encontrar la $\\mathbf{y}$ óptima sería prohibitivamente caro, pero se puede obtener una aproximación útil mucho más barata. Una heurística es elegir $\\mathbf{y}$ como la solución del sistema $\\mathbf{A}^T\\mathbf{c}=\\mathbf{c}$, donde $\\mathbf{c}$ es un vector cuyas componentes son $1$, con los signos elegidos sucesivamente para hacer que la $\\mathbf{y}$ resultante sea lo más grande posible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Precisión de las soluciones"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Residual de una solución"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Intuitivamente, la forma más obvia de verificar la validez de una solución es sustituirla en la ecuación para ver qué tan cerca coinciden los dos lados. El vector residual de una solución aproximada $\\mathbf{\\hat{x}}$ al sistema lineal $\\mathbf{Ax=b}$ de $n\\times n$ se define como \n",
    "\n",
    "$$\\mathbf{r=b-A\\hat{x}}$$\n",
    "\n",
    "En teoría, si $\\mathbf{A}$ no es singular, entonces el error $\\|\\mathbf{\\hat{x}}-\\mathbf{b}\\|=\\mathbf{0}$ si y solo si $\\|\\mathbf{r}\\|=\\mathbf{0}$. En la práctica, sin embargo, estas cantidades no son necesariamente pequeñas simultáneamente. Si la solución calculada $\\mathbf{\\hat{x}}$ satisface exactamente\n",
    "\n",
    "$$(\\mathbf{A+E})\\mathbf{\\hat{x}}=\\mathbf{b}$$\n",
    "\n",
    "entonces\n",
    "\n",
    "$$\\|\\mathbf{r}\\|=\\|\\mathbf{b-A\\hat{x}}\\|=\\|\\mathbf{E\\hat{x}}\\| \\le \\|\\mathbf{E}\\|\\cdot\\|\\mathbf{\\hat{x}}\\|$$\n",
    "\n",
    "llegando a la siguiente desigualdad\n",
    "\n",
    "$$\\frac{\\|\\mathbf{r}\\|}{\\|\\mathbf{A}\\|\\cdot\\|\\mathbf{\\hat{x}}\\|} \\le \\frac{\\|\\mathbf{E}\\|}{\\|\\mathbf{A}\\|}$$\n",
    "\n",
    "relacionando el residuo relativo con el cambio relativo en la matriz. Por lo tanto, un residuo relativo grande implica un gran error hacia atrás en la matriz, lo que significa que el algoritmo utilizado para calcular la solución es inestable.\n",
    "\n",
    "pero cuán grande debería ser $\\|\\mathbf{E}\\|$ en la práctica? Se puede demostrar que para la factorización $\\mathbf{LU}$ en la *eliminación de Gauss simple*, un valor límite es de la forma:\n",
    "\n",
    "$$\\frac{\\|\\mathbf{E}\\|}{\\|\\mathbf{A}\\|} \\le \\rho n \\epsilon_{mach}$$\n",
    "\n",
    "donde, $\\rho$ es llamado *factor de crecimiento*, que es la relación entre la entrada más grande de $\\mathbf{U}$ y la entrada más grande de $\\mathbf{A}$. Sin pivotar, $\\rho$ puede ser arbitrariamente grande y, por lo tanto, la eliminación gaussiana sin pivotar es inestable, como ya se ha visto. Con pivotamiento parcial, el factor de crecimiento aún puede ser tan grande como $2^{n-1}$ (ya que en el peor de los casos el tamaño de las entradas puede duplicarse en cada etapa de eliminación), pero tal comportamiento es extremadamente raro. En la práctica, hay poco o ningún crecimiento en el tamaño de las entradas, por lo que\n",
    "\n",
    "$$\\frac{\\|\\mathbf{E}\\|}{\\|\\mathbf{A}\\|} \\approx n \\epsilon_{mach}$$\n",
    "\n",
    "Esta relación significa que resolver un sistema lineal por eliminación gaussiana con pivoteo parcial seguido de sustitución hacia atrás casi siempre produce un residuo relativo muy pequeño, independientemente de cuán mal condicionado pueda estar el sistema. Por lo tanto, un residuo relativo pequeño no es necesariamente un buen indicador de que una solución calculada está cerca de la solución \"verdadera\" a menos que el sistema esté bien acondicionado.\n",
    "\n",
    "El pivote completo produce un factor de crecimiento aún menor, tanto en la teoría como en la práctica, pero el margen adicional de estabilidad que proporciona generalmente no vale la pena el gasto adicional."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplo numérico"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usando aritmética de tres decimales, resolver el sistema\n",
    "\n",
    "\\begin{align*}\n",
    "\\left[\\begin{array}{cc}\n",
    "  0.641 & 0.242 \\\\\n",
    "  0.321 & 0.121 \\\\\n",
    " \\end{array}\\right]\n",
    "\\begin{Bmatrix}\n",
    "  x_{1}  \\\\\n",
    "  x_{2}  \\\\\n",
    "\\end{Bmatrix}\n",
    "= \\begin{Bmatrix}\n",
    "  0.883  \\\\\n",
    "  0.442  \\\\\n",
    "\\end{Bmatrix}\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "El proceso de eliminación con pivoteo parcial conduce a un sistema triangular de la forma\n",
    "\n",
    "\\begin{align*}\n",
    "\\left[\\begin{array}{cc}\n",
    "  0.641 & 0.242 \\\\\n",
    "  0 & -0.000242 \\\\\n",
    " \\end{array}\\right]\n",
    "\\begin{Bmatrix}\n",
    "  x_{1}  \\\\\n",
    "  x_{2}  \\\\\n",
    "\\end{Bmatrix}\n",
    "= \\begin{Bmatrix}\n",
    "  0.883  \\\\\n",
    "  -0.000383  \\\\\n",
    "\\end{Bmatrix}\n",
    "\\end{align*}\n",
    "\n",
    "y el proceso de sustitución regresiva da la respuesta\n",
    "\n",
    "\\begin{align*}\n",
    "x = \\begin{Bmatrix}\n",
    "  0.782  \\\\\n",
    "  1.58  \\\\\n",
    "\\end{Bmatrix}\n",
    "\\end{align*}\n",
    "\n",
    "El residuo para esta solución es\n",
    "\n",
    "\\begin{align*}\n",
    "r=b-Ax = \\begin{Bmatrix}\n",
    "  -0.000622  \\\\\n",
    "  -0.000202  \\\\\n",
    "\\end{Bmatrix}\n",
    "\\end{align*}\n",
    "\n",
    "\n",
    "que es lo suficientemente pequeño, como era lo esperado, usando solo aritmética de tres dígitos. Sin embargo, se ve fácilmente que la solución exacta para este sistema es\n",
    "\n",
    "\\begin{align*}\n",
    "x = \\begin{Bmatrix}\n",
    "  1.00  \\\\\n",
    "  1.00  \\\\\n",
    "\\end{Bmatrix}\n",
    "\\end{align*}\n",
    "\n",
    "de modo que el error es casi tan grande como la solución. La causa de este fenómeno es que la matriz es casi singular (su número de condición es > 4000). La división que determina $x_2$ es entre dos cantidades que están en el orden del error de redondeo y, por lo tanto, el resultado es esencialmente arbitrario. Sin embargo, por diseño, cuando este valor arbitrario para $x_2$ se sustituye en la primera ecuación, se calcula un valor para $x_1$ de modo que se satisfaga la primera ecuación. Por lo tanto, obtenemos un pequeño residuo, pero una mala solución."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Estimación de la precisión"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Además de ser un indicador confiable de casi singularidad, el número de condición también proporciona una estimación cuantitativa del error en la solución calculada para un sistema lineal. Sea $\\mathbf{x}$ la solución del sistema lineal no singular $\\mathbf{Ax=b}$, y sea $\\mathbf{\\hat{x}}$ la solución del sistema $\\mathbf{A\\hat{x}}=\\mathbf{b}+\\Delta\\mathbf{b}$ con un lado derecho alterado. Si definimos $\\Delta \\mathbf{x}=\\mathbf{\\hat{x}}-\\mathbf{x}$, entonces tenemos\n",
    "\n",
    "$$\\mathbf{b}+\\Delta\\mathbf{b}=\\mathbf{A\\hat{x}}=\\mathbf{A}(\\mathbf{x}+\\Delta\\mathbf{x})=\\mathbf{Ax}+\\mathbf{A}\\Delta\\mathbf{x}$$\n",
    "\n",
    "Como $\\mathbf{Ax=b}$, debemos tener $\\mathbf{A}\\Delta\\mathbf{x}=\\Delta\\mathbf{b}$, y por lo tanto $\\Delta\\mathbf{x}=\\mathbf{A}^{-1}\\Delta\\mathbf{b}$, entonces\n",
    "\n",
    "$$\\mathbf{b}=\\mathbf{Ax} \\Rightarrow \\|\\mathbf{b}\\| \\le \\|\\mathbf{A}\\| \\cdot \\|\\mathbf{x}\\|$$\n",
    "\n",
    "y\n",
    "\n",
    "$$\\Delta \\mathbf{x}=\\mathbf{A}^{-1}\\Delta\\mathbf{b} \\Rightarrow \\|\\Delta \\mathbf{x}\\| \\le \\|\\mathbf{A}^{-1}\\| \\cdot \\|\\Delta \\mathbf{b}\\|$$\n",
    "\n",
    "usando la definición del número de condición de una matriz $\\mathbf{A}$, se llega a la estimación\n",
    "\n",
    "$$\\frac{\\|\\Delta \\mathbf{x}\\|}{\\|\\mathbf{x}\\|} \\le cond(\\mathbf{A})\\frac{\\|\\Delta\\mathbf{b}\\|}{\\|\\mathbf{b}\\|}$$\n",
    "\n",
    "Por lo tanto, el número de condición de la matriz determina el posible cambio relativo en la solución debido a un cambio relativo dado en el vector del lado derecho, independientemente del algoritmo utilizado para calcular la solución. Un resultado similar es válido para cambios relativos en las entradas de la matriz $\\mathbf{A}$. Si $\\mathbf{Ax=b}$ y\n",
    "\n",
    "$$(\\mathbf{A}+\\mathbf{E})\\mathbf{\\hat{x}}=\\mathbf{b}$$\n",
    "\n",
    "entonces\n",
    "\n",
    "$$\\mathbf{x}-\\mathbf{\\hat{x}}=\\mathbf{A}^{-1}(\\mathbf{b}-\\mathbf{A}\\mathbf{\\hat{x}})=\\mathbf{A}^{-1}\\mathbf{E}\\mathbf{\\hat{x}}$$\n",
    "\n",
    "así que\n",
    "\n",
    "$$\\|\\Delta \\mathbf{x}\\| \\le \\|\\mathbf{A}^{-1}\\| \\cdot \\|\\mathbf{E}\\| \\cdot \\|\\mathbf{\\hat{x}}\\|$$\n",
    "\n",
    "que conduce al estimado\n",
    "\n",
    "$$\\frac{\\|\\Delta \\mathbf{x}\\|}{\\|\\mathbf{\\hat{x}}\\|} \\le cond(\\mathbf{A})\\frac{\\|\\mathbf{E}\\|}{\\|\\mathbf{A}\\|}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Interpretación geométrica"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una interpretación geométrica en dos dimensiones de estos resultados de sensibilidad es que si las dos líneas rectas definidas por las dos ecuaciones son casi paralelas, entonces su punto de intersección no se define claramente si las líneas son un poco borrosas debido a errores de redondeo u otros fuentes de error. Si, por otro lado, las líneas están lejos de ser paralelas, digamos casi perpendiculares, entonces su intersección está relativamente definida. Estos dos casos se ilustran en la siguiente figura, donde las líneas punteadas indican la región de incertidumbre para cada línea continua, de modo que el punto de intersección en cada caso podría estar en cualquier lugar dentro del paralelogramo sombreado. Por lo tanto, un gran número de condiciones se asocia con una gran incertidumbre en la solución.\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img14_illconditioned.PNG?raw=true\" width=\"500\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://heath.cs.illinois.edu/scicomp/\">Heath, Michael. Scientific Computing. An introductory survey</a> </div>\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Conclusión"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "En resumen, si los datos de entrada son exactos a la precisión de la máquina, entonces una estimación razonable del error relativo en la solución calculada para un sistema lineal viene dada por\n",
    "\n",
    "$$\\frac{\\|\\mathbf{\\hat{x}}-\\mathbf{x}\\|}{\\|\\mathbf{x}\\|}\\approx cond(\\mathbf{A})\\epsilon_{mach}$$\n",
    "\n",
    "Una forma sencilla de interpretar estos resultados es que la solución calculada pierde aproximadamente $log_{10} (cond(\\mathbf{A}))$ dígitos decimales de precisión en relación con la precisión de la entrada. En el ejemplo visto, con un número de condición mayor que $10^3$, perdimos toda la precisión de tres dígitos disponible y obtuvimos una solución arbitraria.\n",
    "\n",
    "El análisis anterior estima el error relativo en los componentes más grandes del vector solución. El error relativo en los componentes más pequeños puede ser mucho mayor, porque una norma vectorial está dominada por los componentes más grandes de un vector. Se pueden obtener límites de error por componentes, pero son algo más complicados de calcular. Los límites por componentes son de particular interés cuando el sistema está mal escalado.\n",
    "\n",
    "El número de condición de una matriz se ve afectado por la escala de la matriz. Un número de condición elevado puede deberse simplemente a un escalado deficiente, así como a una singularidad cercana. Cambiar la escala de la matriz puede ayudar a lo primero, pero no a lo segundo.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ejemplo numérico - Matriz de Hilbert "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se desea evaluar el condicionamiento de una matriz de [$Hilbert$](https://en.wikipedia.org/wiki/Hilbert_matrix). Emplearemos los métodos que trae el lenguaje de programación `python` para realizar la mayoría de los cálculos. \n",
    "\n",
    "Lo primero es cargar las bibliotecas numéricas y científicas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import hilbert, norm\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Seguidamente montamos una matriz de $Hilbert$ de dimensión $n$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "n = 3\n",
    "H = hilbert(n)\n",
    "print(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normalizando la matriz $H$, de tal forma que el mayor elemento en cada fila sea de valor unitario"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = np.array([H[i,:]/H[i,0] for i in range(n)])\n",
    "print(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Empleamos una de la normas vistas, en este caso usaremos la norma $L_{\\infty}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Hinf = norm(H, np.inf)\n",
    "print(Hinf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ahora determinamos la inversa de dicha matriz y calculamos su respectiva norma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Hinv = np.linalg.inv(H)\n",
    "print(Hinv)\n",
    "Hinv_inf = norm(Hinv, np.inf)\n",
    "print(Hinv_inf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Por último, calculamos el [número de condición](#NCOND):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ncondH = Hinf * Hinv_inf\n",
    "print(ncondH)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Qué análisis hace acerca de este resultado?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tipos especiales de sistemas lineales"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hasta el momento solo hemos tratado un tipo de matriz llamada *llena* (*dense matrix*), es decir, que la mayoría de sus elementos son diferentes de cero ($a_{ij}\\neq 0; i, j = 1, 2, 3, \\ldots, n$). Un ahorro considerable en recurso computacional (tiempo, almacenamiento, etc) se puede obtener si se tiene en cuenta el tipo de matriz de coeficientes resultante. Algunos ejemplos de tipos especiales son:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Simétrica"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una matriz se dice que es [simétrica](https://en.wikipedia.org/wiki/Symmetric_matrix), respecto a su diagonal principal, si los elementos $a_{i,j}=a_{j,i} \\text{ para todo } i,j = 1,2,3,\\ldots, n$, entonces se cumple que $\\mathbf{A}=\\mathbf{A}^T$, es decir, la matriz y su transpuesta son la misma.\n",
    "\n",
    "\\begin{bmatrix}\n",
    "a_{11} & a_{12} & a_{13} & \\ldots & a_{1n}\\\\\n",
    "a_{12} & a_{22} & a_{23} & \\ldots & a_{2n}\\\\\n",
    "a_{13} & a_{23} & a_{33} & \\ldots & a_{3n}\\\\\n",
    "\\vdots & \\vdots & \\vdots & \\ddots & \\vdots\\\\\n",
    "a_{1n} & a_{2n} & a_{3n} & \\ldots & a_{nn}\\\\\n",
    "\\end{bmatrix}\n",
    "\n",
    "En este caso, solo sería necesario almacenar una de las dos triangulares (superior o inferior), reduciendo casi a la mitad el tamaño del almacenamiento.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Dispersa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una matriz dispersa ([sparse](https://en.wikipedia.org/wiki/Sparse_matrix) en inglés), es una matriz en la que la mayoría de sus elementos son cero. De forma contraria, una matriz en la que la mayoría de sus elementos son diferentes de cero se denomina *densa*.\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img15_Sparse.png?raw=true\" width=\"250\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a href=\"https://upload.wikimedia.org/wikipedia/commons/8/8a/Finite_element_sparse_matrix.png\">Wikipedia</a> </div>\n",
    "\n",
    "\n",
    "En este tipo de matrices es importante emplear algoritmos y estructuras de datos que permitan tomar ventaja de la estructura dispersa de la matriz, esto para optimizar almacenamiento en memoria y manipulación al realizar operaciones sin necesidad de considerar la gran cantidad de elementos cero, que no aportan información al problema. \n",
    "\n",
    "La cantidad de memoria requerida para almacenar una matriz es proporcional a su tamaño ($n\\times\\ n$, en el caso de matrices cuadradas). Es por esto que, en el caso de una matriz dispersa, se pueden realizar reducciones sustanciales de los requisitos de memoria almacenando solo las entradas distintas de cero. Dependiendo del número y la distribución de las entradas distintas de cero, se pueden usar diferentes estructuras de datos y producir grandes ahorros de memoria en comparación con el enfoque básico. La compensación es que acceder a los elementos individuales se vuelve más complejo y se necesitan estructuras adicionales para poder recuperar la matriz original sin ambigüedades.\n",
    "\n",
    "Los formatos se pueden dividir en dos grupos:\n",
    "\n",
    "\n",
    "- Aquellos que soportan modificaciones eficientes, como DOK (Diccionario de claves), LIL (Lista de listas) o COO (Lista de coordenadas). Estos se utilizan normalmente para construir las matrices.\n",
    "\n",
    "\n",
    "- Aquellos que admiten operaciones matriciales y de acceso eficientes, como CSR (Compressed Sparse Row) o CSC (Compressed Sparse Column).\n",
    "\n",
    "A continuación se presentan algunas referencias bibliográficas donde podrá profundizar en este tipo de estructuras de almacenamiento:\n",
    "\n",
    "\n",
    "- En la página de la [netlib.org](https://www.netlib.org/) encuentra el libro en formato virtual de *[Jack Dongarra](http://www.netlib.org/utk/people/JackDongarra/)* - [Templates for Solution of Algebraic Eigenvalues Problems: A practical Guide](http://www.netlib.org/utk/people/JackDongarra/etemplates/book.html) y un capítulo dedicado a [algoritmos de almacenamiento para matrices dispersas](http://www.netlib.org/utk/people/JackDongarra/etemplates/node372.html).\n",
    "\n",
    "\n",
    "- También podrá encontrar información relevante acerca de estos métodos en el libro [Numerical Recipes - The Art of Scientific Computing](http://numerical.recipes/), editado por la Universidad de Cambridge.\n",
    "\n",
    "\n",
    "- La biblioteca matemática [Scipy.org](https://docs.scipy.org/doc/scipy/reference/sparse.html) cuenta con un paquete de funciones para tratar este tipo de matrices.\n",
    "\n",
    "\n",
    "- En la página [machinelearningmastery.com](https://machinelearningmastery.com/sparse-matrices-for-machine-learning/) podrá encontrar una breve introducción a las matrices dispersas aplicada a ML.\n",
    "\n",
    "\n",
    "- En la página [SuiteSparse Matrix Collection](https://sparse.tamu.edu/) encontrará un conjunto de matrices dispersas recopilados de una amplia gama de aplicaciones y que podrá descargar y usar para realizar sus propias evaluaciones."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Definida positiva"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una matriz [definida positiva](https://en.wikipedia.org/wiki/Definite_symmetric_matrix) es una matriz para la que se cumple que el producto $\\mathbf{x}^T\\mathbf{Ax}>0$ para todo vector $\\mathbf{x} \\neq \\mathbf{0}$. Este tipo de matrices permite tener una factorización $A=\\mathbf{LU}=\\mathbf{LL}^T$, también es llamada como *[Factorización de cholesky](https://en.wikipedia.org/wiki/Cholesky_decomposition)*. El proceso de eliminación para la $k$-ésima fila es:\n",
    "\n",
    "- elementos diferentes a la diagonal principal:\n",
    "\n",
    "$$l_{ki}=\\frac{a_{ki}-\\sum \\limits_{j=1}^{i-1}l_{ij}l_{kj}}{l_{ii}}, \\text{ para } i = 1,2, \\ldots, k-1$$\n",
    "\n",
    "\n",
    "- elemento de la diagonal principal:\n",
    "$$l_{kk}=\\sqrt{a_{kk}-\\sum \\limits_{j=1}^{k-1}l_{kj}^2}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Ejemplo: Factorización de Cholesky"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dada la siguiente matriz\n",
    "\n",
    "\\begin{bmatrix}\n",
    " 6 &  15 &  55 \\\\\n",
    "15 &  55 & 225 \\\\\n",
    "55 & 225 & 979 \\\\\n",
    "\\end{bmatrix}\n",
    "\n",
    "Encontrar su descomposición $\\mathbf{LL}^T$ (factorización de Cholesky)\n",
    "\n",
    "- para $k=1$: (primera fila)\n",
    "\n",
    "$$l_{11}=\\sqrt{a_{11}}=\\sqrt{6}=2.4495$$\n",
    "\n",
    "\n",
    "- para $k=2$\n",
    "$$l_{21}=\\frac{a_{21}}{l_{11}}=\\frac{15}{2.4495}=6.1237$$\n",
    "\n",
    "$$l_{22}=\\sqrt{a_{22}-l_{21}^2}=\\sqrt{55-(6.1237)^2}=4.1833$$\n",
    "\n",
    "\n",
    "- para $k=3$\n",
    "\n",
    "$$l_{31}=\\frac{a_{31}}{l_{11}}=\\frac{55}{2.4495}=22.454$$\n",
    "\n",
    "$$l_{32}=\\frac{a_{32}-l_{21}l{31}}{l_{22}}=\\frac{225-6.1237(22.454)}{4.1833}=20.916$$\n",
    "\n",
    "$$l_{33}=\\sqrt{a_{33}-l_{31}^2-l_{32}^2}=\\sqrt{979-(22.454)^2-(20.916)^2}=6.1106$$\n",
    "\n",
    "con esto, la descomposición de Chlesky queda:\n",
    "\n",
    "\\begin{align*}L=\n",
    "\\left[\\begin{array}{cccc}\n",
    " 2.4495 &   &   \\\\\n",
    "6.1237 &  4.1833 &  \\\\\n",
    "22.454 & 20.916 & 6.1106 \\\\\n",
    "  \\end{array}\\right]\n",
    "\\end{align*}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Implementación computacional algoritmo de cholesky"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Chapra\n",
    "\n",
    "def cholesky(A):\n",
    "    n = A.shape[0]\n",
    "    for k in range(n):\n",
    "        print(\"k:\", k)\n",
    "        for i in range(k-1):\n",
    "            print(\"i: \", i)\n",
    "            suma = 0\n",
    "            print (\"suma: \", suma)\n",
    "            for j in range(i-1):\n",
    "                print(\"k:\", k)\n",
    "                print(\"i: \", i)\n",
    "                print(\"j: \", j)\n",
    "                suma += A[i,j] * A[k,j]\n",
    "            A[k,i] = (A[k,i] - suma) / A[i,i]\n",
    "        suma = 0\n",
    "        for j in range(k-1):\n",
    "            suma += A[k,j]**2\n",
    "        A[k,k] = np.sqrt(A[k,k] - suma)\n",
    "    return A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "A = np.array([[ 6, 15, 55],\n",
    "              [15, 55, 225],\n",
    "              [55, 225, 979]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = cholesky(A)\n",
    "print(L)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "    for k in range(0,n-1):\n",
    "        for i in range(k+1, n):\n",
    "            m = U[i][k]/U[k][k]\n",
    "            L[i,k] = m\n",
    "            for j in range(k, n):\n",
    "                U[i][j] = U[i][j] - m * U[k][j]\n",
    "    return (L, U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Heath\n",
    "\n",
    "def choleskyH(A):\n",
    "    n = A.shape[0]\n",
    "    print(\"n: \", n)\n",
    "    for j in range(0,n):\n",
    "        print(\"columna j: \", j)\n",
    "        for k in range(0,j-1):\n",
    "            print(\"columnas anteriores a j: \", k)\n",
    "            for i in range(j, n):\n",
    "                print(\"resta multiplo col k de col j: \", i)\n",
    "                A[i,j] = A[i,j] - A[i,k] * A[j,k]\n",
    "        A[j,j] = np.sqrt(A[j,j])\n",
    "        for k in range(j+1, n):\n",
    "            A[k,j] = A[k,j] / A[j,j]\n",
    "    return A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Matriz de banda"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se defina una [matriz de banda](https://en.wikipedia.org/wiki/Band_matrix) es una matriz dispersa cuyos elementos en las posiciones $a_{ij}=0$, para todo $|i-j|>\\beta$, donde $\\beta$ es llamada el *ancho de banda* de la matriz $\\mathbf{A}$.\n",
    "\n",
    "<p float=\"center\">\n",
    "  <img src=\"https://github.com/carlosalvarezh/Analisis_Numerico/blob/master/images/C03_Img16_bandMatrix.PNG?raw=true\" width=\"250\" />\n",
    "</p>\n",
    "\n",
    "<div style=\"text-align: right\"> Fuente: <a>Chapra & Canale, 5a ed.</a> </div>\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "casos especiales de estas matrices serían:\n",
    "\n",
    "- $\\beta=0$: matriz diagonal\n",
    "\n",
    "- $\\beta=1$: matriz tridiagonal, para su resolución es empleado el [algoritmo de Thomas](https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Métodos Iterativos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introducción"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A diferencia de los métodos directos, los [métodos iterativos](https://en.wikipedia.org/wiki/Iterative_method) intentan determinar una solución mediante aproximaciones sucesivas hasta alcanzar una tolerancia determinada. \n",
    "\n",
    "Sea un sistema de $n$ ecuaciones lineales con $n$ incógnitas dado por \n",
    "\n",
    "$$a_{1,1}x_1 + a_{1,2}x_2 + a_{1,3}x_3 +\\ldots +  a_{1,j}x_j + \\ldots + a_{1,n-1}x_{n-1} + a_{1,n}x_n = b_1$$\n",
    "\n",
    "$$a_{2,1}x_1 + a_{2,2}x_2 + a_{2,3}x_3 +\\ldots +  a_{2,j}x_j + \\ldots + a_{2,n-1}x_{n-1} + a_{2,n}x_n  = b_2$$\n",
    "\n",
    "$$a_{3,1}x_1 + a_{3,2}x_2 + a_{3,3}x_3 +\\ldots +  a_{3,j}x_j + \\ldots + a_{3,n-1}x_{n-1} + a_{3,n}x_n  = b_3 $$\n",
    "\n",
    "$$\\vdots$$\n",
    "\n",
    "$$a_{i,1}x_1 + a_{i,2}x_2 + a_{i,3}x_3 +\\ldots +  a_{i,j}x_j + \\ldots + a_{i,n-1}x_{n-1} + a_{i,n}x_n  = b_i$$\n",
    "\n",
    "$$\\vdots$$\n",
    "\n",
    "$$a_{n,1}x_1 + a_{n,2}x_2 + a_{n,3}x_3 +\\ldots +  a_{n,j}x_j + \\ldots + a_{n,n-1}x_{n-1} + a_{n,n}x_n  = b_n$$\n",
    "\n",
    "despejando $x_1$ de la primera ecuación, $x_2$ de la segunda, y así sucesivamente, hasta despejar la variable $x_n$ en la última ecuación, pero presentada de dos formas diferentes:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Método de *Jacobi*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El método de [*Jacobi*](https://en.wikipedia.org/wiki/Jacobi_method) es un método iterativo para resolver sistemas de ecuaciones lineales en los que la matriz $\\mathbf{A}$ es [estrictamente diagonalmente dominate](https://en.wikipedia.org/wiki/Diagonally_dominant_matrix). En cada ecuación se despeja la variable $x_i$ que corresponde a la $i$-ésima ecuación, de la siguiente forma:\n",
    "\n",
    "\n",
    "$$x_1 = \\frac{b_1 - (a_{1,2}x_2 + a_{1,3}x_3 +\\ldots +  a_{1,j}x_j + \\ldots + a_{1,n-1}x_{n-1} + a_{1,n}x_n)}{a_{11}}$$\n",
    "\n",
    "$$x_2 = \\frac{b_2 - (a_{2,1}x_1 + a_{2,3}x_3 +\\ldots +  a_{2,j}x_j + \\ldots + a_{2,n-1}x_{n-1} + a_{2,n}x_n)}{a_{22}}$$\n",
    "\n",
    "$$x_3 = \\frac{b_3 - (a_{3,1}x_1 + a_{3,2}x_2 +\\ldots +  a_{3,j}x_j + \\ldots + a_{3,n-1}x_{n-1} + a_{3,n}x_n)}{a_{33}}$$\n",
    "\n",
    "$$\\vdots$$\n",
    "\n",
    "$$x_i = \\frac{b_i - (a_{i,1}x_1 + a_{i,2}x_2 +\\ldots + a_{i,n-1}x_{n-1} + a_{i,n}x_n)}{a_{ii}}$$\n",
    "\n",
    "$$\\vdots$$\n",
    "\n",
    "\n",
    "$$x_n = \\frac{b_n - (a_{n,1}x_1 + a_{n,2}x_2 +\\ldots + a_{n,j}x_j + \\ldots + a_{n,n-1}x_{n-1} + a_{n,n}x_n)}{a_{nn}}$$\n",
    "\n",
    "Genaralizando, el anterior despeje nos conduce al esquema:\n",
    "\n",
    "$$x_i^{(k)}=\\frac{1}{a_{ii}} \\left(b_i-\\sum \\limits_{\\substack{j=1 \\\\ j\\ne i}}^n a_{i,j}x_j^{(k-1)} \\right) \\text{, con } i=1,2,3,..., n$$\n",
    "\n",
    "donde el superíndice $k$ indica la iteración. \n",
    "\n",
    "En el método de *Jacobi*, para encontrar el valor de cada $x_i^{k}$ se usan los valores de $x$ calculados en la iteración anterior, $k-1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Algoritmo del método de Jacobi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ***Pseudocódigo***\n",
    "<div style=\"background-color:rgba(0, 0, 0, 0.0470588); padding:10px 0;font-family:monospace;\">\n",
    "<font color = \"blue\">dados un vector inicial $x^{(0)}$ y un número máximo de iteraciones, $kmax$</font><br>\n",
    "<font color = \"blue\">$k=1$</font><br>\n",
    "mientras <font color = \"red\">$\\|x\\|>tol$  o $k<kmax$</font> haga:<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp; desde <font color = \"red\">$i=1$ hasta $n$ </font>haga:<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <font color = \"blue\">$suma = 0$</font><br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; desde <font color = \"red\">$j=1$ hasta $n$ </font>haga:<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; si <font color = \"red\">$j \\neq i$ </font>entonces:<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <font color = \"blue\">$suma = suma + a_{i,j}x_j^{(k)}$</font><br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; fin si <br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; fin desde <br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <font color = \"blue\">$x_i^{(k+1)} = \\frac{1}{a_{i,i}}(b_i - suma)$</font><br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp; fin desde <br>     \n",
    "&nbsp;&nbsp;&nbsp;&nbsp; <font color = \"blue\">$k=k+1$</font><br>\n",
    "fin mientras <br>     \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import norm\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def jacobi(a, b, tol, kmax):\n",
    "    \n",
    "    n = a.shape[0]\n",
    "    x = np.zeros(n)\n",
    "    k = 1\n",
    "\n",
    "    while k < kmax:\n",
    "        for i in range(n):\n",
    "            suma = 0\n",
    "            for j in range(n):\n",
    "                if j != i:\n",
    "                    suma += a[i,j]*x[j]\n",
    "            x[i] = (b[i] - suma) / a[i,i]        \n",
    "        k += 1\n",
    "    print(\"Número de iteraciones: \", k)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "x = jacobi(A, b, 0.000001, 500)\n",
    "\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert alert-success\">\n",
    "$\\color{red}{\\textbf{Actividad para ser realizada por el estudiante:}}$ Implementar la parte de código correspondiente a la convergencia del método.\n",
    "</div>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Método de Gauss - Seidel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Otra forma de presentar el despeje descrito anteriormente es escribir por separado el bloque de los elementos a la izquierda y a la derecha del elemento de la diagonal principal. Esto conduce al esquema iterativo llamado [*Gauss-Seidel*](https://en.wikipedia.org/wiki/Gauss%E2%80%93Seidel_method):\n",
    "\n",
    "$$x_1 = \\frac{b_1}{a_{1,1}}-\\frac{1}{a_{1,1}}\\left[a_{1,2}x_2 + a_{1,3}x_3 +\\ldots +  a_{1,j}x_j + \\ldots + a_{1,n-1}x_{n-1} + a_{1,n}x_n \\right]$$\n",
    "\n",
    "$$x_2 = -\\frac{1}{a_{2,2}}\\left[a_{2,1}x_1 \\right]+\\frac{b_2}{a_{2,2}}-\\frac{1}{a_{2,2}}\\left[a_{2,3}x_3 +\\ldots +  a_{2,j}x_j + \\ldots + a_{2,n-1}x_{n-1} + a_{2,n}x_n \\right]$$\n",
    "\n",
    "$$x_3 = -\\frac{1}{a_{3,3}}\\left[a_{3,1}x_1 + a_{3,2}x_2 \\right]+\\frac{b_3}{a_{3,3}}-\\frac{1}{a_{3,3}}\\left[\\ldots + a_{3,j}x_j + \\ldots + a_{3,n-1}x_{n-1} + a_{3,n}x_n \\right]$$\n",
    "\n",
    "$$\\vdots$$\n",
    "\n",
    "$$x_i = -\\frac{1}{a_{i,i}}\\left[a_{i,1}x_1 + a_{i,2}x_2 + \\ldots \\right]+\\frac{b_i}{a_{i,i}}-\\frac{1}{a_{i,i}}\\left[\\ldots + a_{i,n-1}x_{n-1} + a_{i,n}x_n \\right]$$\n",
    "\n",
    "$$\\vdots$$\n",
    "\n",
    "$$x_n = -\\frac{1}{a_{n,n}}\\left[a_{n,1}x_1 + a_{n,2}x_2 + \\ldots + a_{i,n-1}x_{n-1}\\right]+\\frac{b_n}{a_{n,n}}$$\n",
    "\n",
    "Genaralizando, el anterior despeje nos conduce al siguiente esquema:\n",
    "\n",
    "$$x_i^{(k)}=-\\frac{1}{a_{ii}} \\left[\\sum \\limits_{j=1}^{i-1} a_{i,j}x_j^{(k)} \\right]+\\frac{b_i}{a_{i,i}}-\\frac{1}{a_{ii}} \\left[\\sum \\limits_{j=i+1}^{n} a_{i,j}x_j^{(k-1)} \\right] \\text{, con } i=1,2,3,..., n$$\n",
    "\n",
    "En el método de *Gauss-Seidel* se usan los valores ya calculados en la iteración actual. Esto, en principio, representa una mejora la método de *Jacobi*.\n",
    "\n",
    "El método de *Gauss-Seidel* no requiere que la matriz de coeficientes, $\\mathbf{A}$, sea estrictamente diagonal dominante para convergir."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Forma matricial de los métodos iterativos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "la matriz de coeficientes $\\mathbf{A}$ se puede expresar como:\n",
    "\n",
    "$$\\mathbf{A=D-L-U}$$\n",
    "\n",
    "donde, \n",
    "\n",
    "- $\\mathbf{D}$: es una matriz que contiene únicamente los elementos de la diagonal principal $\\mathbf{A}$\n",
    "\n",
    "- $\\mathbf{L}$: contiene los inversos aditivos de los elementos que están por debajo de la diagonal principal de $\\mathbf{A}$, y los demás elementos cero, $0$.\n",
    "\n",
    "- $\\mathbf{U}$: contiene los inversos aditivos de los elementos que están por encima de la diagonal principal de $\\mathbf{A}$, y los demás elementos cero, $0$.\n",
    "\n",
    "Para una matriz de $4 \\times 4$, por ejemplo:\n",
    "\n",
    "\\begin{align*}D=\n",
    "\\left[\\begin{array}{cccc}\n",
    "  a_{1,1} & 0       & 0       & 0\\\\\n",
    "  0       & a_{2,2} & 0       & 0\\\\\n",
    "  0       & 0       & a_{3,3} & 0\\\\\n",
    "  0       & 0       & 0       & a_{4,4} \\\\\n",
    "  \\end{array}\\right]\n",
    "\\end{align*}\n",
    "\n",
    "\\begin{align*}L=\n",
    "\\left[\\begin{array}{cccc}\n",
    "  0        & 0        & 0       & 0\\\\\n",
    "  -a_{2,1} & 0        & 0       & 0\\\\\n",
    "  -a_{3,1} & -a_{3,2} & 0       & 0\\\\\n",
    "  -a_{4,1} & -a_{4,2} & -a_{4,3}& 0\\\\\n",
    "  \\end{array}\\right]\n",
    "\\end{align*}\n",
    "\n",
    "\\begin{align*}U=\n",
    "\\left[\\begin{array}{cccc}\n",
    "  0 & -a_{1,2}& -a_{1,3} & -a_{1,4}\\\\\n",
    "  0 & 0       & -a_{2,3} & -a_{2,4}\\\\\n",
    "  0 & 0       & 0        & -a_{3,4}\\\\\n",
    "  0 & 0       & 0       & 0\\\\\n",
    "  \\end{array}\\right]\n",
    "\\end{align*}\n",
    "\n",
    "Se debe cumplir que todos los elementos en diagonal principal de $\\mathbf{A}$ deben ser diferentes de 0. Con esta condición, $\\mathbf{D}$ es no  singular y por lo tanto $\\mathbf{D}^{−1}$ existe. El sistema de ecuaciones $\\mathbf{Ax=b}$ puede transformarse de la siguiente manera:\n",
    "\n",
    "|Jacobi||Gauss-Seidel|\n",
    "|:-----:|:-:|:----------:|\n",
    "|$$\\boldsymbol{Ax}=\\boldsymbol{b}$$|$$\\quad$$|$$\\boldsymbol{Ax}=\\boldsymbol{b}$$|\n",
    "|$$(\\boldsymbol{D-L-U})\\boldsymbol{x}=\\boldsymbol{b}$$|$$\\quad$$|$$(\\boldsymbol{D-L-U}) \\boldsymbol{x} = \\boldsymbol{b}$$|\n",
    "|$$\\boldsymbol{Dx}-(\\boldsymbol{L+U})\\boldsymbol{x}=\\boldsymbol{b}$$|$$\\quad$$|$$(\\boldsymbol{D-L})\\boldsymbol{x}-\\boldsymbol{U}\\boldsymbol{x}=\\boldsymbol{b}$$|\n",
    "|$$\\boldsymbol{Dx}=(\\boldsymbol{L+U})\\boldsymbol{x}+\\boldsymbol{b}$$|$$\\quad$$|$$(\\boldsymbol{D-L})\\boldsymbol{x} =\\boldsymbol{U}\\boldsymbol{x}+\\boldsymbol{b}$$|\n",
    "|$$\\boldsymbol{x}=\\boldsymbol{D}^{-1}(\\boldsymbol{L+U})\\boldsymbol{x} +\\boldsymbol{D}^{-1} \\boldsymbol{b}$$ |$$\\quad$$|$$\\boldsymbol{x}=(\\boldsymbol{D-L})^{-1}\\boldsymbol{U}\\boldsymbol{x}+(\\boldsymbol{D-L})^{-1} \\boldsymbol{b}$$|\n",
    "\n",
    "Ambas expresiones presentan la forma $\\mathbf{x=Tx+C}$. Se puede observar que esta expresión asemeja a la iteración de punto fijo, con:\n",
    "\n",
    "$$\\mathbf{x}^{(0)} \\hspace{3.5cm} \\text{aproximación inicial}$$\n",
    "$$\\mathbf{x}^{(k)}=\\mathbf{Tx}^{(k-1)}+\\mathbf{C} \\hspace{4.5cm} \\text{ }$$\n",
    "\n",
    "donde,\n",
    "\n",
    "- $\\mathbf{T}$: Matriz de iteración\n",
    "\n",
    "- $\\mathbf{C}$: vector constante\n",
    "\n",
    "- $\\mathbf{x}^{(k)}$: $k$-ésima aproximación del vector solución $\\mathbf{x}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convergencia de los métodos iterativos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "La convergencia de los métodos iterativos dependen del tipo de matriz de coeficientes. A continuación se presentarán algunos teoremos que se tienen qué cumplir para garantizar la convergencia:\n",
    "\n",
    "- ***Teorema:*** Para cualquier $𝒙^{(0)} \\in \\mathbb{R}$, la sucesión $\\left \\{ \\mathbf{x}^{(𝑘)} \\right \\}_{𝑘=0}^{\\infty}$ generada por la aplicación de la ecuación de iteración $𝒙^{(𝑘)}=\\mathbf{Tx}+\\mathbf{c} \\forall k \\ge 1$, converge a la solución única de $𝒙=\\mathbf{Tx}+\\mathbf{c}$ si y solo si $\\rho(\\mathbf{T})<1$.\n",
    "\n",
    "\n",
    "- ***Teorema:***  Si $\\mathbf{A}$ es [estrictamente Diagonalmente Dominante](https://en.wikipedia.org/wiki/Diagonally_dominant_matrix), entonces tanto el método iterativo de Jacobi como el de Gauss – Seidel, generan sucesiones $\\left \\{ \\mathbf{x}^{(𝑘)} \\right \\}_{𝑘=0}^{\\infty}$ que convergen a la única solución del sistema $\\mathbf{Ax=b}$ para cualquier aproximación inicial $\\mathbf{x}^{(0)}$.\n",
    " \n",
    "\n",
    "- ***Teorema:***  Si $\\|\\mathbf{T}\\|<1$ para alguna norma matricial y $\\mathbf{C}$ es un vector cualquiera, entonces para cualquier $\\mathbf{x}^{(0)}$ la sucesión $\\left \\{ \\mathbf{x}^{(𝑘)} \\right \\}_{𝑘=0}^{\\infty}$ generada por $\\mathbf{x}^{(𝑘)}=\\mathbf{Tx}+\\mathbf{C}$ converge a un $\\mathbf{x}_v \\in \\mathbb{R}^n$, que es la solución única de la ecuación $\\mathbf{x=Tx+c}$.\n",
    "\n",
    "Además se cumplen los siguientes estimadores (cotas) para el error de la aproximación actual con respecto a la solución del sistema:\n",
    "\n",
    "$$\\|\\mathbf{x}_v-\\mathbf{x}^{(k)}\\| \\leq \\|\\mathbf{T}\\|^{(k)} \\|\\mathbf{x}_v-\\mathbf{x}^{(0)}\\| \\hspace{1.5cm} \\|\\mathbf{x}_v-\\mathbf{x}^{(k)}\\| \\leq \\frac{\\|\\mathbf{T}\\|^{(k)} \\|}{1-\\|\\mathbf{T}\\| \\|} \\mathbf{x}^{(1)}-\\mathbf{x}^{(0)}\\|$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Método SOR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una variación del método de *G-S* es el método de [sobre-relajación sucesiva (SOR)](https://en.wikipedia.org/wiki/Successive_over-relaxation). Después de que se calcula un nuevo valor de $\\mathbf{x}$ por medio del esquema de *G–S*, ese valor se modifica por un promedio ponderado de los resultados de las iteraciones anterior y actual:\n",
    "\n",
    "$$\\mathbf{x}_i^{(nuevo)}=\\lambda \\mathbf{x}_i^{(nuevo)}+(1-\\lambda)\\mathbf{x}_i^{(anterior)}$$\n",
    "\n",
    "donde $\\lambda$ es un factor de ponderación, $0\\leq \\lambda \\leq 2$\n",
    "\n",
    "- si $\\lambda =1$, se obtiene el método de *G-S*\n",
    "\n",
    "- si $0 \\leq \\lambda < 1$, se tiene el promedio ponderado de los resultados actuales y previos, llamado subrelajación\n",
    "\n",
    "- Si $1 < \\lambda \\leq 2$, ponderación extra sobre el valor actual. Suposición implícita de que el nuevo valor se mueve en la dirección correcta, pero con una velocidad demasiado lenta."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Algoritmo método SOR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- ***Pseudocódigo***\n",
    "<div style=\"background-color:rgba(0, 0, 0, 0.0470588); padding:10px 0;font-family:monospace;\">\n",
    "<font color = \"blue\">dados un vector inicial $x^{(0)}$, un número máximo de iteraciones, $kmax$, una tolerancia, tol, y un valor adecuado del parámetro $\\lambda$</font><br>\n",
    "<font color = \"blue\">$k=1$</font><br>\n",
    "mientras <font color = \"red\">$\\|x\\|>tol$  o $k<kmax$</font> haga:<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp; desde <font color = \"red\">$i=1$ hasta $n$ </font>haga:<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <font color = \"blue\">$suma = 0$</font><br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; desde <font color = \"red\">$j=1$ hasta $n$ </font>haga:<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; si <font color = \"red\">$j \\neq i$ </font>entonces:<br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <font color = \"blue\">$suma = suma + a_{i,j}x_j^{(k)}$</font><br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; fin si <br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; fin desde <br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <font color = \"blue\">$x_i^{(k+1)} = (1-\\lambda)x_i+\\frac{\\lambda}{a_{i,i}}(b_i - suma)$</font><br>\n",
    "&nbsp;&nbsp;&nbsp;&nbsp; fin desde <br>     \n",
    "&nbsp;&nbsp;&nbsp;&nbsp; <font color = \"blue\">$k=k+1$</font><br>\n",
    "fin mientras <br>     \n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Volver a la Tabla de Contenido](#TOC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open('./nb_style.css', 'r').read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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": "428px",
    "left": "91px",
    "top": "139.708px",
    "width": "530.662px"
   },
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}