A.1 Solution of linear equations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical algorithms by examples\n",
    "\n",
    "## Contents\n",
    "   - [LU decomposition](#sec1)\n",
    "   - [Jacobi method](#sec2)\n",
    "   - [Gauss-Seidel method](#sec3)\n",
    "   - [SOR method](#sec4)\n",
    "   - [Thomas algorithm](#sec5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as scp\n",
    "from scipy.linalg import solve_triangular\n",
    "from functions.Solvers import Thomas\n",
    "\n",
    "from IPython.display import display\n",
    "import sympy; sympy.init_printing()\n",
    "\n",
    "def display_matrix(m):\n",
    "    display(sympy.Matrix(m))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem we want to solve in this notebook is very simple:\n",
    "$$ A x = b$$\n",
    "where A is a $n\\times n$ square matrix and b is a vector of dimension $n$. **Our goal is to find the $x$!**\n",
    "\n",
    "It is well known that the solution of this problem is simply:\n",
    "$$ x = A^{-1} b $$\n",
    "however, inverting a matrix is not (in general) a good idea.  \n",
    "Matrix inversion is a slow operation and there are plenty of algorithms that permit to solve the matrix equation with no matrix inversion.\n",
    "\n",
    "I will review the main methods by examples.   \n",
    "Let us choose the following values of $A$,$x$ and $b$: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matrix A:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHgAAABkCAMAAAB6rgQAAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRCK7ie/dzWZ8bFsxDvgAAAAJcEhZcwAADsQAAA7EAZUrDhsAAARwSURBVGgF7VrtcqMwDDQ2ob3wfef3f9aTREyQQNK0TJjONP6RgBd7LdnAshCqTCWGi8q48IVQ5VRDuV3EG2ZkixmJq6s4nzwzI25T6jp9FFPdhCqpmWnquk/3Z998qxpj6rEs1Yy4wWnus9p2wMl5NOS90l7C36iNa15mNj+WEiOODTYdBvw9KlMf6/YIoLrbjH/tpByQlpZlCTPiYcRGUZ3y0ui4677D+ob6ODhiyVRfJpIRTxlD/i7xPXfQvKe4D4ip6l4XhBEvlSPRlwO2/3Ge61SGvAWW7ZjzbPMGSgodvSe+68sHI2r1tQfrMk+0TPaDWmo2w9oTj/ZMhk6bxHCL9zEP6jkB8w/nbik74kQnRYEP/lNWFvaMI67zoMdcb1a8JK4N3kjt1PN8oAFVWTuRQxifUxwE8Q15WyWkgYiTsvZKHqO6rJty8cA8cuKK4tUuEksyxk2+sIe1LBGHThk23YzWgzlxO8KlNEVt9bTIPD/6f3ZRtu60ouf1TC316/8tL0OnChbxtFxONeLQpj4mfe1UcBfQbxIwhcNm+hnxOrYLNt7EFyR5oXin+p3ql2Xgty8uR7E7et+Bud7nqbYVu6P3HTjQrWnV+5zYVuyO3ndgofc5sa3zHL3vwELvf4XY0fsOLPS+ILYVO53Uut53YK73ObGn2KFrQ+8jsQUzvc+Jacy6YkfY0fsGzPX+AbGq2JHX0fsGLPQ+En/kD+wTnhNtxQ5HWHrfgYXe/7u1ImzFDh2bet+Gpd5nqXYUO/gfQB40ve/AQu8zYkexe3rffhwQep8RO4rd0fsOHLje58SYyYvKm/iiRMvn48tofytxA6bs5aUF5+AH+NUXxv0zLiDCRN/F7+h9sG77evWFd605zCIWJvquqa33QzOBj5k0G0zAjFiY6DtiW++HDm/X3cY25B1wmBELE523gz1b7981d3XpR8CMmI54muhfJO7UdxnUkYD3xGqqIGJT7+ep6mvd2RPwjnhjou8iNvV+kyf0bbe24bYDCUviIga3bfi2pvebTH5yUoxyCUvirYnOCcueqvczma+z9s5CwEj88fmn9MpM9LXyseHo/cXInrNilAv43ye7STATXRI7ej/SqlYjFrBINbzRlXTPfUfvV485fjZgWwIWxMxEZ+1gx9H7IcIFqNFWNZyMDBbEzESXxI7eh5tASlF/EcdhQbzjelnFm/hlqZUdv1MtM/Ky/XeqX5Za2fEPSbWr2OEeaty+ICoD/r5DT9navm6X6YN9A6YRf8uhR556MFSoCZ9w6IG3nSeL2IJPOPRA3AeT2IJPOfRzaxLb8AmHvgH7wIjYgfkXOQfnsabYMdEmsQN/36G/oZ7SI3bgA4d+NV9sxd5gRDqxAwfh0LNU24r9FrHkIR5/YuLA5aGsfJHDiB3FjgEH+wJiwCcceiLO1gUkBB0+4dADcd/lXJJF4+A/Jvx26NljKs/cK/fYqn4lkez7txLTl2u2bydTdWJ//aIcVBgW40vOEywHTemL8roO/wHWgUEZoFrzzAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}2 & 5 & 8 & 7\\\\5 & 2 & 2 & 8\\\\7 & 5 & 6 & 6\\\\5 & 4 & 4 & 8\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡2  5  8  7⎤\n",
       "⎢          ⎥\n",
       "⎢5  2  2  8⎥\n",
       "⎢          ⎥\n",
       "⎢7  5  6  6⎥\n",
       "⎢          ⎥\n",
       "⎣5  4  4  8⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Vector x:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABoAAABkCAMAAACvvNBNAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRO/NZiK7id18bOsX9RIAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAGRSURBVEgN7ZbbcoQgDIaR025FEdq8/7OWEHCNhJledKa92Fw4ykci/CRRtUA1rU7LNKLUAsYWcydROz5rQLS8Rl93+8+RC777MS+vjQYZlfn2jU7ljQ6w6Z1UZBp2Yd+IK4FPXKhozHH0VGbIY0UkWCkCQ7qmZwgCChkHdTsj5rXV1BURBcott5lXRSskmjKi3At3QMaQD235AY/+qOxJ1CcvWIckxjqTB1yqjxVQzCklo+vGb/Ju1EkkdC5GeBdjfBl/irw9+9J1HdH+Rj+8RORbjkXeo/dghiKeiQOpA1hAD5CUX0PJaQ8bvZAFrENyQEQ+tzS8FZGKKZsuzRgw4QvRRqRC+wogejw/aB4l4FG3oNTX8yIv1BrpiAUMB3pnEFJ0R4kclHNCY15qNaUFtAXeEE1vV+71RpdDqWJQgyi3o1C5noyELB2agOK+zbySmqE9zpAvZz/xwpyWkcMPg4h8LQQROY0GQVO2jWpMt1zqdbrlA6B96ceAZZFk/wpN/yo9/kNa20uuLL3+VVqrvgFO7SD82wBfAQAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}1\\\\2\\\\3\\\\4\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡1⎤\n",
       "⎢ ⎥\n",
       "⎢2⎥\n",
       "⎢ ⎥\n",
       "⎢3⎥\n",
       "⎢ ⎥\n",
       "⎣4⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Vector b:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACUAAABkCAMAAAA2RAn4AAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRInNuyLvZt18bE31V9wAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAKKSURBVEgN7VfRkp0gDI2A3O0FFbv8/7c2CQSNgjPtS7ed5cF1b44kxJycCFPmZaC39mIEmLJ1uOYeCBKZTCbU1AUcP6YLKgQX5JnJCk6jfHwD2FiN+9JHLfT0Uo1u66PeeZWnAdYU+6hlO0AQYIDKcQrOYmi40jpA+RwT2jdKnHcwRGWPALvhNeDNyOOONkj5DTPlbIDaOVMpJ09bjVCGz4h7zYZW3gxGR0vlfqpxFRMMsgoGPXk+IwFzP6t4NGtNfdlhydlQZnApj+WnzvUbJUn5Cplgrk67sYFWieweF3M1lR6Sa4e5oQpXbSGmtKErqnK1eGod44o6c/Vd6/lWX4qrrVQvVai4Gmqh4jG1RwpHWOixs8lSKMVVJ20MoYR65Rc/pLl69DiAn6eOqbjqJaO0gfLIO1auYuPmf/lyR1Wuzk+oxtW1Ubzr8fBz3N09Hrbj7t9HeUcN/GmtWI+/r4/9Hf8kXxEjnCzrEJLWNB3Ue21EViprv2HRr5vANCoG49hSVMRIiWmUcBkyF30QgvRRvpSikxHigkrJWdKOre5VdPdS9wu+hBXlEYqvmMcKs5CWRnwnc8h1qtEeOfOW5wWcb6Zwjuv18YPNYHQ0YFgv0fb5cXrbJWYrNoiSGeWxJJHUO1HapzaqKNRKMAa46MHzdMGxKBSsNhjLJYmiW4cUgmkUP9i5fKMkKX87E421j1rbWPuotY21j1ortcmDKDJJpnOdiYbiRI201ghrGTXS2sZaQg21tmzBEzAMtZZBUFgLV61t3Vex9qy16oxn1o619mAtUVYaIcah9jpYiw1niDqxFlvv6VtL7VXO2Ln+DyieFXVxyUnbt6inz03nqgKIuf7lb1Hn4Bcb4ytHWSx1cQAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}64\\\\47\\\\59\\\\57\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡64⎤\n",
       "⎢  ⎥\n",
       "⎢47⎥\n",
       "⎢  ⎥\n",
       "⎢59⎥\n",
       "⎢  ⎥\n",
       "⎣57⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])\n",
    "x = np.array([[1],[2],[3],[4]])\n",
    "b = A @ x \n",
    "print(\"Matrix A:\"); display_matrix(A)\n",
    "print(\"Vector x:\"); display_matrix(x)\n",
    "print(\"Vector b:\"); display_matrix(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix A of this example is invertible!  \n",
    "We can check it is a full rank matrix and calculate (just for curiosity) its determinant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rank: 4\n",
      "Determinant: 194.0\n"
     ]
    }
   ],
   "source": [
    "print(\"Rank:\", np.linalg.matrix_rank(A))\n",
    "print(\"Determinant:\", np.linalg.det(A).round(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first method I'm presenting, is forbidden!!   \n",
    "So... read it and then forget it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The inverse of A is:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAABkCAMAAABqx9+wAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRM3diWYi77t8bNSLqrQAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAn1SURBVHgB7V3RoqMoDLVqndlatd31/791CRECIcEw8zZDXy7BcAiniJajucPj9J9x6B/PwBv5GIbHOc3u8+y8IAM7sDGeQMyjc8IZ2KvELNO2bwlrZC/zvE0vjpba5JvWQtkd2caD14o2Bzmm6fv1Aa3zMjymX5zlBey2bd8cq07M6pyXN42B7AnGMeZY+dDIN9Q/PwsUl88+DMeHUMNx4S8DWWAd3E74Qj6wEGxCE0sVgz1gNM/ThUWfKjH7Bxy3NbiT/fQgRzwSPOgv+WLdMk7j6YkZPeroqSV/scRBRg/wAYB1G2cTtwIwh51P+IbPd+paJebrL1RPHI9rRfb2BZAlg0phU99YPyPQ6enc3Dl8+6EO0fXjOxxhTfydayiHfX3cHFwwrhBTlZjT9/7yfEIDsl/n1315Wzb5AiT+Jd9Yj8Qspyd1tqz4HGT11P42MRzWR9hwKi2nn+6vc8axpfZ4nnuNl9Q3MHPNmM81Y6pLt28jgbgDb6Bn3Pd5Sq4LoRPDXxF2eecLVm3GHLi4PcIal9nbea7+jJcjyXwvl4sYPInWfK0TUSSQYXj5eGDGHn4ZFpvWKgXYY3tP+WjqxPgZkxBD9nN8vU84NZXPgbMttgW3i5hhdTPwucUTVEFw1RKImzC0vHwri1wr7JaPJiNmWd/xsz7ccqSeSjsEN594/ZVC4G3BJxAzuFvLx2ZYYySQYfJBYZ/Tabow3Yzriv+TjSYj5vKIf3CResbvlmy8C3nEI7FJLJBvrIrEQM117Y4HxYIEgryM5oVKAOawh6f3m42mSszXXz92ulwHe7mutaN+WeJtXXwZMSudEELoV1UJ8gRe3EBwBZ9iaDqGcITDnh6mgRi8EZribRzZ133rV5/J5Bsju4jZofHDdBIUIA8/X9ydHc6bdwwtdmIpcNiP/8LfWUjVGTO84SfBx10Ur7ufaL/8FWm/ruNiMNE33jldxMyu7bLqcy1FYyDH2/2omUa35Pq7eE9x6m4tM9gdgnmGuxIEqROzbO73t79ZwEsB2Y9xuvsRydpO4+dc/bm3TVP99yeNjzr0Aay4UwLXomPaRnaFpWZ3JQY7vFxE3/wSWyfmroM/+HgnRvlyOzGdGIUBpbrPmE6MwoBS3WdMJ0ZhQKnuM6YTozCgVPcZ80vEFMJU0Lseb/dbCT4KLK/mQM62am4ExUFIfSMfc4kay9phfcYwYYr0rh1/zeGOjyEYDtSguRE6B4EdHVTfyKdSuvQ+9KCh4A5GoR1WicF9CxLcSO+acCPGstkEgRRAds2NRlqA+N1rr76Rk1YivQ89aCiKdlglhgtTpHfhKZTK2lpEvp4DNWhuhMtBKBryqZWyDURqrGiHVWLwTCHBLdG7IIJXbZ8qC5EBtWhuhMNABhYNOSqljBhqrGiHNWLETfpL74LO/YagEkVWXQDZNTfCKUDwkFffyKtSyohJGsvaYY0YQZhyeKh3uUJNiMzjK4DsmhsBFSD+UIyGHLVSSczVWNQO68T4LedMNCO9KygFWhxJfSmcmTU3QilB4FiivpGrXCqJwcaydpgRYxGmot41V3bo74EMmts9yJCpbzIdVFsQg0NRtEMg5sf5g5qnJVzuSHBzx2bULVzpbV5iwmMSGZBDMGluFE81GnJTS5yYayiKdvhv7VEzLkw5jQF48cLdgnGqYWQHSiA4bNLcCKcEoWjISy8xYq7GYUXg2mF2KnFQLky5h97AxT/J5B715O66zYEaNDcC5SBpNOSll3Ji4lAU7bBKDBfcSO8CfaqBGA7UornRSJlMlkZDTnrpIgb1P2qsaId1Ypgwlehd7unC2pOJPDwGNDRobgTFQNJoyEkrkd7HpTtZO6wTo/XyF9R3YpQvuRPTiVEYUKr7jOnEKAwo1X3GVIhZ3Isc/cMYONwuXH9fiZECZj+VBFI6MQopOjGqtNUqtMHbbNk7ctAn/rCtBMUOcQzSypijxWRgyYBy4U0+lVRpq1VoGxiSD71liwsaMIxEK7MwwXwYWDIgv1sQhTeRGF3aahTaCqUNopzxeWMWsG7q0eht1CMcjAbEhDeRGF3aahTaknfiYqjHvjZsirpmejQR017gYDQgJryJxNxIW3ahLXknLsa+DY3E3EQTgU0FDuYb+QEx4U0i5k7asn/hAtJ+NBIjYMBo7EJbSpgMhgPKhTeJmBtpyy60uTex/FRNlKnF3VC2zZgSA0baILSlxIhgYUCZ8CYT4xfoZDiAHaStsK2e9qeVS5EMiGolphaN1rNSXwYEL5Cgcy68ScSI880ktPFwCqQnvLHRRkyBAZ3EaHiP3DaodpdyyIQ3IObHz39yPFyhMn3s14S2a/GNSIs/s9qI4RgQKkWTB35vCUO7bquY8PbfT+FHZE3aahLaBob0HOFzfkbzAyRuqAzD1bQJbRlbJdg1oHBCBeFNOpWu2zJ6sy2VtpqEthIJovylGzw5mmzQFqNQ7SBJjG/IhDeRGK6PkTrVKrRxJB8CvqpvGQb6/KbQlnfEwGhATHiTialIW21C28CQXJDb9zzDfM1j1iyG0Sa0cVAG5pT4oBzmwptMDEf7C+1OjPKld2I6MQoDSnWfMZ0YhQGlus+YTozCgFLdZ0wnRmFAqe4zpokYpkolspk9kyPHoAAO92YcS+BIB+WSDmYNiCOk9rbNxRtG8oxhqlQiedkzOXIM9zsWk0BJCRxlOqhWBTOnluQIZC+rS5CS7GpgryIxXJVKbHMmx6SN74leMZMSOBIFYkkHs6aW5AiJ/YX9GNzASjoXieGqVGJbX/YrZDLX5/UEspTAMYlIKiYBhMPhOW9jQByB7FeWgSnAy4+BcFUqsY1xhPcq6O0412MYC3TOEjjGeORCEkBwCGDGgDgC2V//fmaAjX+lGcP35VPbmskxbRM6C2NxNk/gGFzkvzUwW0AcIbHP9bHNZWooiRiuSqW2NZNj2iaMNhJTJnAMLvLfGpgtII5Atnu1AFJ3hW28GIBMTC5xFSqVIZNj0cb1GIlxZZbAMcYjFu7AhtuAOALZCybBm7K0iS4KiZhknvk4uT0YMjkWbRxSSsyQJ3BkfFhksjQ34G1APJzExvzGO89nKhFTSFy4UoFsZs/kSG3ioC9ipASO0Ucu6GDWgDgC2ZhkcOfJYkViuCpFtj2TI7WJY72IkRI4Rh+5oINZA+IIZGOmbtuMwbsfuhkkG7UpSyZHahPHehEjJXCMPnJBB7MGxBHIfvhzcuIXbXHGFDJZVKkaMjnGNjylpJTAUeaDalUwc0AcIdrD6AT1xXRVcn4sHSTZ9kyO1AYfIKFXzIQEjkSBXNLBrAFxBLLhdTvMm5n2Lc+Y1OMvLXdilC++E9OJURhQqvuM6cQoDCjVOGN8Ki7jzoYC9AdVx//I5V67gE+eGvoPGmjrUPx/5Jrn4X84mZNoXNd3OwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}-0.08 & 0.19 & 0.25 & -0.3\\\\-0.27 & -0.9 & -0.2 & 1.28\\\\0.27 & 0.4 & 0.2 & -0.78\\\\0.05 & 0.13 & -0.15 & 0.06\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡-0.08  0.19  0.25   -0.3 ⎤\n",
       "⎢                         ⎥\n",
       "⎢-0.27  -0.9  -0.2   1.28 ⎥\n",
       "⎢                         ⎥\n",
       "⎢0.27   0.4    0.2   -0.78⎥\n",
       "⎢                         ⎥\n",
       "⎣0.05   0.13  -0.15  0.06 ⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The value of x, given by matrix multiplication 'x = inv(A) @ b' is:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACoAAABkCAMAAADHT1J1AAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRO/NZt2JIrt8bAICnK4AAAAJcEhZcwAADsQAAA7EAZUrDhsAAAJYSURBVFgJ7ZjRlqMgDIYjIDMrCrib93/WJUEEPRLaPXsxF+VCW/mbxiR8RGFCHgq6w2cFwITapDF3lbDRvEKSTn1VndkkqdVhC6eZm3R2tpqBNflkfTyutFKrtMJGujkShfVBmi6ZVrpzTOZyqbV6lyJLFzxCI0gtavrrBQ2dAARpxECKKZ8GUrb6ivQNByDf1vzCbcG+57s5Qi3cFuQUaDkFFnneU2JdKYKLVa0crmpLf+s5+jak0ivKe1zJud64WO2J8vWP9GdFIGq972fOU6l0kWGpRgIuZ377yFBcw45BQXIBGc6TQJ2wE5CxMkeqVEZGsuoLeAYrlmDCoEi/GSCjLBVyOmYQ9ZABmolCSnhy4Au/eC4dTFXCHRm/L4CfSRkLpiVkTGzTFKmAjOhDCFqlRAyRseaNjHL2QUbJXnt+c3Fb02ysrZ3r55i2uv/RD1ytvunr6UBM5bI3LUwfGZFqcEbaYfLoI8PwFolMA9IKyFhcYtBRqiQVkEHTrQMDZFhfMPC4Ys8IQAxen5kbIiOQwzxGyABwpSt7RMb3r2yGV/Veuoo7Mv5813JBZmCVCshw3IB4fAEZG6V0pg5siAxYdNpiOAAfZOQquB7/dRlerdy//TyreUPKfvaRwfOeKzFL+8igeZOLlqUCMtJ83NZqVUZGgEYqImOLjfRpxZ7IsGmtVKsiMohBrZQ38scuY6bGrUoFByyDrUrvyGjKZVY00Knj+WrwYJKQVeMqdBmcIiTpGBmpz9wR6cnkgwyO2+2Qs/XyWwdL7xSMKZvQzRh95bcOxsBfLMkuJFdQyOYAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}1.0\\\\2.0\\\\3.0\\\\4.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡1.0⎤\n",
       "⎢   ⎥\n",
       "⎢2.0⎥\n",
       "⎢   ⎥\n",
       "⎢3.0⎥\n",
       "⎢   ⎥\n",
       "⎣4.0⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "inv_A = scp.linalg.inv(A)\n",
    "print(\"The inverse of A is:\")\n",
    "display_matrix(inv_A.round(2))\n",
    "print(\"The value of x, given by matrix multiplication 'x = inv(A) @ b' is:\") \n",
    "display_matrix((inv_A @ b).round(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec1'></a>\n",
    "## LU decomposition\n",
    "\n",
    "For more info about the theory of the LU decomposition, have a look at the wiki page [link](https://en.wikipedia.org/wiki/LU_decomposition)\n",
    "\n",
    "Let us have a look at the LU decomposition.  \n",
    "The matrix A is decomposed into the product of three matrices:\n",
    "$$ A = P L U $$\n",
    "where L and U are lower and upper triangular matrices, and P is a permutation matrix, which, reorders the rows of L@U. \n",
    "\n",
    "If you are not familiar with these concepts, they will be immediately clear after seeing them in practice.  \n",
    "Scipy offers 3 modules:\n",
    "- [scipy.linalg.lu](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu.html)\n",
    "- [scipy.linalg.lu_factor](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu_factor.html)\n",
    "- [scipy.linalg.lu_solve](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu_solve.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lower triangular matrix:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAABkCAMAAACcqw0BAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRO/NZt2JIrt8bAICnK4AAAAJcEhZcwAADsQAAA7EAZUrDhsAAAgRSURBVHgB7VzdorI4DERAzreKoru8/7Nuf0kyTQro1ceBG5qQmXYoICRgc5nD0ja/ZHlGvU1zmbveLddfort5eLXt7IVffotm0vlYFz5002PatWsQgTZ1v9ZCJNo1PMZKWxV+vQ+McXTnwPB8Mc9aExFor+FpOyLRpsiyhbHSLoUPbdfOTPjj7jmnsWS2PIhA28KVfkSiXSLIg7Fgl8IdtOfC3+Fqf+UuYldbiEBbBalORKKtgpITY8FeFz4H4bd5+0UfEWjXhiu3IRJtGS0tjAV7Vfgwd57wNveS17YQgbaNxC2IRBvjuY2xaK8Kf82T57vEFae22ohA28KVfkSiXSLIg7FobxAeZnyXcIl4xWNmB0MePiLRznHaGmPRXhWOh4jWifQhAm0ZXbMQifY32FXhTbwoXHdf3AixnyErQiTaOU5bYyzY68Lfb0/72PNzBoj9DFkIItHOcdoaY8FeFx5/+LvdNzCE2M+QhSAS7RynrTEW7JrwYQ5qn/6W9b7nZn1BfMyQlXzDVMcqwrv2Po/tw/X9DPcuw+Qe4vbobgjxKUMW/g1THasIz50ee30KP/b8lurOGS/3ybE954wfe35LdeeMl/vk2J5zxo89v6U6OeMy5940r657v8NtutsytXty67wrpOXb6m1Eol1DY6y0pXCZc28G/5AyzTf/dOYeWl73zcq/LUlkQTCgBu0cp60xVtpCeHxkpdpBG+oKd1dRaENVoQ3JNK0T4fu+JJHpcEBo5zhtjbFgC+GQc2/uT8/YuppifDKfXKFt2/JlSSJ3ggNCO8dpa4wFWwiPaSmqHYwh3+SED3PIJvWby6pCONJqw9R9iERbR0UvxoLNhRtJzKeTfw+5mHC61zpbtnHhBu0SazcQibaNdFclKISgzYVjzj3y3nwpIR7k4+zzMlsWLlyn3cKCSLRrHBiLthQuKwGRN2aPRldBuk6bc8xSuEZbG3TehkUAtHOctsZYtLlwPBoCX5eu5O79icv00Tmu0mpjLXyIRLsAMAfGou2F/8w/ERHPf6oEOG/Pf8FE3Zx1Ujb5jH9QksiEOCC0c5y2xliw/+WvgkDO3dFdve5Xvm8ZQ9ZV6wV9QnhJi+GWjUi0LZz3YyzY/FBv4m88VQKaS5jv/uW2OPGXOe+BWodhmxBe0K7CcwAi0c5x2hpjwRbCG8jBv57TNHWtu43px6EZxq3X9OWdir+moAA5+DG+Cefv36au69w9+7bl+5JE7gcGxEoVOcJe17Fyxm2Ww205hR9uSlcEnTO+soMOt/mc8cNN6YogP+NDz15dXYk/zOaXe9w831c/zHRuEHJe3MJOkjn35vJsXR3BLRv2YAoBBgak4gRz7m3a9CUTxkpbzrjMubu3+uISHsNljaDsJ3qAwTkTjooTFlT3y25Leh3lvRgrbSE8PrJSQaGLz99ON9YIrO6QgXBUnLCwip/gcSPSK5DFhbFgC+GQc2/iIZ4+xBG5hYUeGsjgNiccFScAsmKKbhV6E46xYAvhMS1FBYVAektv6IsRWN0pDAm3FCcsqOEX3Sr0BqxJmT4SA1guHBORkTPUUFxTjMDoTmMQOF+c2LdwuEZvsWEs2lw45twD55TzTXwEVm8aA8eF4oQF1v0crtHrKJchhS8r0JbCQ25RfEowLHVCPgK7t5KB42JxwkKrfg7HooAKSE6MRZsLx6PBU/TL69p8BFaPKgMd3bk4YcE1P+9Wo9cw3oexaHvhP3/+ifB4/ouCwjOf4pvOca12QCMXxQlrwOgnuNuiDBDjFxtjwf7vD3tIgZy74xhiuGcTI1jooVEyEE4WJwBomqJbhd4EYizY/FAvCwruwS2ctJ5djMDqLt4lsJIE4ZbihIXV/aJbhV5HOS/Ggi2EY0HB3W/uFF4y5B1GxQlzqOqGJPyD0gRUR2hooR8pHHPw7n2f9GUl1QjU4S1OZFhwrDixBK83FvgHX0vgUMj2/Urh6yM5TMQp/DBTuVHIOeMbd9Rhws4ZP8xUbhRyzvjGHXWYsHPGDzOVG4XIGZc5908KCqlbIAre+Hjmm9RK4ZVVjWma+h3/xSSZpHCZc/+koJAkAFHwUk6DWhXFq0zD6F7DEk/AyFYrRgjh8ZH1m4JC6hqJvLu/52QOtXCkpV1hevtUQUwvlLiyBgJMQjjk3D8oKKQRIJFzvx5jEk4tbbzgs5lu6+9ZiiQGMAnhMc9EOfgwiF0FhTRshWhqsnBqgUjNtJne4TMZDbP4hHBg4sIxERkJlgOUsqULs9FQiB6vLJxaBpq7K0zzeJn6+tuWXDgyceGYcw8j2FVQSGMuiQZXh4ozTi0u0GrbTC4X5UsdOUOkEnDhyCSFl+WAfQWF1D0m792rsG5LFE4tdazgtJmGORyBnfjPPUBL4VIcF45Hg6fZVlAYxueyjJcimd9c/deKQTi1YJDJ3M40h2/DHv5zQGvhwlEcF67l6+knl7NYPWV/vJAslYkh1Ju9cGrl0JW1ydQ8Q43nUfs+SAwZmIRwyLm7Me0tKCQZQHRt/TLf255aK4LzZpMpffy4ecbxiwUhPP7G85uhvQWFNN6SyG1YbmBYK+sz1zbTJZ3jJpRqGSEEmIRwyrmn/P3+gkIaBibzvTt+rShbKbyyspladwINW67qajFCCqeceyro7i4oJAkFUTO95zn8jxJrVfQumypM7quJ2l8x1YsRUvjS3fEbp/Djz7FUeM643B/Ht375jIdXVjd/OPq3Hw7LP+oP/h/m+75yt/+3S5XjD/+o3/fN/7JjgV+jf1/TAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}1.0 & 0.0 & 0.0 & 0.0\\\\0.29 & 1.0 & 0.0 & 0.0\\\\0.71 & 0.12 & 1.0 & 0.0\\\\0.71 & -0.44 & -0.46 & 1.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡1.0    0.0    0.0   0.0⎤\n",
       "⎢                       ⎥\n",
       "⎢0.29   1.0    0.0   0.0⎥\n",
       "⎢                       ⎥\n",
       "⎢0.71  0.12    1.0   0.0⎥\n",
       "⎢                       ⎥\n",
       "⎣0.71  -0.44  -0.46  1.0⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Upper triangular matrix:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOgAAABkCAMAAAC7BYzpAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRN2JZiK7ze98bHWeL8UAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAjQSURBVHgB7Vxd26I4DEVAZlY+hF3+/2/dtKFpchoUZx4v5JULbUpO0kNLgUO1uqxxq6uTbj3zq6rL2rS0XU/KsxoDu3oNRC9n5Zh5jZropa+bIWxpf9cM4/DkKMxtV12aPByOYFL89D0MrUrzlgiG6MiDeZWzdab2d/2U2uN+LwEkh6aqjmBsoG6+VVUzS+VbIhiiDVMSnuMSkg+5CdIWVZiHulWH4hBGwal4b8LHPVW+J4Ihyv2SB9E9Ur6uXWqD9y2HhXcewpg4t1UdJ9rzngiGaMx/a6UZPIZvaz4BZVcuANFDmIym0j2Om1z1ngglURlCVbeGMVXd1kw9N0dK9Ti2jUxYxzACDoV1vgxtQ6cpb2+KUBAdxpSxmniOueipRnZK4U4De1pTO49hBEyFbp1DyiUNmzdFQKIdXW7SNnGPPiEave/9BnoBsyG6Nc4BzbLNBG+KgERbNcW+MIiaNKG8gEnHc40HaUyD4k0RkGifT1E6e7ZZN42q1DT9XccjM6RmHsJofFX1McK4pnPmSNbXIwDRjrNscfjaNj68vCyxmY34HMHYZtZx1pUe3a6oj7O+HgGI0v29isGXbnXPovalIvtzr4S6I5iE5e/Ldo6m2vdEAKLXjShNhSFvH24BF7l2pKbo7ykwHRe65h/HaDyVa7pR6cKs+84IQHRKs3wfT89uoOebhzzp0tIMdRNnzOMYYDo0DWd5YwQgCi04kfkleqLOjFS+Pfrt0Q89At+h+6Edt9vsb4/uHpoP3fFDexSlY7ShNyfSuu/qYVW07EIJB6AyrXY9Nc39Hm+uKfVQW3lQoXJRcm5VexFsj6J0jHaOH0rxweUqD8yk+yQtu1DCLTBboF134VEiPsR3Cz2HT+Gh6MkmOdlvN4IhitIx2pCzjTooKyFxl2jZqIQDMJugXdfxKWihJ3F+Gq/103FG6ZLk5MrdCIYoSsdo6wRUvi2k/W3PkHGXSLyohANQTNSulygf1fTSi5+HByXVCcgWJCdX70YwRFE6RtumYEsPXZtUKeEeMNShdj1HsYGIdmvUrtrn7/lszmo3giaK+hvaXnO7Xr1fslq2ltk8KNWBdr159USXpaisue0EoDFu9PP9CJooSsdol8mmoWdxgXcZLVsp4SWQa2jYG+2aa29BMOdBO6uZbieIybn5eBEs0Xjui2B9SEoewomqt03L1kq43q3LqF3zPtZTZnoPch0ev/aRWKKf70fQRHGooi1xTWFJEvtWu2nZWgk3/toA7TruaraZll7IX4bn5yhjkn6+HyEQ/bX+4uwoHaPNXvI5xavcXY660bKNEi4QKKB2TbtbfUWpRS4GoJgmJ9f6Ef7Vr/ZRfEZbwnOBX5pkolrLtko4AMVE7ZpGa+DJR5AKM8ypApSCzsmVOxH00C3E5ydS8hLn1V6GDfcF95JVwqVdUEDtmhZDBA96gx6V4ouEBlw2dc5YuxfBEM2C9SEpeQxT5jW8PWX3rGWHam5CbpJbAu166unuuanptqGdu6qLc7KLk8qcc2vDXgRLNAvWh6TkW0P34HHSZfesZdNtqrrZl2aVBatdz7xaJNwf0Y78crjESU3OGduwG8ESFfj5Cl+iZ+vTb49+e/RDj8CPGrodrc88/TbRfc3PW6975m79UedoHrooWKPtdTn6oO1hntbZIPxE4oNQRM/Irm2HfLtsexQFa7S9XOiDtodx6q5GqLBBHjzEFyJ6RsbHpzo9WxiiKFij7bRve4bNq7SPYIo4Xd0YMcEGafnBt0CFChTRM/Ial9xNcb0UORqiKFij7aVCH7Q9jFfXatXEBJnGeV85RRE9I4cI6qIkjkRRsEbbayD6oO1hvDpD1AQZqgdEYyglomfkbQ1KqIiuukdR9UPbax/6oO1h3DpN1AQZp2dElYiukfW6jsLTDF0UrNH2Gog+aHsYt04T1UE6uqN52KNGRNdIEsFXEmS2TfcoCtZoJ4z+Rh+0te/DsiUaZ8yopIc3Hg+J0v4souv01/rWryKva6K630Oj0A51uKEP2ujPdjf3ss28qFITVUGuYe8zopWI6Ao5Bqm0XdNVKxD99fsfzs9n8lUUabTZy36iD9rWe9/SRPMy7i6+wnpAFEX0nJ5fIl8Smf9+q5t6FKzR9pqJPmh7GK/OEJUg1zps61LTmepuKKILMr37qePlFK6jfLXNK67R9lKhD9oexqszRCHIgxsGFNEzclsWcN9WB+hzVK24PiRgc3tllfYLGIfpRtQNwm+FHRB1FIjoWYO/xRl3TEPBEn1RwI6p/wRTtLmpl3WOo6xUzof7uqYBWAArFNFzcy70G9G9m/oyzmlqbI+ehlZJ5Eu0PCafXfPt0c/uv7L13x4tj8ln13x79LP7r2z9D+3RrP7yIUG7PFD0eA5/M4K2h7F1iEDbelvLXeudBG+9utv2aFZ/ORzaNonvcwRj4yACbettLXetNwvednW3IcoPc6+J0X+CgbbCH5VgROsNlrfWexO87epuQzSrvxwObUgSTfRB28PYOkSgbb3B4tXC+V9PaPcmeMPqbkOUBZf8lyFoQ5Joog/aHsbWIQJt6+1Zdq33JnjD6m5NVEloMRzaXg70QdvD2DpEoG29Xcu8sUiCN6zu1kSN+ksR0faSoA/aHsbWIQJt6+1ZWY6nvUnwJk0mqCyyUs8SzbpxCKjV4GB7G/qg7WFsHSLQtt6OlfQ+3pUEb1zdrYnioEHbSVKI3EcwNg4i0LbejmXWemfBG1Z3a6JZN97C8bSQBW0nyx9hbBzMgrb1Li39mlgJ3rC62xAV9XeLhnaZJP2FVv7DjyMYGwcRaFvvwjJrvZXgDau7DdGs/nI4tIskVIE+aHsYW4cItK13YTlrvaOqDau7DdGs/ro6cpEjVvy9gP13EdJa763JoU0seMvq7thMSzSrv6WOHN2djz/B2DB/F0HWenOT6SViErxldXdIZ4naBpzK+hI9VXd+h+7ZuvPn9Wj8Xc3T34F9aj/LPybTMsiwwe9AP5VW2e74j8ltW/0PyX2FgtFvYusAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}7.0 & 5.0 & 6.0 & 6.0\\\\0.0 & 3.57 & 6.29 & 5.29\\\\0.0 & 0.0 & -1.04 & 3.08\\\\0.0 & 0.0 & 0.0 & 7.46\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡7.0  5.0    6.0   6.0 ⎤\n",
       "⎢                      ⎥\n",
       "⎢0.0  3.57  6.29   5.29⎥\n",
       "⎢                      ⎥\n",
       "⎢0.0  0.0   -1.04  3.08⎥\n",
       "⎢                      ⎥\n",
       "⎣0.0  0.0    0.0   7.46⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Permutation matrix:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALgAAABkCAMAAAACEQqhAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRIlmzd0i77t8bBwggJIAAAAJcEhZcwAADsQAAA7EAZUrDhsAAASLSURBVHgB7Zxhd6IwFESjIu1Wq3aX//9bNyEtZG5aEg4BtSd+aBnf8Lg+qMpINbuuv+3Nk9wunteYXXdo7O34JNzm7Gj3nQPfPQvzyHkW8PZwOp+CR0E9rvftUsLOMjV6sqxawd/t4dJerkMH6qHwtXC8tV+L9nfCzjJ10MktsqxawM83t8Lp3f10N2p/7/Cz3R/2XQCesLNMPfT1CyxDC/hH/8RyHGCo0drKZvBakbCzTI3uLEMLeNeDv3Vfzy/UaE3whJ1lanRnGToEb7uDW/mta3wPanR2Mpx4ws4yNbqzTB2CX7uTW3vnfxlDjdZOhuAJO8vU6M4ytYL3Ew/AVaN1DD5pv/r9mds9ZQ/BuTeoE+AJO8vU6M4ydQhu/PF/xB/nqNHbyvBQiVaHfWb3hF3APz7cts7DUxw1SKwU8ISdZWp0ZxlawP1z/AEvQKNGa4JzddhZpp5nF3BzcS/5N/tmpe16+kGj6SA/J55nH7qVsCt4e7LvF/s3WZf+tWjUA2q4cNjfuvf92d6VZR+7FbAreEj14MsVfOsdVCdeJ545gXqoZA6qmK1OvNgoMxv9kolr5mLfayEgwjRYpl5mx9qA0Ylr5hJnMui1rt2YqbhJwP075OxAaF17Im4ScGQu2yY82JtO6vmVplUCjszl8yRyDIjQe107wbG1EJwn0tTgZpl6md2tHU6c3UNwZi7UIGGZepmd4Oyu4JroMJOJSNa0x+C6tRCce4Ma4CxTL7MTnN0d+Ev34jeSyGBAEgVAXB1+lqlhtzI8xrm1v+FHKchcDDVas0y9zE5wdA8Plc9PIMYAaFlkA+653QkOGAGfGwgVTXj4OEfwb/MjBS8a2UQkM7tPx00KHm3rce+o4FvvmzrxOvHMCTz1odI2wQfymQ/47rar/RT5F1yvcvc55gM89TH+86FSNuFhN2qMm2XVnPhUBoPGcV7EgAgrsEw9yy7giQwGjaMLiBgQwc8y9Ty7gNtV5WyJARFas0y9qn0S3J8V3ikQ4sahp8B5Yo0Jsky9rn0KnBkMSFimXtc+Da4ZTESi5W3zoynwbfc9xsKNUzvwl9c/w1ryrJKIbFimHpr6BZap59n/veqbLAFHBoPOUV60qX3qUIkiHJAjotnW/gP4txkMuO3lNfMuKCpqV/DpDCYCn5nwFLUreIT2uHdU8K33TZ14nXjmBOqhkjmoYrY68WKjzGz0SyaumUvpK4Q4S24NdZZV68QZ0VCjNcvUsE9f8hOZ2U21gDOioUZvlqlhLxs3CTgTHWqQsEwNu5V6fqWX/NDNbtAC7s8CxwCIGr1Zpoad4Ak7y9AhOE+kqQHCMjXsToYTT9hZpg7BmehQg4RlatgJnrCzTK3gayY8MbhuDY+T8RJ1CM69QY3WLFPDTvCEnWXqEJwX4UQaKP7PZfyPMmrYrQyP8YXdBZyJDjVIWKaGneAJO8vQAu5fQda6Qojg3BoeJ8vQAr5ywjMcKiXiJgUvGtlggsYUjZsUPNrW495RwbfeN3XideKZE/CHypN+o03rviGmad4yH+vdbf032jSN+Q8nTHUiyppjHwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}0.0 & 1.0 & 0.0 & 0.0\\\\0.0 & 0.0 & 0.0 & 1.0\\\\1.0 & 0.0 & 0.0 & 0.0\\\\0.0 & 0.0 & 1.0 & 0.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡0.0  1.0  0.0  0.0⎤\n",
       "⎢                  ⎥\n",
       "⎢0.0  0.0  0.0  1.0⎥\n",
       "⎢                  ⎥\n",
       "⎢1.0  0.0  0.0  0.0⎥\n",
       "⎢                  ⎥\n",
       "⎣0.0  0.0  1.0  0.0⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Have we obtained the correct decomposition? True\n"
     ]
    }
   ],
   "source": [
    "P, L, U = scp.linalg.lu(A)\n",
    "print(\"Lower triangular matrix:\")\n",
    "display_matrix(L.round(2))\n",
    "print(\"Upper triangular matrix:\")\n",
    "display_matrix(U.round(2))\n",
    "print(\"Permutation matrix:\")\n",
    "display_matrix(P)\n",
    "print(\"Have we obtained the correct decomposition?\", np.allclose(P@L@U, A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok... But what is the Permutation matrix?    \n",
    "For more information on this subject have a look at the [wiki link](https://en.wikipedia.org/wiki/Permutation_matrix)\n",
    "\n",
    "We can see that the product L@U is almost equal to A, but the rows have a different order.  For this reason we need to multiply it by P which produces the desired permutation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALgAAABkCAMAAAACEQqhAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRN2JZiK7ze98bHWeL8UAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAWRSURBVHgB7ZzrcqQgEIURHJNVvO3y/s+6dCMgTLTpUlMxJT8yk8yZ9tggwmdXRGWwSXGT1jm/QlRG1ba9buJb9OBWGjBe3cVz9NmvjVedVBqa/7xRutfbZzXVjahU7ClCLoTW9SrcIXlivHeDx4TRPllTTTf488hfZ9CHsxSCkDdTK4SaQpRD8sS4chaD736Gg+h4qHBM92bSsl6dFSUflf3aOPogx+SJcZe82JkjnsLLNP5Y2Ws4Q/d3Qt6a1UnabxyTJ8bx+G0d7Lkx05o4isNH+CYzTshH7L8Y4Zj83XjoStEY6FvRmngq8bDwTvZ9rcK1S8nNVOla2WHu2kH5m3Hd+8hicBdetb7+wofwZrRjaDDeCiFvzAShZ999B+W58cZOj74NLuPbxlE4douekDcGrxU1L1fMQXluvF5NIVRnOsPKX3OU3OAZ9r6HDspz410c4kK4y+e1dXFKPEntnVDyDuW98WORiE7IM+ONi7Z0vptz+63pcEYnKnxMyCXOKiHjy4y+GZ2QZ8btemsxDS/uFrG6160+s2+d1GWmQF4tY9wHIaIT8sz4azFupwCI38Etfw4Tnj/k8jqA8362t5UiuZD2BtfArHKGPDM++Nmqw5tLo+36ccu3nQqVlgoniSK50Eq5aCfIM+NZTn/wr4/x7+6cJ+NPxgsz8AyVwkSdJnsyfloqCwP9kowPSo3jam1CIBsmEGrqWsc9p11r7eImgk4lGW9gZRU3BiThYQIhXAZLv+ekohN0KjEucaU3B4xAIRseEHrh1mfA9TIMZCI6QacS4zPuCmVgoASyERlXIeQad4WN31tTQIigU4nxCfco0TiBbHLjhLw1gDMi/iDkOLts06nEOGpF5zeR1D6cCYSENKaPvsno4CZu3XP5u/E28B8C2TCBkLWhjZkWqgIkCQfDLrWJp/kmfzfutlVwugSyAYltpUBIiJdsOzNH8IWzzJ7xPTr1Zly5vTtYynsH/vbeioFQD9dybTzJKoi+R6fA+If5CH7q6PtsIAQ4AJ7b+IncXZybuMlq9+jU3/WjFNub4HvwHJsgPDwg5PtdepJFRIceX8+3mTwdKhXmOzxlIJCN65xiIOQyLkafFiI6Pg2EPlpaJk+MD519cKWkvUUUIRsmEGpxRuktbC+Kbrt/j04lxif38ArubUXIhgmEKvtMzy2yiqKLXTqVGPe9cofXx/h399KT8SfjhRm49VBpbBHE7dpgb2T3r1e5U9pvPcbjUGESHught6B0fbVPeFDDkO/zozTjTMIDVtaLfaLkhynHVfMmP0qN8wgPGKnnuBEnCA9TTvCj1Ph6x2GPQxAeqxj6KRo/V07wo13jNLLRYmX8XDnBjzLjrJIfi/+GlXF6286T7/Oj1Dir5MduwezdK2acJDxM+T4/So3D9cMgPICi1sZxHtgmPEz5Pj/6wngx4XnBI4BonBoqTDnBj8D4x+cfyLStakN2Hcn+PrJpkP1F40SFEFMuCH7073O1yOIRnpeEZmbpq/syZOOSEX4y5RQ/SoYKk/Cgp7cb0GZBEVNO8KPEOJPwoBMDN6BCwmOV5XKCHyXG2SU/Qo/GAAwsIzw8+T4/So1jDu/x4zH+3f30ZPzJeGEGnqFSmKjTZE/GT0tlYaBfknEmEGIWFDHl1wEhZkERU+7q6i8BQsyCIqb8QiDELChiyi8EQsyCIqb8UiAEU1d5QRFOdOXyC4EQOGEUFHHluwVFX8zj5SU/4IRZUMSQXweEwDezoIghLwBC4YkEDwiBb05BEVNOAKFkqPCAkDXCKijiyS8FQryCouVxUWn90YVAiFtQxKs/uhAIMQuKmHLxACGYRX5OS2aVn2OLdvIYp3N0ruLmGcdCyeyx8rkJOjNa+I82dlcKzRd2n3mMS2Lhf7Spa/EfMYpzvimkHiQAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}7.0 & 5.0 & 6.0 & 6.0\\\\2.0 & 5.0 & 8.0 & 7.0\\\\5.0 & 4.0 & 4.0 & 8.0\\\\5.0 & 2.0 & 2.0 & 8.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡7.0  5.0  6.0  6.0⎤\n",
       "⎢                  ⎥\n",
       "⎢2.0  5.0  8.0  7.0⎥\n",
       "⎢                  ⎥\n",
       "⎢5.0  4.0  4.0  8.0⎥\n",
       "⎢                  ⎥\n",
       "⎣5.0  2.0  2.0  8.0⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display_matrix(L@U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a deeper look.  \n",
    "The coordinates of the ones in the matrix P indicate the relation between the rows of A and the rows of L@U:  \n",
    "- (0,1) The first row of A corresponds to the second row of L@U; \n",
    "- (1,3) the second row of A correspond to the fourth of L@U; \n",
    "- (2,0) the third row of A corresponds to the first of L@U;\n",
    "- (3,2) the fourth of A corresponds to the third of L@U.  \n",
    "\n",
    "This is clear by looking at the columns of the matrix ```np.where(P==1)```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHQAAAAyBAMAAAB/teJlAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhCZZs3dIrurRHbLQ9+lAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACK0lEQVRIDe2Wvy8EQRTHv27dHcdxgkJEyEUoFHca9UlQLxGJRGKra1yh00q0EkKikrgoJDQUOoWT+BGV8xfoFAriV1CwZudmxps5ezglW83Me599b99uPlm0u3co40q4t2gcGCyDxOhAL5rKAT2mn6Nj3Z2lbrAxtUvDkXhPUqFzGM/SoNVMd8EcEnlysAHriW151XAMgVkSm8g8kh1qbVTEyMFNEi9sy9GaHCwtuU7bVSyi9oGg5/ngs0Trc4jekxh0tOpRR0EannQQffNHWaRaawPjHeyMN7zdgUrtu9CrsrQLL1dd6ys2WxdQ5yt0VWF8EWiT6JcNhxwdRQMry6uyMVklxgRcauQWsJ+UL2cHYW0OxrNaDli6ulz7A2WfRGhWRdjCQIeAKxJuBdayoiqWMZwiMQONLMUziyQ8A8t7lfxZMXJ9SEIIHTyfkH2V67oUjaan8goled9fFqp+P59k/qNkGP7LPz8m4Wg1IcPwkfReisXEmHptlccWwtHqyDD8BAKvEg0en2mocLRETcMfAQ0SBW40VDhaoqbhW4CzlGpYR4WjJWoafsEuhXJHS7RImEwwrEsxJqNqwdESLTI8op7wP0fBHa1Q0/Cozvmj3NESLW447YU+rSocLdFiwzu+qHC0RGtMw/chkvWpKhwtUdPwlQ7q/FDhaImaht+Md52yWOFZM+3z0yoREI5WJ4bhF1xX/UuonB8txIR/xIjkv4X+4qe2/F/pd2FfrOCPai47AAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}0 & 1 & 2 & 3\\\\1 & 3 & 0 & 2\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡0  1  2  3⎤\n",
       "⎢          ⎥\n",
       "⎣1  3  0  2⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "After the permutation the resulting matrix is equal to A:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALgAAABkCAMAAAACEQqhAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRCK7ie/dzWZ8bFsxDvgAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAXgSURBVHgB7VztlqsoEEQxzl2/vcv7P+vSjYCNkabPxjnjHPmROEnZlk2DpihHVQZbrW7SZsdXqcroxrbXTXirFdjWBohXd+Ecea6EeK/1MOzOotXjOu7+jvvh1tK0qtKxpxh42zSj7mKMPLyaaz1C8ztQOCHeQpmPJsZeLKl27v2u6fsExRYCK8XANexfx/PMw1dXySYMPQonxOsWQk8TvEJbcWtc3F/H12Wsm91ZMfDXChH6EI2Baxc58E7ghPg0Q+g6lPyAe70Mng98lbQQ1H3OwMcBYC0eA7Y4OGBULNQETogvSDESd93Umdi7GCy8JMQZeGcGm4ER8w4hGDgepWvCwRI4Ie5As89wa7AqOxP3DmFwo17XRoexy8ONWSNvFg6HwE7CY6XwI/EuDLfebVXhA4ywe4EM9mEss3A77s0Syo6H77vHHgdngUDmSHwOBdC7jAfsjvJuc/BFy8JfdTebyc9ZLNyOBztX+5bCD8Q1lgfC097xMei7NtvEwsFXyEhjpi3nHBzAYQaCk6B1mxJvIm8/fF5ng7PGuHHed8PnFD7hCVYhGgO3xOdY4gcyCfEX8O793DzgjqsfrPYr0iYkrsPXebjv99pPK3m4PVLrTm07ZgKnxCvMd7iouDlf7zqMEHedM4evGbjLuBp8Whg43vztjpfACfF+trcGurajrTVIZ4ZL/hQmvF0Y2OyB+Qp8iuAdziirnVuL4Eq9XFm/hxPii7s9gGnCzS3taO8fz3gD87HWONaK4JW9a3I3WUVw1U/ble8dnBCHNN6lPcS/u6eejD8ZL8zAUyqFifoY7Mn4x1JZGOiXZFyo8Aj1IyE8rx/RjMsUHqF+JIQrvGs+1Y8ocZnCI9SPhHBGP6LEw+9kN0ISDSYdNkL9SAhn9KMs8USDSYkL9SMhnNGPEuIyhQfPpFw/EsLrrH5EiUsVHktFoB8BcQk8qx9R4piTcoUH4O5nFWylkg18dmgCeF4/ekO8WOEBVkL9SABn9CMg/mW+XGaECo/dSaIfCeGMfvR3v5QiU3gsEZF+JINz+hEpFaHCY9d/oKdK9SMhnNGPCHGhwiPVj2RyE6MfEeJChUeoHwnhKq8fUeLQ8zdpD/Hv7qgn40/GCzNw61JprQnidq23Kxu/wK9yo7zfusZDqTCmnEN/CPUjuzY/NtHHoajl5xA9DycZZ0w5h9Ay/Ui1i13H362aUsvPIXoeTogzppxDaJl+pAa4fXcLxBDKrbie+4/ycELc2atiZzKCkJLpR513K2wJYKIzcEIcI56bcrYDxreEOKMfDcHr5SL8P/iReHQspE6LyHjbkjmEzFKNTbTvcdEZ+IF49B4dTDkH4iL9yC7gg2/CL3Nz0Tl4Stz/uAaSRQqPKtWPWoP+EO2NNkx0Dp4Sz5lyDhnHD4r1I4MeqNV7uNhSycOB+NeffwKnnCkngLYNoX7kjC2r8UYbNzhPDUUM/N8/5CYra8pJiQv1oxpnlZDxbUY/9R8x8KRUrCN7Ry8x5ey+wU2hflRtNe7DMNEZeEI8a8rxR/TvQv1I1fYC18KsUuYQysMT4llTjicc3oUOITVq7fxG7yw/IWrYyMIT4mGnH7/xEP/uLnoy/mS8MANPqRQm6mOwJ+MfS2VhoF+ScbHCY++Jd7eTrMIjgl/nEMJO3f/yYBQewAvgmJFLHEJApJmiKsApPDL4hQ4hS6Rfl0icUXiE8AsdQpbJqHbEGYVHCL/UIbT2O+Lcz3arForgFzqEWrucETOePtRlM0ybEJ5/wuzNPF6q8EChUOI4D5w/YSaEX+cQesHzKjHjXKkI4QUOobAiIVN4WlSlI/HtAc0zhUcIV4xDiJSKTOF51dDMVPvHPZOHumiBKyHci5hnT5gR4kKFB5kdLkC7tZKEOvxZDr/QIYS8DFyAyhQewJfDL3QIWR7jYAx0ZpnCI4M/DiHo6Z/TyOD8ObR4Jg9xPkefRdw84/jkbLLc+tkEfTJa+I829lcpNP+k/yePcUks/I82TaP+AwKJdYmcSJNDAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}2.0 & 5.0 & 8.0 & 7.0\\\\5.0 & 2.0 & 2.0 & 8.0\\\\7.0 & 5.0 & 6.0 & 6.0\\\\5.0 & 4.0 & 4.0 & 8.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡2.0  5.0  8.0  7.0⎤\n",
       "⎢                  ⎥\n",
       "⎢5.0  2.0  2.0  8.0⎥\n",
       "⎢                  ⎥\n",
       "⎢7.0  5.0  6.0  6.0⎥\n",
       "⎢                  ⎥\n",
       "⎣5.0  4.0  4.0  8.0⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display_matrix(np.where(P==1))\n",
    "print(\"After the permutation the resulting matrix is equal to A:\")\n",
    "display_matrix( (L@U)[[1,3,0,2]] )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous code is useful to obtain the decomposition.  \n",
    "Now we can solve the problem \n",
    "\n",
    "$$ PLU x = b $$\n",
    "\n",
    "using both forward (for L) and backward (for U) substitution.  \n",
    "Since permutation matrices are orthogonal matrices (i.e., $PP^T = I$) the inverse matrix exists and can be written as \n",
    "$$ P^{-1} = P^T $$\n",
    "and there is no need to perform a matrix inversion.\n",
    "\n",
    "The method **solve_triangular** performs the forward/backward substitution (it is a wrapper of the LAPACK function *trtrs*).  \n",
    "We can write:\n",
    "\n",
    "$$ L y = P^T b $$\n",
    "and\n",
    "$$ U x = y $$\n",
    "\n",
    "and solve it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACoAAABkCAMAAADHT1J1AAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRO/NZt2JIrt8bAICnK4AAAAJcEhZcwAADsQAAA7EAZUrDhsAAAJYSURBVFgJ7ZjRlqMgDIYjIDMrCrib93/WJUEEPRLaPXsxF+VCW/mbxiR8RGFCHgq6w2cFwITapDF3lbDRvEKSTn1VndkkqdVhC6eZm3R2tpqBNflkfTyutFKrtMJGujkShfVBmi6ZVrpzTOZyqbV6lyJLFzxCI0gtavrrBQ2dAARpxECKKZ8GUrb6ivQNByDf1vzCbcG+57s5Qi3cFuQUaDkFFnneU2JdKYKLVa0crmpLf+s5+jak0ivKe1zJud64WO2J8vWP9GdFIGq972fOU6l0kWGpRgIuZ377yFBcw45BQXIBGc6TQJ2wE5CxMkeqVEZGsuoLeAYrlmDCoEi/GSCjLBVyOmYQ9ZABmolCSnhy4Au/eC4dTFXCHRm/L4CfSRkLpiVkTGzTFKmAjOhDCFqlRAyRseaNjHL2QUbJXnt+c3Fb02ysrZ3r55i2uv/RD1ytvunr6UBM5bI3LUwfGZFqcEbaYfLoI8PwFolMA9IKyFhcYtBRqiQVkEHTrQMDZFhfMPC4Ys8IQAxen5kbIiOQwzxGyABwpSt7RMb3r2yGV/Veuoo7Mv5813JBZmCVCshw3IB4fAEZG6V0pg5siAxYdNpiOAAfZOQquB7/dRlerdy//TyreUPKfvaRwfOeKzFL+8igeZOLlqUCMtJ83NZqVUZGgEYqImOLjfRpxZ7IsGmtVKsiMohBrZQ38scuY6bGrUoFByyDrUrvyGjKZVY00Knj+WrwYJKQVeMqdBmcIiTpGBmpz9wR6cnkgwyO2+2Qs/XyWwdL7xSMKZvQzRh95bcOxsBfLMkuJFdQyOYAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}1.0\\\\2.0\\\\3.0\\\\4.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡1.0⎤\n",
       "⎢   ⎥\n",
       "⎢2.0⎥\n",
       "⎢   ⎥\n",
       "⎢3.0⎥\n",
       "⎢   ⎥\n",
       "⎣4.0⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "y = solve_triangular(L,P.T @ b, lower=True)\n",
    "x_solution = solve_triangular(U,y)  # by default it considers upper triangular matricies as input\n",
    "display_matrix(x_solution.round(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Another LU solver\n",
    "A faster method to solve the linear equation using the LU decomposition is to use the **lu_factor** method (it is a wrapper of the LAPACK function getrf).   \n",
    "It returns two matrices\n",
    "- The first, LU,  is a matrix from which it is possible to derive easily the L and U matrices\n",
    "- The second, piv, corrsponds to the pivot indices of the permutation matrix P:  \n",
    "   In the matrix P, changing the row i with the row piv[i] (for all $0\\leq i\\leq 3$), produces the identity    matrix. The same permutation transforms A into L@U.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Construction of the matrix L:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAABkCAMAAACcqw0BAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRO/NZt2JIrt8bAICnK4AAAAJcEhZcwAADsQAAA7EAZUrDhsAAAgRSURBVHgB7VzdorI4DERAzreKoru8/7Nuf0kyTQro1ceBG5qQmXYoICRgc5nD0ja/ZHlGvU1zmbveLddfort5eLXt7IVffotm0vlYFz5002PatWsQgTZ1v9ZCJNo1PMZKWxV+vQ+McXTnwPB8Mc9aExFor+FpOyLRpsiyhbHSLoUPbdfOTPjj7jmnsWS2PIhA28KVfkSiXSLIg7Fgl8IdtOfC3+Fqf+UuYldbiEBbBalORKKtgpITY8FeFz4H4bd5+0UfEWjXhiu3IRJtGS0tjAV7Vfgwd57wNveS17YQgbaNxC2IRBvjuY2xaK8Kf82T57vEFae22ohA28KVfkSiXSLIg7FobxAeZnyXcIl4xWNmB0MePiLRznHaGmPRXhWOh4jWifQhAm0ZXbMQifY32FXhTbwoXHdf3AixnyErQiTaOU5bYyzY68Lfb0/72PNzBoj9DFkIItHOcdoaY8FeFx5/+LvdNzCE2M+QhSAS7RynrTEW7JrwYQ5qn/6W9b7nZn1BfMyQlXzDVMcqwrv2Po/tw/X9DPcuw+Qe4vbobgjxKUMW/g1THasIz50ee30KP/b8lurOGS/3ybE954wfe35LdeeMl/vk2J5zxo89v6U6OeMy5940r657v8NtutsytXty67wrpOXb6m1Eol1DY6y0pXCZc28G/5AyzTf/dOYeWl73zcq/LUlkQTCgBu0cp60xVtpCeHxkpdpBG+oKd1dRaENVoQ3JNK0T4fu+JJHpcEBo5zhtjbFgC+GQc2/uT8/YuppifDKfXKFt2/JlSSJ3ggNCO8dpa4wFWwiPaSmqHYwh3+SED3PIJvWby6pCONJqw9R9iERbR0UvxoLNhRtJzKeTfw+5mHC61zpbtnHhBu0SazcQibaNdFclKISgzYVjzj3y3nwpIR7k4+zzMlsWLlyn3cKCSLRrHBiLthQuKwGRN2aPRldBuk6bc8xSuEZbG3TehkUAtHOctsZYtLlwPBoCX5eu5O79icv00Tmu0mpjLXyIRLsAMAfGou2F/8w/ERHPf6oEOG/Pf8FE3Zx1Ujb5jH9QksiEOCC0c5y2xliw/+WvgkDO3dFdve5Xvm8ZQ9ZV6wV9QnhJi+GWjUi0LZz3YyzY/FBv4m88VQKaS5jv/uW2OPGXOe+BWodhmxBe0K7CcwAi0c5x2hpjwRbCG8jBv57TNHWtu43px6EZxq3X9OWdir+moAA5+DG+Cefv36au69w9+7bl+5JE7gcGxEoVOcJe17Fyxm2Ww205hR9uSlcEnTO+soMOt/mc8cNN6YogP+NDz15dXYk/zOaXe9w831c/zHRuEHJe3MJOkjn35vJsXR3BLRv2YAoBBgak4gRz7m3a9CUTxkpbzrjMubu3+uISHsNljaDsJ3qAwTkTjooTFlT3y25Leh3lvRgrbSE8PrJSQaGLz99ON9YIrO6QgXBUnLCwip/gcSPSK5DFhbFgC+GQc2/iIZ4+xBG5hYUeGsjgNiccFScAsmKKbhV6E46xYAvhMS1FBYVAektv6IsRWN0pDAm3FCcsqOEX3Sr0BqxJmT4SA1guHBORkTPUUFxTjMDoTmMQOF+c2LdwuEZvsWEs2lw45twD55TzTXwEVm8aA8eF4oQF1v0crtHrKJchhS8r0JbCQ25RfEowLHVCPgK7t5KB42JxwkKrfg7HooAKSE6MRZsLx6PBU/TL69p8BFaPKgMd3bk4YcE1P+9Wo9cw3oexaHvhP3/+ifB4/ouCwjOf4pvOca12QCMXxQlrwOgnuNuiDBDjFxtjwf7vD3tIgZy74xhiuGcTI1jooVEyEE4WJwBomqJbhd4EYizY/FAvCwruwS2ctJ5djMDqLt4lsJIE4ZbihIXV/aJbhV5HOS/Ggi2EY0HB3W/uFF4y5B1GxQlzqOqGJPyD0gRUR2hooR8pHHPw7n2f9GUl1QjU4S1OZFhwrDixBK83FvgHX0vgUMj2/Urh6yM5TMQp/DBTuVHIOeMbd9Rhws4ZP8xUbhRyzvjGHXWYsHPGDzOVG4XIGZc5908KCqlbIAre+Hjmm9RK4ZVVjWma+h3/xSSZpHCZc/+koJAkAFHwUk6DWhXFq0zD6F7DEk/AyFYrRgjh8ZH1m4JC6hqJvLu/52QOtXCkpV1hevtUQUwvlLiyBgJMQjjk3D8oKKQRIJFzvx5jEk4tbbzgs5lu6+9ZiiQGMAnhMc9EOfgwiF0FhTRshWhqsnBqgUjNtJne4TMZDbP4hHBg4sIxERkJlgOUsqULs9FQiB6vLJxaBpq7K0zzeJn6+tuWXDgyceGYcw8j2FVQSGMuiQZXh4ozTi0u0GrbTC4X5UsdOUOkEnDhyCSFl+WAfQWF1D0m792rsG5LFE4tdazgtJmGORyBnfjPPUBL4VIcF45Hg6fZVlAYxueyjJcimd9c/deKQTi1YJDJ3M40h2/DHv5zQGvhwlEcF67l6+knl7NYPWV/vJAslYkh1Ju9cGrl0JW1ydQ8Q43nUfs+SAwZmIRwyLm7Me0tKCQZQHRt/TLf255aK4LzZpMpffy4ecbxiwUhPP7G85uhvQWFNN6SyG1YbmBYK+sz1zbTJZ3jJpRqGSEEmIRwyrmn/P3+gkIaBibzvTt+rShbKbyyspladwINW67qajFCCqeceyro7i4oJAkFUTO95zn8jxJrVfQumypM7quJ2l8x1YsRUvjS3fEbp/Djz7FUeM643B/Ht375jIdXVjd/OPq3Hw7LP+oP/h/m+75yt/+3S5XjD/+o3/fN/7JjgV+jf1/TAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}1.0 & 0.0 & 0.0 & 0.0\\\\0.29 & 1.0 & 0.0 & 0.0\\\\0.71 & 0.12 & 1.0 & 0.0\\\\0.71 & -0.44 & -0.46 & 1.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡1.0    0.0    0.0   0.0⎤\n",
       "⎢                       ⎥\n",
       "⎢0.29   1.0    0.0   0.0⎥\n",
       "⎢                       ⎥\n",
       "⎢0.71  0.12    1.0   0.0⎥\n",
       "⎢                       ⎥\n",
       "⎣0.71  -0.44  -0.46  1.0⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Construction of the matrix U:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOgAAABkCAMAAAC7BYzpAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRN2JZiK7ze98bHWeL8UAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAjQSURBVHgB7Vxd26I4DEVAZlY+hF3+/2/dtKFpchoUZx4v5JULbUpO0kNLgUO1uqxxq6uTbj3zq6rL2rS0XU/KsxoDu3oNRC9n5Zh5jZropa+bIWxpf9cM4/DkKMxtV12aPByOYFL89D0MrUrzlgiG6MiDeZWzdab2d/2U2uN+LwEkh6aqjmBsoG6+VVUzS+VbIhiiDVMSnuMSkg+5CdIWVZiHulWH4hBGwal4b8LHPVW+J4Ihyv2SB9E9Ur6uXWqD9y2HhXcewpg4t1UdJ9rzngiGaMx/a6UZPIZvaz4BZVcuANFDmIym0j2Om1z1ngglURlCVbeGMVXd1kw9N0dK9Ti2jUxYxzACDoV1vgxtQ6cpb2+KUBAdxpSxmniOueipRnZK4U4De1pTO49hBEyFbp1DyiUNmzdFQKIdXW7SNnGPPiEave/9BnoBsyG6Nc4BzbLNBG+KgERbNcW+MIiaNKG8gEnHc40HaUyD4k0RkGifT1E6e7ZZN42q1DT9XccjM6RmHsJofFX1McK4pnPmSNbXIwDRjrNscfjaNj68vCyxmY34HMHYZtZx1pUe3a6oj7O+HgGI0v29isGXbnXPovalIvtzr4S6I5iE5e/Ldo6m2vdEAKLXjShNhSFvH24BF7l2pKbo7ykwHRe65h/HaDyVa7pR6cKs+84IQHRKs3wfT89uoOebhzzp0tIMdRNnzOMYYDo0DWd5YwQgCi04kfkleqLOjFS+Pfrt0Q89At+h+6Edt9vsb4/uHpoP3fFDexSlY7ShNyfSuu/qYVW07EIJB6AyrXY9Nc39Hm+uKfVQW3lQoXJRcm5VexFsj6J0jHaOH0rxweUqD8yk+yQtu1DCLTBboF134VEiPsR3Cz2HT+Gh6MkmOdlvN4IhitIx2pCzjTooKyFxl2jZqIQDMJugXdfxKWihJ3F+Gq/103FG6ZLk5MrdCIYoSsdo6wRUvi2k/W3PkHGXSLyohANQTNSulygf1fTSi5+HByXVCcgWJCdX70YwRFE6RtumYEsPXZtUKeEeMNShdj1HsYGIdmvUrtrn7/lszmo3giaK+hvaXnO7Xr1fslq2ltk8KNWBdr159USXpaisue0EoDFu9PP9CJooSsdol8mmoWdxgXcZLVsp4SWQa2jYG+2aa29BMOdBO6uZbieIybn5eBEs0Xjui2B9SEoewomqt03L1kq43q3LqF3zPtZTZnoPch0ev/aRWKKf70fQRHGooi1xTWFJEvtWu2nZWgk3/toA7TruaraZll7IX4bn5yhjkn6+HyEQ/bX+4uwoHaPNXvI5xavcXY660bKNEi4QKKB2TbtbfUWpRS4GoJgmJ9f6Ef7Vr/ZRfEZbwnOBX5pkolrLtko4AMVE7ZpGa+DJR5AKM8ypApSCzsmVOxH00C3E5ydS8hLn1V6GDfcF95JVwqVdUEDtmhZDBA96gx6V4ouEBlw2dc5YuxfBEM2C9SEpeQxT5jW8PWX3rGWHam5CbpJbAu166unuuanptqGdu6qLc7KLk8qcc2vDXgRLNAvWh6TkW0P34HHSZfesZdNtqrrZl2aVBatdz7xaJNwf0Y78crjESU3OGduwG8ESFfj5Cl+iZ+vTb49+e/RDj8CPGrodrc88/TbRfc3PW6975m79UedoHrooWKPtdTn6oO1hntbZIPxE4oNQRM/Irm2HfLtsexQFa7S9XOiDtodx6q5GqLBBHjzEFyJ6RsbHpzo9WxiiKFij7bRve4bNq7SPYIo4Xd0YMcEGafnBt0CFChTRM/Ial9xNcb0UORqiKFij7aVCH7Q9jFfXatXEBJnGeV85RRE9I4cI6qIkjkRRsEbbayD6oO1hvDpD1AQZqgdEYyglomfkbQ1KqIiuukdR9UPbax/6oO1h3DpN1AQZp2dElYiukfW6jsLTDF0UrNH2Gog+aHsYt04T1UE6uqN52KNGRNdIEsFXEmS2TfcoCtZoJ4z+Rh+0te/DsiUaZ8yopIc3Hg+J0v4souv01/rWryKva6K630Oj0A51uKEP2ujPdjf3ss28qFITVUGuYe8zopWI6Ao5Bqm0XdNVKxD99fsfzs9n8lUUabTZy36iD9rWe9/SRPMy7i6+wnpAFEX0nJ5fIl8Smf9+q5t6FKzR9pqJPmh7GK/OEJUg1zps61LTmepuKKILMr37qePlFK6jfLXNK67R9lKhD9oexqszRCHIgxsGFNEzclsWcN9WB+hzVK24PiRgc3tllfYLGIfpRtQNwm+FHRB1FIjoWYO/xRl3TEPBEn1RwI6p/wRTtLmpl3WOo6xUzof7uqYBWAArFNFzcy70G9G9m/oyzmlqbI+ehlZJ5Eu0PCafXfPt0c/uv7L13x4tj8ln13x79LP7r2z9D+3RrP7yIUG7PFD0eA5/M4K2h7F1iEDbelvLXeudBG+9utv2aFZ/ORzaNonvcwRj4yACbettLXetNwvednW3IcoPc6+J0X+CgbbCH5VgROsNlrfWexO87epuQzSrvxwObUgSTfRB28PYOkSgbb3B4tXC+V9PaPcmeMPqbkOUBZf8lyFoQ5Joog/aHsbWIQJt6+1Zdq33JnjD6m5NVEloMRzaXg70QdvD2DpEoG29Xcu8sUiCN6zu1kSN+ksR0faSoA/aHsbWIQJt6+1ZWY6nvUnwJk0mqCyyUs8SzbpxCKjV4GB7G/qg7WFsHSLQtt6OlfQ+3pUEb1zdrYnioEHbSVKI3EcwNg4i0LbejmXWemfBG1Z3a6JZN97C8bSQBW0nyx9hbBzMgrb1Li39mlgJ3rC62xAV9XeLhnaZJP2FVv7DjyMYGwcRaFvvwjJrvZXgDau7DdGs/nI4tIskVIE+aHsYW4cItK13YTlrvaOqDau7DdGs/ro6cpEjVvy9gP13EdJa763JoU0seMvq7thMSzSrv6WOHN2djz/B2DB/F0HWenOT6SViErxldXdIZ4naBpzK+hI9VXd+h+7ZuvPn9Wj8Xc3T34F9aj/LPybTMsiwwe9AP5VW2e74j8ltW/0PyX2FgtFvYusAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}7.0 & 5.0 & 6.0 & 6.0\\\\0.0 & 3.57 & 6.29 & 5.29\\\\0.0 & 0.0 & -1.04 & 3.08\\\\0.0 & 0.0 & 0.0 & 7.46\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡7.0  5.0    6.0   6.0 ⎤\n",
       "⎢                      ⎥\n",
       "⎢0.0  3.57  6.29   5.29⎥\n",
       "⎢                      ⎥\n",
       "⎢0.0  0.0   -1.04  3.08⎥\n",
       "⎢                      ⎥\n",
       "⎣0.0  0.0    0.0   7.46⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "piv vector:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABoAAABkCAMAAACvvNBNAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRCK7ie/dzWZ8bFsxDvgAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAFBSURBVEgN7ZbReoMgDIUjIHYioFve/1lHEqFG2b7eKxet8nMShXAEBuRmoLUsPQADWlfa2AhsdG+Q0PDufV9tCgVrl6WOU8hTuoizKBUynjpT6qCUqdPsqZVqRZJ1kQTKPABAqRjNGGXIFeU6KxdkrWiuAV0j51wjkRBYqAMOrHEdFHKM0Rp+8VOuVZaph+qzyb/Opdjt0W1KdMKprfyhRL/VhoWnRP9zNu/Ym9o87hfBfeiHSvjx1gvFAZbqwUoVyDZG3CSuQg5JgT1zmFOxXY9rR8Vd/YCEfN6d8uxsIWZbp0Y9BgeMlJDaFRWfFx2h6fUl48QjF34FgJ/XYXqRbbwiFTAtpM7YcdGNpmjEsk7UlApmW75S+wOekAzff7XqQYdFgduUaDv26BJVtfGUKH9T/jxVejpDOle3HG1I6YFfTcYj3Y81t/kAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}2\\\\2\\\\3\\\\3\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡2⎤\n",
       "⎢ ⎥\n",
       "⎢2⎥\n",
       "⎢ ⎥\n",
       "⎢3⎥\n",
       "⎢ ⎥\n",
       "⎣3⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "LU, piv = scp.linalg.lu_factor(A)\n",
    "print(\"Construction of the matrix L:\")\n",
    "display_matrix((np.tril(LU, k=-1) + np.eye(4)).round(2) ) \n",
    "print(\"Construction of the matrix U:\")\n",
    "display_matrix(np.triu(LU).round(2))\n",
    "print(\"piv vector:\")\n",
    "display_matrix(piv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Product of L@U:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALgAAABkCAMAAAACEQqhAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRN2JZiK7ze98bHWeL8UAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAWRSURBVHgB7ZzrcqQgEIURHJNVvO3y/s+6dCMgTLTpUlMxJT8yk8yZ9tggwmdXRGWwSXGT1jm/QlRG1ba9buJb9OBWGjBe3cVz9NmvjVedVBqa/7xRutfbZzXVjahU7ClCLoTW9SrcIXlivHeDx4TRPllTTTf488hfZ9CHsxSCkDdTK4SaQpRD8sS4chaD736Gg+h4qHBM92bSsl6dFSUflf3aOPogx+SJcZe82JkjnsLLNP5Y2Ws4Q/d3Qt6a1UnabxyTJ8bx+G0d7Lkx05o4isNH+CYzTshH7L8Y4Zj83XjoStEY6FvRmngq8bDwTvZ9rcK1S8nNVOla2WHu2kH5m3Hd+8hicBdetb7+wofwZrRjaDDeCiFvzAShZ999B+W58cZOj74NLuPbxlE4douekDcGrxU1L1fMQXluvF5NIVRnOsPKX3OU3OAZ9r6HDspz410c4kK4y+e1dXFKPEntnVDyDuW98WORiE7IM+ONi7Z0vptz+63pcEYnKnxMyCXOKiHjy4y+GZ2QZ8btemsxDS/uFrG6160+s2+d1GWmQF4tY9wHIaIT8sz4azFupwCI38Etfw4Tnj/k8jqA8362t5UiuZD2BtfArHKGPDM++Nmqw5tLo+36ccu3nQqVlgoniSK50Eq5aCfIM+NZTn/wr4/x7+6cJ+NPxgsz8AyVwkSdJnsyfloqCwP9kowPSo3jam1CIBsmEGrqWsc9p11r7eImgk4lGW9gZRU3BiThYQIhXAZLv+ekohN0KjEucaU3B4xAIRseEHrh1mfA9TIMZCI6QacS4zPuCmVgoASyERlXIeQad4WN31tTQIigU4nxCfco0TiBbHLjhLw1gDMi/iDkOLts06nEOGpF5zeR1D6cCYSENKaPvsno4CZu3XP5u/E28B8C2TCBkLWhjZkWqgIkCQfDLrWJp/kmfzfutlVwugSyAYltpUBIiJdsOzNH8IWzzJ7xPTr1Zly5vTtYynsH/vbeioFQD9dybTzJKoi+R6fA+If5CH7q6PtsIAQ4AJ7b+IncXZybuMlq9+jU3/WjFNub4HvwHJsgPDwg5PtdepJFRIceX8+3mTwdKhXmOzxlIJCN65xiIOQyLkafFiI6Pg2EPlpaJk+MD519cKWkvUUUIRsmEGpxRuktbC+Kbrt/j04lxif38ArubUXIhgmEKvtMzy2yiqKLXTqVGPe9cofXx/h399KT8SfjhRm49VBpbBHE7dpgb2T3r1e5U9pvPcbjUGESHught6B0fbVPeFDDkO/zozTjTMIDVtaLfaLkhynHVfMmP0qN8wgPGKnnuBEnCA9TTvCj1Ph6x2GPQxAeqxj6KRo/V07wo13jNLLRYmX8XDnBjzLjrJIfi/+GlXF6286T7/Oj1Dir5MduwezdK2acJDxM+T4/So3D9cMgPICi1sZxHtgmPEz5Pj/6wngx4XnBI4BonBoqTDnBj8D4x+cfyLStakN2Hcn+PrJpkP1F40SFEFMuCH7073O1yOIRnpeEZmbpq/syZOOSEX4y5RQ/SoYKk/Cgp7cb0GZBEVNO8KPEOJPwoBMDN6BCwmOV5XKCHyXG2SU/Qo/GAAwsIzw8+T4/So1jDu/x4zH+3f30ZPzJeGEGnqFSmKjTZE/GT0tlYaBfknEmEGIWFDHl1wEhZkERU+7q6i8BQsyCIqb8QiDELChiyi8EQsyCIqb8UiAEU1d5QRFOdOXyC4EQOGEUFHHluwVFX8zj5SU/4IRZUMSQXweEwDezoIghLwBC4YkEDwiBb05BEVNOAKFkqPCAkDXCKijiyS8FQryCouVxUWn90YVAiFtQxKs/uhAIMQuKmHLxACGYRX5OS2aVn2OLdvIYp3N0ruLmGcdCyeyx8rkJOjNa+I82dlcKzRd2n3mMS2Lhf7Spa/EfMYpzvimkHiQAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}7.0 & 5.0 & 6.0 & 6.0\\\\2.0 & 5.0 & 8.0 & 7.0\\\\5.0 & 4.0 & 4.0 & 8.0\\\\5.0 & 2.0 & 2.0 & 8.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡7.0  5.0  6.0  6.0⎤\n",
       "⎢                  ⎥\n",
       "⎢2.0  5.0  8.0  7.0⎥\n",
       "⎢                  ⎥\n",
       "⎢5.0  4.0  4.0  8.0⎥\n",
       "⎢                  ⎥\n",
       "⎣5.0  2.0  2.0  8.0⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The matrix A after the permutations indicated by piv:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHgAAABkCAMAAAB6rgQAAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRN2JZiK7ze98bHWeL8UAAAAJcEhZcwAADsQAAA7EAZUrDhsAAARBSURBVGgF7ZrrcqMwDIWFTWg33Hf9/s+6lkgTS1jSNJmwO1P8oyWc4GMLYz6phSZRC3BQ6zc/gCbFNrfLQb4woltIaNwc5fnwGUvjpg9xwPbQ2dHSdtBEIzLD0A7GJEqZGY9b8JN2t1fUtVEBdMsVIC5srMUHLjPjONH3NF9YhtBuXyn6exzOMR/P8+MEP+IyM94mowdLHRE5XJMxKAAhM+Pt+pYPtPhkG89r8dX9oZD3xmqoAMI4tlFdPGlphjbm21xvQt4ZD2P9Ojw7dwBTUrru0oKXrsqil7I07vLjZbe5r+tdSnlcEFf8uW9Slsat+jR89RW1NZRoRKMWESFL4966xTSoQeu5J3lMyr0SsjDu1M0jT3ilniNF9Gv+xe9Aq1qdsZCFcX5fFF2Jw03aRi4k/Njc7nFFqsjC+GIZT+g8ruo2EfIG1GmrOj+MTBbGk35dNp3iEGJ90dI0hxiD+pgDMFkY0/WH/DiNDwkzmpyhPkP9tgj89MU1xTjP+m7rAH3XtoPOXCJbYKHuECPVF31+IdtAT6/NoDBX3jG2dkNVZhzozbOqmGoD/YXQY9LYSWQLzHglagpqDmdz9UDU1CksCCJbYMYLMcSzxteE+GvhcU4n7tkCM94e2l6DKgfoIaQ02r7wQMm98VXPB02gp3WZFgNQWDj2xr19J0EDeoBLuPZpVRINjGaZLeyMo4GZdCtUoB9xxG1SMgm8tswW0PgjfVCXpBm+wQb6DT+bpD3IAGW28LssReRooe+kAKwN9F9xDEomgZEu7iIPdUPz1bL+LRgq0N+Ae1aGTdWle2Q5+kx9LrzEoG0BDtBfaUWP9yf1YXI7YtkCm/Gy7aaasQf0Ta4Z6S+JfAvLbIEZ78b4xhOn8RuDy7s+Q83j8cZPGOouF2UPb1PeaP59vfrIaf8nq9ohdozI9gJTY2PInPf5jB1iR7vyZV6xN2R6q955nxvbxI4+7foAxYqvIQve58YFIVR6zaemcbGMLVnw/jeNBzCNLVnwvjA2S/C5VDOZxrbMeZ8bO8Te5d3GmLEj5zy04H1uTDdWJ3bMuwxjR+a8XzFWif2CKbtu7MiC99H44/MXTRVsYu8o0VSNHRkE7//5LF4SNrFfAra0hjrAOrLkfRZqh9gpLPYGAroseJ8ZO8ROxsnaQAB0WfA+M/aIPSe4c0p6cmTLnPe58bbGDvl5Gh8SZjQ5Q32G+m0R+OmLywF6p4DvyM8DvVPAd2R4HuidAr4jvwD0TgHfkV8AeqeA78gvAT0+1XoBn555XX4B6LFno4Dvyc8DPfbsFPAN+XmgR1+ngG/IFaC/F19soEffdiNgPKw1SxZAz/ZqG+izk1nAt+WXgN4u4N/KI1p9/wWg9wr4dn3/BaB3CviODCfQF2lq7Ql51zn2OL3LpNbvTzWmP3S5Zb1axJ45d/+P8kxh2Iw//D7TvX4N/Ud528JfJSw/znLD3LMAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}7 & 5 & 6 & 6\\\\2 & 5 & 8 & 7\\\\5 & 4 & 4 & 8\\\\5 & 2 & 2 & 8\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡7  5  6  6⎤\n",
       "⎢          ⎥\n",
       "⎢2  5  8  7⎥\n",
       "⎢          ⎥\n",
       "⎢5  4  4  8⎥\n",
       "⎢          ⎥\n",
       "⎣5  2  2  8⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Permutation of A in order to obtain L@U:\n",
    "temp_A = list(A)\n",
    "for i,j in enumerate(piv):\n",
    "    temp_A[i], temp_A[j] = temp_A[j], temp_A[i]  # swap row i with row j=piv[i]\n",
    "    \n",
    "print(\"Product of L@U:\")\n",
    "display_matrix(L@U)\n",
    "print(\"The matrix A after the permutations indicated by piv:\")\n",
    "display_matrix(np.array(temp_A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solution\n",
    "\n",
    "The solution of the linear system of equation can be obtained by feeding the output of lu_factor into the method **lu_solve** (which is a wrapper of the LAPACK function getrf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The value of x is:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACoAAABkCAMAAADHT1J1AAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRO/NZt2JIrt8bAICnK4AAAAJcEhZcwAADsQAAA7EAZUrDhsAAAJYSURBVFgJ7ZjRlqMgDIYjIDMrCrib93/WJUEEPRLaPXsxF+VCW/mbxiR8RGFCHgq6w2cFwITapDF3lbDRvEKSTn1VndkkqdVhC6eZm3R2tpqBNflkfTyutFKrtMJGujkShfVBmi6ZVrpzTOZyqbV6lyJLFzxCI0gtavrrBQ2dAARpxECKKZ8GUrb6ivQNByDf1vzCbcG+57s5Qi3cFuQUaDkFFnneU2JdKYKLVa0crmpLf+s5+jak0ivKe1zJud64WO2J8vWP9GdFIGq972fOU6l0kWGpRgIuZ377yFBcw45BQXIBGc6TQJ2wE5CxMkeqVEZGsuoLeAYrlmDCoEi/GSCjLBVyOmYQ9ZABmolCSnhy4Au/eC4dTFXCHRm/L4CfSRkLpiVkTGzTFKmAjOhDCFqlRAyRseaNjHL2QUbJXnt+c3Fb02ysrZ3r55i2uv/RD1ytvunr6UBM5bI3LUwfGZFqcEbaYfLoI8PwFolMA9IKyFhcYtBRqiQVkEHTrQMDZFhfMPC4Ys8IQAxen5kbIiOQwzxGyABwpSt7RMb3r2yGV/Veuoo7Mv5813JBZmCVCshw3IB4fAEZG6V0pg5siAxYdNpiOAAfZOQquB7/dRlerdy//TyreUPKfvaRwfOeKzFL+8igeZOLlqUCMtJ83NZqVUZGgEYqImOLjfRpxZ7IsGmtVKsiMohBrZQ38scuY6bGrUoFByyDrUrvyGjKZVY00Knj+WrwYJKQVeMqdBmcIiTpGBmpz9wR6cnkgwyO2+2Qs/XyWwdL7xSMKZvQzRh95bcOxsBfLMkuJFdQyOYAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}1.0\\\\2.0\\\\3.0\\\\4.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡1.0⎤\n",
       "⎢   ⎥\n",
       "⎢2.0⎥\n",
       "⎢   ⎥\n",
       "⎢3.0⎥\n",
       "⎢   ⎥\n",
       "⎣4.0⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = scp.linalg.lu_solve((LU, piv), b)\n",
    "print(\"The value of x is:\")\n",
    "display_matrix(x.round(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec2'></a>\n",
    "## Jacobi method\n",
    "\n",
    "This is an iterative method.   \n",
    "The idea is to write the matrix $A = D + R$ as the sum of a diagonal matrix D and a \"reminder\" R (a matrix with only zeros in the diagonal).\n",
    "\n",
    "$$ (D+R) x = b $$\n",
    "\n",
    "$$ Dx = b - Rx$$\n",
    "\n",
    "Using this idea, we can start with a guess solution $x^0$ and then iterate until $||x^{n+1}-x^n|| < \\epsilon$.\n",
    "\n",
    "$$ x^{n+1} = D^{-1} (b - Rx^{n})$$\n",
    "\n",
    "The Jacobi and Gauss–Seidel methods converge if the matrix A is strictly diagonally dominant:\n",
    "\n",
    "$$ |a_{i,i}| > \\sum_{i \\not= j} |a_{i,j}| $$\n",
    "where $a_{i,j}$ are the entries of A.\n",
    "\n",
    "For more info on the Jacobi method have a look at the wiki page [link](https://en.wikipedia.org/wiki/Jacobi_method).\n",
    "\n",
    "#### New example:\n",
    "Since the matrix A in the previous example is not diagonally dominant, let us define a new A, and therefore a new b as well. I like x=[1,2,3,4]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matrix A with Rank:  4\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJcAAABkCAMAAAC8R1L8AAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRO/NZondIrt8bFiOv0QAAAAJcEhZcwAADsQAAA7EAZUrDhsAAATpSURBVGgF7ZvrmpwgDIZR1N3qeGq5/2stCYpEkISdbmfaR3/MKsbkNTDKR2ZVZXCr1Xtss8NRqjK6sVv7HlhqBJjaAFf1JkgHxnjB1elpnDjcR9OpSgvyvGi9rpy3A8runbnaocPzDxutmxdiGx0MMBKmqDlq6GDoTqaPTpwa9uC2mXB1ta4Nco0DXDI9TheeDh9T3TDoeEWNPgf0eXJxHB7BoY1w2ePGca345WzdwXHtaU/6DR5muLBmh/EWHIwvuAyG7E1+7Ei5Hnh/f4CrMxqge9PAn8utHsdGi4fzzGTfdxbES+drcaO5Ygb1asfNwg9nd18948xasf24uHxxXBhwxcHjYuc+Z77XWS5hPyKGNpKvpNI4MnLggnwpN+7b/Liv8TEieC5ZmkaAxfejWle4szE/Ugfk0nkjl6EWsBYusWw/Kvdc1fnnqkvBnDdCrgpN2Wcwz6VmeA8N+WfAAtHGgcuCTdQ8TZOu2S/IJZeuB/OoRxuum+xcI49l4+mp1u59imm5+ni4GRXDdQSH272YT1yF+FvtN1dZpu983fkqy0CZ9T2+7nyVZaDMOjG+CkWoRNyKXC725b7uQifmKhWhAnErcomzk9bApMFuMVepCBWIW5HLBifHZpt0xFylIpSXE0rksh/sOkFntllmzFUqQgVcUpdKZfrR9a9chIrFLe/SrtTsqzBxvpCrQIRKxS3vcplmP/u94CoToUoJxK3ApV1BgkEGW5qrUIQqxYtbiUvLM2zrb8D1YT6Q0X+UiFChuOVdOnW5blL6Z0J3FIlQmbgVuDSoj3euRD+WiVCRuJW4HFDiz9tqR8xVKEIl4lbkcoQ3ULuvuMVcpSJUIG5lLnttV633teGYy4/+l+7cXGXpv/N156ssA2XWML46WxZ7t22xFY1/qf746vz9l8+Jrmkmvb/R2AwT3cpZJ/MV1E2z1+Mcp94l8pXp5o3q1qRxEDfmonXT5PVbY4vaeMnWFQ5vVLfGfg9LOBdz2cZgfT++/miZcCbXMcvyuzeqWw8vwV4Q9xmu3oBGm7YVhcA/3Q2iHbqVmmxHgeUzXLZibUYWa88XxD5067dy2d8emAf7rvBZCHXrt3K1dT+bXYkmI0Gj57L7XrcmrQPLZ/pxhCWTxmxKNBkJGoNoh25NWgeWwPXx+YOaBefpCXrkKnxVvqjruahupZ62oyDur8/Eezs4n7zeNXb2lmDDuqDbTX5u3qhuzVnCuWf6cauIrkxhdOOiuvU7uXr8Lo753374fqS69QtcpG6avH5vrGrNvrcPb0S37i6Cv4elbUz2Y2D8qt2bqyzzd77ufJVloMz6Hl93vsoyUGadGF+i4ipGKVKqLBdRyTGXqLjqsEDX+spcPm4gWa8MiUqOuUTFVfTNKVUPQCWrb6Y7VCXHXKLiKroUKFUfmp8CU5Ucc8mLqxhU2I88F1XJMZe7RUFxFQ0ZpVqQL6qSL7j44ipGZJVqCRdRyRdcouIqBs0r1RIuopLTXMLiKkbdK6weIbnDjy+qkpNcfHEVYwuUqmfkuahKTnEJiqsYT6BU5VwnlZzgkhRXMZ5Aqcq5FFXJMZeouIrxBEq1gIuq5JhLVlzFgJxS9VREsvrW0w5RyTHXyfpFhzdXWeLvfP1P+cL/JxD88K3snr9o7f9f1E74YRMXxr4YT3oZ/r9o06jfM65EVPt6NXsAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}10 & 5 & 2 & 1\\\\2 & 15 & 2 & 3\\\\1 & 8 & 13 & 1\\\\2 & 3 & 1 & 8\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡10  5   2   1⎤\n",
       "⎢             ⎥\n",
       "⎢2   15  2   3⎥\n",
       "⎢             ⎥\n",
       "⎢1   8   13  1⎥\n",
       "⎢             ⎥\n",
       "⎣2   3   1   8⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACUAAABkCAMAAAA2RAn4AAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRCKJu2bd7818bEWMfqcAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAJySURBVEgN7ZfbcuMgDIZlwLRrjt7q/Z91JZEE2xC125nO9CJcuCl8FiD4JRkWlGZg1kobBFjQOmrrDILMQwaZWqZA78wnKsQYk9j0Nub4ePlEBUuvr5jpuRPrS7iZO1EO2Q4WgFx5PO4zaqsbWUAaS7LhFX3DTrakS2ZEoTaxTb0D5Utke7xC2NDJi1cqxGJploCEsovkz8QWRFpcaLYUCmr16oxB/JNozW3163T1KDtnKiVeV556ospYwXDzqp16NfPZrLL9widU7wd59tdmrU3kfyIi3Zc7dPUXA5N2tjUBpOtF3T3zc57YnYfFstwU1ULl0CF3XVEt7NE4ua+aauERnFTVPihVtSZnZ+nuqRqCJIplPX6mWkjlK6q1GGYzvuFbOzAjyoq4XVX79xgxq1CWtKqpVsIRFGKbV+eqlbiaK3tfU22w0XD8eqm2HdT4ZKV50oTeAmWH/82Pzyx+T9sxOkmxmmr9TgFa7oum2sQ3jG+gptqNskZrmmqTZGvmNNXivkRnJXc/z7WU1zkR1VVVrW+Zz1av5lquJDgrblPVvv/hURKO6DFTcXLJtR/vh9M2skeypap2kQxsidVUC4YCkqc96qqFaG3LsK9cKyc1PL6nocGMdPysrUWufK+VecpxxsKl2qFWnlJOCrpeKzM02Ap5Z1u9Vp5SEYSSsVad88/LunLolNTKgl8oT5HxZutWK88ozsN9Rq6VJ9TKRWWnuFYeKd+qEvHXo1YW7Lj61XDDahz0WnmgmnXxaq+Vn1DIMx5qZcaOM/L/MSEaimK9Vp5R3De2q62R4J7fTX3lW9Tz56Zz9wtw3qZ8izoH/wADiy5ZROF4qAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}30\\\\50\\\\60\\\\43\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡30⎤\n",
       "⎢  ⎥\n",
       "⎢50⎥\n",
       "⎢  ⎥\n",
       "⎢60⎥\n",
       "⎢  ⎥\n",
       "⎣43⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANcAAABkCAMAAAAi/VVcAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRO/NZondIrt8bFiOv0QAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAZ4SURBVHgB7VzduqI4EERAzqwK6mze/1k3Cfnr6k7AvthPPXAxpqCqk6KVoyVOdzJ+67vv2O6rna47mWG02/k7bHWzM9Mb5+v0JZayjbniaxqWeSnsIs4VihGSEBdUOkQiYsoOCEkEo6/zdfKym31WTvdHKog4HSgHSEJccskYiYgJOQIkEUx8Tf3QG+9rvjr1cos1EMf95BFJiAm5BEhEXHLTGEkUE19WM66+nv7ieF6B3Y04VS8HSEJccskYiYgJOQIkUVzxZbyvi4nXSMSxOHlEEmJCLgESEZfcNEYSxbKvyQxOfzHjWgZxKl4OkIS45JIxEhETcgRIAiz7epjF6U/rQ9chjsXJI5IQE3IJkIi45KYxkgDXfPl+Fb4oTtXLwWNt8msiV0Al3BDJvqCpHeLSThojCXEi4gCJiJHvMZIAy7669UV4hutGxuJUKpGrpBK2RRVfz6ebb87XecDuKNtUIldFJWyLKr7WP3ID/F3OmHlyO1QitbA9W8VXd3fvo672DeJkvLmERUdhZyK9InJalbApor6G/mpu/Wxnmhb7Xt+/7737P9EZBw/SQya9IHKFVMKmiPqSFvuZ+w5fn9W3o19Hv97hDBzPw3fowv41HP3af67egfmL+kXyRXvuEW+04zZO3WmIgc8+OU6BWJwSSQQL/SL5oq2IWJyk62KienWxv49HVuIeOXIQi1MiiWDui+aL8VNVjkjFOXKi2t2Wfsw5cfhQ1parpmyLuC+aL+6MRK3ZkKh28HUTlpNOC3IQSxoW1VIR97XmBi9GonVfWE5aI3IQS5oQitTWyXxBrrMvinITp37N8zikb2KwnLRE5CCWNGxdIGK+IF/cF4mWvp72i4uHuYTFYDlpjchBLGnYukAk+KIRKOaP4iRuZ+yXJzzvgbdHjhzEoRR9QBJg5gv6yfpNqxeI+BpMuCRiuUKQhshBnIjlAEmAna8f81Mo1hdtjkARF1QyDL56H18t6Ym4R44cxGSeCJBE8V/2PWw7b4xV+WPwdfW+hpcSVdWUbRF7HqrTzeDLvzq7+0uJajvi5KfQ72mLuC9dSJmuGw9nbL7al9cL4Wgz4qz4aq9T8NXMG2uT5ET1MSz94L+k3h+OqqZsigRftaV/1P7D10e1y37F9bvuI/qs7vDVHv3i5+Sd97h+TTZp+bbtYe+o+U33H356/37RdYPki7ZviOutRCbiutIdeSzL8sx56ua8WJ1goV8kX7TzIa6vDpmIZWUIVP0HgbNxNyOs25YajxPMfbXzxjip9KhR5kB19DctmRiMbAaq7dm4L5ov7s9FORMrSefC7gsfSC9Xm2GFz2yOuaXG4xRzX2tOUMsbK2vzu5XK4MuXKJ6HWA1nxuMUM1+Q6+zPoxgTK+HCIi58Tff0hcWWGo8DZr4gX2T5Y1wNf9Qqk6/Hcl8/aLvaWA3nw+OABV/KXNSuRKdMvuzSF/ci8xtWC7vTAx4HzHxBP9mzKxVmA62y9NVdw+8SNufdmM35+vnzT7HI9fX3ei6qvO0zXQ8fPiF+7r5Ftb3Of//g+9523licADZUKkO/jI9Ssy+shtPhcYrZ81Cdi8ZvNvNNpZhc4soCDr6u/lbbewz22TpQjdUp5r7aeSNWJ1gVb8a/y7N7B3V2v0TYGag2ZxN8NfNG4gOBRpkD1cswDE9/OdwXqDZnE3zhaj8SH74+q21Hv45+vcMZOJ6H79CF/Ws4+rX/XL0D8xf1i+SL9twjrrcDmYirSiQiFoXTOC5D/HiN6xT6RfJFy0csTuJ3IhNxVYlExKLQpw59DoiJiPtq543iDGGnWqkSnn0y/IAbYOKNqdwXzRe388nsVK1UCRf/OXRK+TAtwn3RfDGmFjknzT5wpFaqhBfj7nRcUp5PizBfGzkPeimwWqkU9sbM2RYUYb4gX9zMJ7MxtVIrXIy5pW+QoYjgS5duKn8O704LRJoM51NHRuf+cje1HJX5gn5u5pN5KrVSJ5zdDe2jqeSozJfuR+3e3frKVSSqKqG7FdBln/EPGC3CfdF80V7n9/yE3vtCJmJPkv5BImJJM9mg2m3+19ZuQEXcF80XedrpasibWqkSrv3qnuEGachRua//PRe1J6kZcconsbv4a+FcyVEFX828sTLJulutVAlP/RDe9wo5quCrufRPOXj4+pROres8+nX06x3OwPo8/M7/r9IGIG4rIpB3OOPqNfj/r3Icu/8A/7uAeZ2lhNgAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}10.0 & 0.0 & 0.0 & 0.0\\\\0.0 & 15.0 & 0.0 & 0.0\\\\0.0 & 0.0 & 13.0 & 0.0\\\\0.0 & 0.0 & 0.0 & 8.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡10.0  0.0   0.0   0.0⎤\n",
       "⎢                     ⎥\n",
       "⎢0.0   15.0  0.0   0.0⎥\n",
       "⎢                     ⎥\n",
       "⎢0.0   0.0   13.0  0.0⎥\n",
       "⎢                     ⎥\n",
       "⎣0.0   0.0   0.0   8.0⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALgAAABkCAMAAAACEQqhAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRIlmzd0i77t8bBwggJIAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAYRSURBVHgB7Zzpkps6EIVlwExiMMu9ev9njVpCSx+BlirVTJyCHzFtHVrfNPuZnoiH1EsnPmRZDK8QD9kPanl+CLeYibaTBP74FGbPOTPwsZ/mKfgpMPbb6bX3MIpH7/dURr72/b6XZ4fJBM/Owd8KYlxWtwnGbsCsbHSwTf7LtHykk2iSL6dPy5XsuY1OLASXM/B5I930tmqM7ff28z11g/8pRUbeaYpNz0EZMvKx6zsZgIOcge/6wvJ0cowtsP2E61BGvi20XedOqIxcaQdHogKQM3CpSV7SHrUY08ThAuAZ+VtjePCMXE3EwEEego+yJ6yXHAwdxubb4N9unofenW1Zud5ysVUskIfgKA/BV3OiPez5hnGAbFZ3dQiu7mzLymmjl02uNtSntZstys4rjnIOrivuUq1mD7j4JLX6ateHrlopki/u8CqQhxVHeQiOewPjc/BeHheWEnmvS6MTFchDcJSH4MIc/084OX0M5J2+bvorM24OchUOnltEs53I7flAQ5Cdge87KWYnx5hGg2XT4H2pXN1QiHu1V/5MdiUNKy5AzsDNNb6HG5CPA2haNeVbSuXq6YC2crcsnI0G+cLAQc7AxUK3/E1d4EapcVzMEx7RSiDzpipYJl+maeo7dS4XyVXuA/xUzsHHST0v6guzOft9fEE+db2+K5fI3+Yxmi5CJXLRd5t8d/OFnIOf4v2dX97g371f7orfFS+swH2oFBaqmeyueLNSFib6Rypeadlwi0Y9PIGfBMWrzL6qZ7LdvrerVDw7q3itZcMtGrRsAFtUZtfPnk9JD1lm4bMx8ErLBiyanMNTmX3QL2LSvtFidgZeadmARYOWja2U/azM/tqUWXc8ilMKmI2BV1o2YNEcb4XeT7LEx2dldr1VcKjAbAzczFBq2eB7N8bAfYSl2bV8XJylitlj8GLLBi0ajE/Bi7OrrddpMe9XlAmzx+DFlg1aNBifghdnN1tPdKDrBbNH4OWWDe48jI8Z2Ud59mOzzVrkmJ3Av+SXy15j2YBFg5aNy+lXarIb/2W/cqf+Y79KqbNswKJBy8bz2rUqQ0hqo8mDw2z8UKmzbMCiOW4Rl/5RpSG0aVttsc4kZmfga6Vl4/yiU8vG1tl+Vmaf6Wb/JK/+NDsDr7VsvF9U4vDUZn/16rd0+qJylp2B29p8wucN/t176a74XfHCCnz0oTKqJoiPW1Z1Z/oH+lU+qO4ffYxfHyrcg4n2xzgMU29fUdRoTg5+UUaO03E5VjzVlIOZzC86O282ccsmkuMwxtEGKRgGnmnKwcRP7TKt8Ata32AE+kr/KAPDwNVM7Je54MEAiJj0o/7ovKaMHIcxxvRpmCQ4eDCY+SWpYWVy7l5GjsMYY3oAB3kKHF+so8ydlLPnzshxGOMoOwdHeQocPZg49STl2911M3IcxjjOzo5blKfBdbvDdYfQs3st8tKyARJ0dDAGOYXhCYfyFDjuHUw9U5vSIK8sG5BjNoxBTmEIjnIC//r1220VanMODzV80JOOvZCb0+eyoQiHMXYMbiUF8/8v/pDFtODBuIRmZVQ/My26M4NW0vJoOCNXCVMwqUMFPRiCCxdTcbHbXiU0iEKtWsdhjEGuQgYO8gvwUw8GM7/0FWW+smxQXtl/5MFPYTh4uiknInl0/fGQdWbZRPI6/+juEIoK+LNf8EPlZ1mqZr/Bq8rVQHxXvEERq1LcFa8qVwPxXfEGRaxKwSvetIcny8EdnkiOwzxm4G17eCKS9N+ARXL0i3jMwNv28CBJxuFBuXn+9vYSxAy8bQ8PkqiYvRlkDCEchpiBt+3hyYGDw4NyHIaYgZttW/XwIImKw4rjazvIcRjjGLxZDw+QUBiCo8MDchzGOAZv1sMDJDF40m5CAwjjCLxdD08GHPc9yHEYYwRv2MMDJFjxnN2EfhHEAN6yhycHnjGEcBhiDl7XIZTp4cmBg8ODchyGmIG37eFBEhUfV5VThyeSp/uPGHjbHp6IpNJuSvtHDDya6i/+4gb/7p1zV/yueGEFzKGi/7bV/cV84aY/JnP/o41qgqAlaIP4MaaiifX/aDMM4g+j8XLi1Xt9PwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}0.0 & 5.0 & 2.0 & 1.0\\\\2.0 & 0.0 & 2.0 & 3.0\\\\1.0 & 8.0 & 0.0 & 1.0\\\\2.0 & 3.0 & 1.0 & 0.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡0.0  5.0  2.0  1.0⎤\n",
       "⎢                  ⎥\n",
       "⎢2.0  0.0  2.0  3.0⎥\n",
       "⎢                  ⎥\n",
       "⎢1.0  8.0  0.0  1.0⎥\n",
       "⎢                  ⎥\n",
       "⎣2.0  3.0  1.0  0.0⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = np.array([[10, 5, 2, 1], [2, 15, 2, 3], [1, 8, 13, 1], [2, 3, 1, 8]]) \n",
    "print(\"Matrix A with Rank: \",np.linalg.matrix_rank(A))\n",
    "display_matrix(A)\n",
    "x = np.array([[1],[2],[3],[4]])\n",
    "print(\"b:\")\n",
    "b = A @ x\n",
    "display_matrix(b)\n",
    "\n",
    "D = np.eye(A.shape[0]) * np.diag(A)       # diagonal\n",
    "D_inv = np.eye(A.shape[0]) * 1/np.diag(A) # inverse of D \n",
    "R = A - D                                 # remainder\n",
    "print(\"D:\")\n",
    "display_matrix(D)\n",
    "print(\"R:\")\n",
    "display_matrix(R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comment\n",
    "Notice that we do not need to invert the matrix D.  \n",
    "Since D is a diagonal matrix, we can just invert the elements of the diagonal. This is a much more efficient method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Convergence in 56 iterations:\n",
      "x = \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACoAAABkCAMAAADHT1J1AAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRO/NZt2JIrt8bAICnK4AAAAJcEhZcwAADsQAAA7EAZUrDhsAAAJYSURBVFgJ7ZjRlqMgDIYjIDMrCrib93/WJUEEPRLaPXsxF+VCW/mbxiR8RGFCHgq6w2cFwITapDF3lbDRvEKSTn1VndkkqdVhC6eZm3R2tpqBNflkfTyutFKrtMJGujkShfVBmi6ZVrpzTOZyqbV6lyJLFzxCI0gtavrrBQ2dAARpxECKKZ8GUrb6ivQNByDf1vzCbcG+57s5Qi3cFuQUaDkFFnneU2JdKYKLVa0crmpLf+s5+jak0ivKe1zJud64WO2J8vWP9GdFIGq972fOU6l0kWGpRgIuZ377yFBcw45BQXIBGc6TQJ2wE5CxMkeqVEZGsuoLeAYrlmDCoEi/GSCjLBVyOmYQ9ZABmolCSnhy4Au/eC4dTFXCHRm/L4CfSRkLpiVkTGzTFKmAjOhDCFqlRAyRseaNjHL2QUbJXnt+c3Fb02ysrZ3r55i2uv/RD1ytvunr6UBM5bI3LUwfGZFqcEbaYfLoI8PwFolMA9IKyFhcYtBRqiQVkEHTrQMDZFhfMPC4Ys8IQAxen5kbIiOQwzxGyABwpSt7RMb3r2yGV/Veuoo7Mv5813JBZmCVCshw3IB4fAEZG6V0pg5siAxYdNpiOAAfZOQquB7/dRlerdy//TyreUPKfvaRwfOeKzFL+8igeZOLlqUCMtJ83NZqVUZGgEYqImOLjfRpxZ7IsGmtVKsiMohBrZQ38scuY6bGrUoFByyDrUrvyGjKZVY00Knj+WrwYJKQVeMqdBmcIiTpGBmpz9wR6cnkgwyO2+2Qs/XyWwdL7xSMKZvQzRh95bcOxsBfLMkuJFdQyOYAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}1.0\\\\2.0\\\\3.0\\\\4.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡1.0⎤\n",
       "⎢   ⎥\n",
       "⎢2.0⎥\n",
       "⎢   ⎥\n",
       "⎢3.0⎥\n",
       "⎢   ⎥\n",
       "⎣4.0⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x0 = np.ones_like( b )          # initial guess\n",
    "eps = 1e-10                     # tolerance error  \n",
    "N_max = 100                     # max number of iterations\n",
    "\n",
    "for i in range(1,N_max+1):\n",
    "    x_new = D_inv @ ( b - R @ x0 )\n",
    "    if scp.linalg.norm(x_new - x0) < eps:\n",
    "        print(\"Convergence in {} iterations:\".format(i))\n",
    "        break\n",
    "    x0 = x_new\n",
    "    if i==N_max:\n",
    "        print(\"Fail to converge in {} iterations\".format(i))\n",
    "print(\"x = \")    \n",
    "display_matrix(x_new.round(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec3'></a>\n",
    "## Gauss-Seidel method\n",
    "\n",
    "This method is very similar to the Jacobi method. It is again an iterative method.   \n",
    "The idea is to write the matrix $A = L + U$ as the sum of a lower triangular matrix L and a strictly upper triangular matrix U.\n",
    "\n",
    "$$ (L+U) x = b $$\n",
    "\n",
    "$$ Lx = b - Ux$$\n",
    "\n",
    "Like in the Jacobi method we look for an approximate solution $x^{n+1}$, and iterate until $||x^{n+1}-x^n|| < \\epsilon$.\n",
    "\n",
    "$$ x^{n+1} = L^{-1} (b - Ux^{n})$$\n",
    "\n",
    "Again, given the specific form of the matrices L and U, we don't need to invert any matrix.  \n",
    "We will just use the forward substitution.\n",
    "\n",
    "For more info on the Gauss-Seidel method have a look at the wiki page [link](https://en.wikipedia.org/wiki/Gauss%E2%80%93Seidel_method).\n",
    "\n",
    "Let us consider the same A,x,b as in the previous example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {},
   "outputs": [],
   "source": [
    "U = np.triu(A, k=1)             # Matrix U       \n",
    "L = np.tril(A)                  # Matrix L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Convergence in 15 iterations:\n",
      "x = \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACoAAABkCAMAAADHT1J1AAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRO/NZt2JIrt8bAICnK4AAAAJcEhZcwAADsQAAA7EAZUrDhsAAAJYSURBVFgJ7ZjRlqMgDIYjIDMrCrib93/WJUEEPRLaPXsxF+VCW/mbxiR8RGFCHgq6w2cFwITapDF3lbDRvEKSTn1VndkkqdVhC6eZm3R2tpqBNflkfTyutFKrtMJGujkShfVBmi6ZVrpzTOZyqbV6lyJLFzxCI0gtavrrBQ2dAARpxECKKZ8GUrb6ivQNByDf1vzCbcG+57s5Qi3cFuQUaDkFFnneU2JdKYKLVa0crmpLf+s5+jak0ivKe1zJud64WO2J8vWP9GdFIGq972fOU6l0kWGpRgIuZ377yFBcw45BQXIBGc6TQJ2wE5CxMkeqVEZGsuoLeAYrlmDCoEi/GSCjLBVyOmYQ9ZABmolCSnhy4Au/eC4dTFXCHRm/L4CfSRkLpiVkTGzTFKmAjOhDCFqlRAyRseaNjHL2QUbJXnt+c3Fb02ysrZ3r55i2uv/RD1ytvunr6UBM5bI3LUwfGZFqcEbaYfLoI8PwFolMA9IKyFhcYtBRqiQVkEHTrQMDZFhfMPC4Ys8IQAxen5kbIiOQwzxGyABwpSt7RMb3r2yGV/Veuoo7Mv5813JBZmCVCshw3IB4fAEZG6V0pg5siAxYdNpiOAAfZOQquB7/dRlerdy//TyreUPKfvaRwfOeKzFL+8igeZOLlqUCMtJ83NZqVUZGgEYqImOLjfRpxZ7IsGmtVKsiMohBrZQ38scuY6bGrUoFByyDrUrvyGjKZVY00Knj+WrwYJKQVeMqdBmcIiTpGBmpz9wR6cnkgwyO2+2Qs/XyWwdL7xSMKZvQzRh95bcOxsBfLMkuJFdQyOYAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}1.0\\\\2.0\\\\3.0\\\\4.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡1.0⎤\n",
       "⎢   ⎥\n",
       "⎢2.0⎥\n",
       "⎢   ⎥\n",
       "⎢3.0⎥\n",
       "⎢   ⎥\n",
       "⎣4.0⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x0 = np.ones_like( b )          # initial guess\n",
    "eps = 1e-10                     # tolerance error  \n",
    "N_max = 100                     # max number of iterations\n",
    "\n",
    "for i in range(1,N_max+1):\n",
    "    x_new = solve_triangular(L, (b - U @ x0), lower=True)\n",
    "    if scp.linalg.norm(x_new - x0) < eps:\n",
    "        print(\"Convergence in {} iterations:\".format(i))\n",
    "        break\n",
    "    x0 = x_new\n",
    "    if i==N_max:\n",
    "        print(\"Fail to converge in {} iterations\".format(i))\n",
    "print(\"x = \")    \n",
    "display_matrix(x_new.round(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comment:\n",
    "The Gauss-Seidel algorithm converges much faster than the Jacobi method (15 iterations agains 56)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec4'></a>\n",
    "## SOR (successive over-relaxation) \n",
    "\n",
    "\n",
    "This is another improvement of the previous methods.   \n",
    "Let's write the matrix $A = L + U + D$ as the sum of a strictly lower triangular matrix L, a strictly upper triangular matrix U and a diagonal matrix D.\n",
    "\n",
    "$$ (L+U+D) x = b $$\n",
    "\n",
    "$$ w (L+U+D) x = w b $$\n",
    "\n",
    "$$ w L x = w b - w D x - w U x \\pm D x $$\n",
    "\n",
    "$$ (w L + D) x = w b - \\bigl( w U + (w - 1) \\bigr)  x $$\n",
    "\n",
    "As in the previous methods, given the specific form of the matrices, we don't need to invert any matrix.  The iterative method is the following:\n",
    "\n",
    "$$ x^{n+1} = (w L + D)^{-1} \\biggl( w b - \\bigl( w U + (w - 1) \\bigr)  x^n \\biggr) $$\n",
    "\n",
    "For more info on the SOR method have a look at the wiki page [link](https://en.wikipedia.org/wiki/Successive_over-relaxation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "metadata": {},
   "outputs": [],
   "source": [
    "D = np.eye(A.shape[0]) * np.diag(A)       # diagonal\n",
    "U = np.triu(A, k=1)                       # Strict U       \n",
    "L = np.tril(A, k=-1)                      # Strict L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Convergence in 34 iterations:\n",
      "x = \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACoAAABkCAMAAADHT1J1AAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRO/NZt2JIrt8bAICnK4AAAAJcEhZcwAADsQAAA7EAZUrDhsAAAJYSURBVFgJ7ZjRlqMgDIYjIDMrCrib93/WJUEEPRLaPXsxF+VCW/mbxiR8RGFCHgq6w2cFwITapDF3lbDRvEKSTn1VndkkqdVhC6eZm3R2tpqBNflkfTyutFKrtMJGujkShfVBmi6ZVrpzTOZyqbV6lyJLFzxCI0gtavrrBQ2dAARpxECKKZ8GUrb6ivQNByDf1vzCbcG+57s5Qi3cFuQUaDkFFnneU2JdKYKLVa0crmpLf+s5+jak0ivKe1zJud64WO2J8vWP9GdFIGq972fOU6l0kWGpRgIuZ377yFBcw45BQXIBGc6TQJ2wE5CxMkeqVEZGsuoLeAYrlmDCoEi/GSCjLBVyOmYQ9ZABmolCSnhy4Au/eC4dTFXCHRm/L4CfSRkLpiVkTGzTFKmAjOhDCFqlRAyRseaNjHL2QUbJXnt+c3Fb02ysrZ3r55i2uv/RD1ytvunr6UBM5bI3LUwfGZFqcEbaYfLoI8PwFolMA9IKyFhcYtBRqiQVkEHTrQMDZFhfMPC4Ys8IQAxen5kbIiOQwzxGyABwpSt7RMb3r2yGV/Veuoo7Mv5813JBZmCVCshw3IB4fAEZG6V0pg5siAxYdNpiOAAfZOQquB7/dRlerdy//TyreUPKfvaRwfOeKzFL+8igeZOLlqUCMtJ83NZqVUZGgEYqImOLjfRpxZ7IsGmtVKsiMohBrZQ38scuY6bGrUoFByyDrUrvyGjKZVY00Knj+WrwYJKQVeMqdBmcIiTpGBmpz9wR6cnkgwyO2+2Qs/XyWwdL7xSMKZvQzRh95bcOxsBfLMkuJFdQyOYAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}1.0\\\\2.0\\\\3.0\\\\4.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡1.0⎤\n",
       "⎢   ⎥\n",
       "⎢2.0⎥\n",
       "⎢   ⎥\n",
       "⎢3.0⎥\n",
       "⎢   ⎥\n",
       "⎣4.0⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x0 = np.ones_like( b )          # initial guess\n",
    "eps = 1e-10                     # tolerance error  \n",
    "N_max = 100                     # max number of iterations\n",
    "w = 1.4                         # relaxation factor\n",
    "\n",
    "for i in range(1,N_max+1):\n",
    "    x_new = solve_triangular( (w*L + D), (w*b - w*U@x0 - (w-1)*D@x0), lower=True)\n",
    "    if scp.linalg.norm(x_new - x0) < eps:\n",
    "        print(\"Convergence in {} iterations:\".format(i))\n",
    "        break\n",
    "    x0 = x_new\n",
    "    if i==N_max:\n",
    "        print(\"Fail to converge in {} iterations\".format(i))\n",
    "print(\"x = \")    \n",
    "display_matrix(x_new.round(2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By changing the *relaxation parameter* we can see how the number of iterations changes.  \n",
    "The relaxation must satisfy $0<w<2$.  \n",
    "It seems that in this case the optimal value of $w$ is 1 and the minimum number of iterations is 15.  \n",
    "The Gauss-Seidel algorithm is a special case of the SOR algorithm, with w=1!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sec5'></a>\n",
    "## Thomas algorithm\n",
    "\n",
    "The Thomas algorithm is a simplified form of Gaussian elimination that is used to solve tridiagonal systems of equations.   \n",
    "Have a look at the wiki page: [TDMA](https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm)\n",
    "\n",
    "The idea is simply to transform the matrix A into a triagonal matrix and then use backward substitution.\n",
    "\n",
    "Here we consider the single coefficients of the linear system instead of performing matrix operations like in the previous algorithms.  \n",
    "\n",
    "$$ a_i x_{i-1} + b_i x_{i} + c_i x_{i+1} = d_i $$\n",
    "\n",
    "The algorithm we implement is the following (taken from wikipedia):\n",
    "\n",
    "$$ w_i = \\frac{a_i}{b_{i-1}} $$\n",
    "\n",
    "$$ b_i = b_i - w_i c_{i-1} $$\n",
    "\n",
    "$$ d_i = d_i - w_i d_{i-1} $$\n",
    "\n",
    "and then we perform backward substitution:\n",
    "\n",
    "$$ x_n = \\frac{d_n}{b_n} $$\n",
    "\n",
    "$$ x_i = \\frac{d_i - c_i x_{i+1}}{b_i} $$\n",
    "\n",
    "Let us introduce a new matrix A wich is tridiagonal. Of course we like to keep the ```x = [1,2,3,4]``` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matrix A with Rank:  4\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJcAAABkCAMAAAC8R1L8AAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRO/NZondIrt8bFiOv0QAAAAJcEhZcwAADsQAAA7EAZUrDhsAAATvSURBVGgF7ZvtmpwoEIVpUCfR9ivh/q91obRbS6DqMLtJJvvIj8TWY/FaoHLAMQ9PxZqvUeYNx5iHd00o7dfAMmOEsT5yPb4I0oExFrg6N42ThvtsOvNwQJ6haEx05WqHjqCfobZuXo4LyG0NsSdMuSOXfVA0JmJcnXXWE9c4xMjT8xL/8vM52UZBpzOgaFzEuEKQZuNa6eZstx8XmuMnegdD0biowOWpyt7LfQflgqJxUZ6r8y5mpffNkZzMlh3Hxmk3R+inSLSLKM+1bL35oXTqNfTFxfcZYrYLinYRlbgoXxoX1b7ODCLzY9nyJUe7iPJcl6Rm6jp2Oa/dklC0iyjPZbZO2Mr93tJjZNIbEorGRQWudY0JGeXnxEBcThbFOFA0LipwbQ85Jz9XqQ+aWRbR9dFTWonGqyxwmTm+hwb5GbBEsHHQulcQIdG4iHM5O/inHUOkbgpjDRkrPCPcZN32Po1pKRcoGhNxrnLk333k5qrL+J2vO191GahT3/3rzlddBurUmf61OLeu2qvxqAUxt1BI0deGN3Z0OcBg7+WAAXOLhSz7WkqCpQHCQEOmIyuXrcMBG8DcQiFlX2vMQD7CqnMpuwM2gImEQsq+1pgnjYv/Uy4opOxrX+01q6P2d74wcxsGrWJIxQ/tXL3iaIPsxQWa22DexZkfxdfuXLPea15cdIZubo0SUvG1G5fbrM4Omf+PcenmVguZaccP/8GrbgCsVzuC5lYPyX3tj8w8ZhuxFs1+7fnCzC0QUva1gehB2VInAneuLbWKuUVCyr42JGqepslZdZZm50LMLRaSmd/0vf3cpvUVrsMBA+YWC3n7Wn6XV/1K27Hq9F8mvrnqUnvnqz5fXVgW+2plCSsaf9P645/O31/U75m/FPPWNc3k1OWhV4gljAdWeYHuJY1LB0n/Yv7yUGa2aIxjtap2A0wDj9bHye5SeVnlcDzl4v6yFCLub6mORVxXOAxwQ4s6vjhOOZQxdsrF/WXUlMpEizVdsab9vH2g1g+hxTsvXcWujOelXNxflpji/t5HjzZJLRNlp9qM3I4nZcJ18SUxcLFY70cV68zVzaKHlLgu/rLIRAcm75/qu+Jd2zLN8qrNW5lpx4u/FLla288+dhuxnGozk6g+Kf9NO47RlDd+/xKkyHaqLcwWSeqTMnJ9fPt+Csr95elAsrmt8D3kRd13v98M6SqpT1w/v12fq9xfJjDHji5cUiy0LrhtZv/da/M0nYNyJe0YFjpjdGVxlQD2FdFVceY710BPu1lapT/lK+Xi66bZDOw7e7oXR/nbj3c7jvE514pfishczF9KXOZhnfrePgxw78JEt3DzHspQaSZfIsrvOnhz1WX6ztedr7oM1Knv/nXnqy4DdepM/6rwtcjHyjAPc8kZLtzX4kpzsqwlUOaSUy7c18JKblkLXNwlp1y4r8WV70FYgSnu5i455cJ9La5EuLhLTrhwX4srQzZOQ9FSzphLTrhwX4srMa7wccThkjNcdF/IHw3TJdc4YCRfzCUnXHjr4EooX9wlJ1zYJ8iUMNwBI/2Lu+SUC/a12OfFdAUA18Ulp1y4r8WVAJfhLjnlwn1thRLo99wlZ7hwXwsrmWXdWzb5j7nkDFdywp/YcXPVZf3O1/8pX/Qph/5pVd0lf1b9/nvRMOCPRZgy+2wVnzqP/l60acw/CIBH0Jz3a3kAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}10 & 5 & 0 & 0\\\\2 & 15 & 2 & 0\\\\0 & 8 & 13 & 1\\\\0 & 0 & 1 & 8\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡10  5   0   0⎤\n",
       "⎢             ⎥\n",
       "⎢2   15  2   0⎥\n",
       "⎢             ⎥\n",
       "⎢0   8   13  1⎥\n",
       "⎢             ⎥\n",
       "⎣0   0   1   8⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACUAAABkCAMAAAA2RAn4AAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRCK7ie/dzWZ8bFsxDvgAAAAJcEhZcwAADsQAAA7EAZUrDhsAAALqSURBVEgN7VfRktsgDMSAydUYY1r9/7dWK8BEDplrH+6hndNMiAOLJIRWcsxCItbMJNdFYxZynmWdgUzBkiWglilgTBaFis7tu2wJLpV0bVaoANcSbTyebD/k2NQplA2YPQ5jCn94w4mRRaGOjCnLbu5y4JVk2w11yixQJKiN2rmVLmgyJlMI5PC0kcfXTZfMbJRM5A/LUr9mqMzGYtX1HuVg6zOLXhxq3q9vvF8BitHsO/wq00iYRTT52KLqplGNOaXkLGKbcUNHv0gVr7PmElAhcb50kL4h+DIVpWuKwOQ3qofm6yIR+bZ3YUTwPjkQU0RZjEiclQqPkkK2UUjfkJfcJM6JFUgTp/m1HWwjEK8lydUgLGa4soj91eJGO9M6icYZKmShoiUqF+iuK6bsam1IRGerEjOLCc6Z1W6Z5AlevPrFlSmYgmriiZ9EgHrQo/5gIrLsHI9DnpbO2p/PFZOEpYwKPAux7ZDK4iFRyhSbLrO3kqlQBVtXFK1Nzlda+bp5vzku0nLHi3Xvblt8mQ3K4gwgc/8+KvgrRd6cMnJ4/7Y/vlF1u6E/Q53s4eLAQu61yfYmetN1oEiDHeHg9Ig1F/mnvqEzWW4JLFZ6rRWGv6CkdwIF5qKFY2TRujoqkGSt768QN1Qp3qFjHE1XqzoaBdpHvClUW6eUnxeL8MHsqCEn38maOtO0LkE55hA77P2Snv16fPyQZWO1N/xm0dLg18fTbVefXV8zZz+zsliDmFljQdgXMQ0zCiV1VQCeWRvOXuQ0ykSXbK1fifnbgnXXVc8wGZXFyXqd+kb10HxdJEavHfyFVWXxqdde/BXXFGr0WnPx9xU1eq3pWSogbVGmanf/BNV7befvq0WmbO+1F39nKJ6rvVYWhb94UmeUJem19any944avfbGX6Xr6rWt5lz8VajRawd/xbZCjV47+PuKMqPXDv4CpnTJvtnwP6DkP4JOqX7U678ov3VCRlHrCHzLf1HvzW8+EC8UhTX+GwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}20\\\\38\\\\59\\\\35\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡20⎤\n",
       "⎢  ⎥\n",
       "⎢38⎥\n",
       "⎢  ⎥\n",
       "⎢59⎥\n",
       "⎢  ⎥\n",
       "⎣35⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = np.array([[10, 5, 0, 0], [2, 15, 2, 0], [0, 8, 13, 1], [0, 0, 1, 8]]) \n",
    "print(\"Matrix A with Rank: \", np.linalg.matrix_rank(A))\n",
    "display_matrix(A)\n",
    "print(\"b:\")\n",
    "x = np.array([[1],[2],[3],[4]])\n",
    "b = A @ x\n",
    "display_matrix(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## TDMA:  Tri-Diagonal-Matrix-Algorithm solver\n",
    "def TDMA(A,B):\n",
    " \n",
    "    N = len(B) # number of equations\n",
    "    a = A.diagonal(offset=-1).copy()  # need to copy because the array is not writable\n",
    "    b = A.diagonal(offset=0).copy()   # see A.diagonal(offset=0).flags \n",
    "    c = A.diagonal(offset=1).copy()\n",
    "    d = B.copy()\n",
    "    x = np.ones_like(d)  # initialize x \n",
    "    \n",
    "    # Overwright coefficients\n",
    "    for i in range(1, N):\n",
    "        w = a[i-1]/b[i-1]\n",
    "        b[i] = b[i] - w * c[i-1]\n",
    "        d[i] = d[i] - w * d[i-1]\n",
    "        \n",
    "    # backward substitution \n",
    "    x[-1] = d[-1]/b[-1]\n",
    "    for i in range(N-2, -1, -1):\n",
    "        x[i] = (d[i] - c[i]*x[i+1] )/ b[i]\n",
    " \n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x= \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABoAAABkCAMAAACvvNBNAAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRO/NZiK7id18bOsX9RIAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAGRSURBVEgN7ZbbcoQgDIaR025FEdq8/7OWEHCNhJledKa92Fw4ykci/CRRtUA1rU7LNKLUAsYWcydROz5rQLS8Rl93+8+RC777MS+vjQYZlfn2jU7ljQ6w6Z1UZBp2Yd+IK4FPXKhozHH0VGbIY0UkWCkCQ7qmZwgCChkHdTsj5rXV1BURBcott5lXRSskmjKi3At3QMaQD235AY/+qOxJ1CcvWIckxjqTB1yqjxVQzCklo+vGb/Ju1EkkdC5GeBdjfBl/irw9+9J1HdH+Rj+8RORbjkXeo/dghiKeiQOpA1hAD5CUX0PJaQ8bvZAFrENyQEQ+tzS8FZGKKZsuzRgw4QvRRqRC+wogejw/aB4l4FG3oNTX8yIv1BrpiAUMB3pnEFJ0R4kclHNCY15qNaUFtAXeEE1vV+71RpdDqWJQgyi3o1C5noyELB2agOK+zbySmqE9zpAvZz/xwpyWkcMPg4h8LQQROY0GQVO2jWpMt1zqdbrlA6B96ceAZZFk/wpN/yo9/kNa20uuLL3+VVqrvgFO7SD82wBfAQAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}1\\\\2\\\\3\\\\4\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡1⎤\n",
       "⎢ ⎥\n",
       "⎢2⎥\n",
       "⎢ ⎥\n",
       "⎢3⎥\n",
       "⎢ ⎥\n",
       "⎣4⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_sol = TDMA(A,b)\n",
    "print(\"x= \"); display_matrix(x_sol)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Ok, I know the previous code is not very efficient! \n",
    "I just wanted to prove that the algorithms work. \n",
    "\n",
    "In order to show how slow is python code compared to routines that call optimized C/fortran, I wrote a wrapper for the LAPACK function **dgtsv** which implements the Thomas algorithm. I called this function *\"Thomas\"*.  \n",
    "You can see that it is about **8 times much faster** than the function I wrote above.\n",
    "\n",
    "On this topic, I wrote an entire notebook about code optimization and speed up. (notebook **A2**)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution of Ax=b with the LAPACK function dgtsv. \n",
      " x= \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACoAAABkCAMAAADHT1J1AAAAP1BMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADFBd4eAAAAFHRSTlMAMquZdlQQQO0wRO/NZt2JIrt8bAICnK4AAAAJcEhZcwAADsQAAA7EAZUrDhsAAAJYSURBVFgJ7ZjRlqMgDIYjIDMrCrib93/WJUEEPRLaPXsxF+VCW/mbxiR8RGFCHgq6w2cFwITapDF3lbDRvEKSTn1VndkkqdVhC6eZm3R2tpqBNflkfTyutFKrtMJGujkShfVBmi6ZVrpzTOZyqbV6lyJLFzxCI0gtavrrBQ2dAARpxECKKZ8GUrb6ivQNByDf1vzCbcG+57s5Qi3cFuQUaDkFFnneU2JdKYKLVa0crmpLf+s5+jak0ivKe1zJud64WO2J8vWP9GdFIGq972fOU6l0kWGpRgIuZ377yFBcw45BQXIBGc6TQJ2wE5CxMkeqVEZGsuoLeAYrlmDCoEi/GSCjLBVyOmYQ9ZABmolCSnhy4Au/eC4dTFXCHRm/L4CfSRkLpiVkTGzTFKmAjOhDCFqlRAyRseaNjHL2QUbJXnt+c3Fb02ysrZ3r55i2uv/RD1ytvunr6UBM5bI3LUwfGZFqcEbaYfLoI8PwFolMA9IKyFhcYtBRqiQVkEHTrQMDZFhfMPC4Ys8IQAxen5kbIiOQwzxGyABwpSt7RMb3r2yGV/Veuoo7Mv5813JBZmCVCshw3IB4fAEZG6V0pg5siAxYdNpiOAAfZOQquB7/dRlerdy//TyreUPKfvaRwfOeKzFL+8igeZOLlqUCMtJ83NZqVUZGgEYqImOLjfRpxZ7IsGmtVKsiMohBrZQ38scuY6bGrUoFByyDrUrvyGjKZVY00Knj+WrwYJKQVeMqdBmcIiTpGBmpz9wR6cnkgwyO2+2Qs/XyWwdL7xSMKZvQzRh95bcOxsBfLMkuJFdQyOYAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}1.0\\\\2.0\\\\3.0\\\\4.0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡1.0⎤\n",
       "⎢   ⎥\n",
       "⎢2.0⎥\n",
       "⎢   ⎥\n",
       "⎢3.0⎥\n",
       "⎢   ⎥\n",
       "⎣4.0⎦"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"Solution of Ax=b with the LAPACK function dgtsv. \\n x= \")\n",
    "display_matrix(Thomas(A,b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "40.2 µs ± 735 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit \n",
    "TDMA(A,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.58 µs ± 91.7 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit \n",
    "Thomas(A,b)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}