Linear_systems_3.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "toc_visible": true,
      "authorship_tag": "ABX9TyOXPvbhaXScY2V34w1nuBQy",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/mwelland/ENGPYHS_3NM4/blob/main/Linear_systems_3.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Goals:\n",
        "- Show how to compute matrix decompositions\n",
        "- Apply decomposed matrix to solve linear systems\n",
        "- Explain matrix types and decompositions\n"
      ],
      "metadata": {
        "id": "ItDCQPlMWWyY"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# LU decomposition\n",
        "\n",
        "Commonly we will have to repeatedly solve\n",
        "$Ax = b$ for multiple $b_i$. Gauss Elimination for each $b_i$ would be grossly inefficient. If you knew all the $b_i$ in advance you could do this in parallel by forming the augmented matrix:\n",
        "\n",
        "$[A|b_1 \\ b_2 \\  b_3 \\ ...]$\n",
        "\n",
        "but this is seldom the case.\n",
        "\n",
        "It is much more efficient to decompose the matrix $A$ into a form that is easier to solve.\n",
        "\n",
        "> There are other reasons to do this for special matrix types and distributed computing which we will discuss later.\n",
        "\n",
        "We have actually already seen this efficiency boost with back-substitution. The equation $U x = b$ solves in $O(n^2)$.\n"
      ],
      "metadata": {
        "id": "KwZZirxA8hJv"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Any square matrix can be decomposed,\n",
        "\n",
        "$A = LU$\n",
        "\n",
        "where:\n",
        "\n",
        "$L$ is a lower triangular matrix\n",
        "\n",
        "$U$ is an upper triangular matrix\n",
        "\n",
        "Now, the linear system becomes:\n",
        "\n",
        "$$\n",
        "\\begin{eqnarray*}\n",
        "Ax &= b \\\\\n",
        "LUx &= b\n",
        "\\end{eqnarray*}\n",
        "$$\n",
        "\n",
        "Now let $y = Ux$, such that\n",
        "\n",
        "$$\\begin{eqnarray*}\n",
        "Ly &= b \\\\\n",
        "Ux &= y\n",
        "\\end{eqnarray*} $$\n",
        "both of which solve in $O(n^2)$.\n",
        "\n",
        "NOTE: L and U are generally *not unique*.\n",
        "\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "CThm-lm88mwc"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Example: Return to the previouis example:\n",
        "\n",
        "\\begin{eqnarray*}\n",
        "4x_1 + 3x_2 - 5x_3 &=& 2 \\\\\n",
        "-2x_1 - 4x_2 + 5x_3 &=& 5 \\\\\n",
        "8x_1 + 8x_2  &=& -3 \\\\\n",
        "\\end{eqnarray*}\n",
        "\n",
        "Through Gaussian Elimination, we found\n",
        "\n",
        "$$ U=\n",
        "\\begin{bmatrix}\n",
        "4 & 3 & -5 \\\\\n",
        "0 & -2.5 & 2.5 \\\\\n",
        "0 & 0 & 12 \\\\\n",
        "\\end{bmatrix}$$\n",
        "\n",
        "by clearing the first column by multiplying the first row by $-0.5$ for the second row, and  $2$ for the third. The second column was cleared with the second row multiplied by $-0.8$. These coefficients turn out to be the elemements of the $L$ matrix (with 1's along the diagonal)!\n",
        "\n",
        "$$ L=\n",
        "\\begin{bmatrix}\n",
        "1 & 0 & 0 \\\\\n",
        "-0.5 & 1 & 0 \\\\\n",
        "2 & -0.8 & 1 \\\\\n",
        "\\end{bmatrix}$$\n",
        "\n",
        "Let's verify:"
      ],
      "metadata": {
        "id": "QmuC4gf-vBeD"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# prompt: Do decomposition on the above matrix\n",
        "\n",
        "import numpy as np\n",
        "\n",
        "# Define the matrix A\n",
        "A = np.array([[4, 3, -5],\n",
        "              [-2, -4, 5],\n",
        "              [8, 8, 0]])\n",
        "\n",
        "print(\"The original matrix A is:\\n\", A, \"\\n\")\n",
        "L = np.array([[1,0,0],\n",
        "               [-.5, 1,0],\n",
        "               [2,-.8,1]])\n",
        "\n",
        "U = np.array([[4,3,-5],\n",
        "              [0,-2.5,2.5],\n",
        "              [0,0,12]])\n",
        "\n",
        "print(\"The reconstructed matrix is:\\n\", L@U)\n",
        "\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "CP0Af5Ztv35r",
        "outputId": "6f01e34e-9b3f-4cf2-9cde-18bd66a919f0"
      },
      "execution_count": 23,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "The original matrix A is:\n",
            " [[ 4  3 -5]\n",
            " [-2 -4  5]\n",
            " [ 8  8  0]] \n",
            "\n",
            "The reconstructed matrix is:\n",
            " [[ 4.  3. -5.]\n",
            " [-2. -4.  5.]\n",
            " [ 8.  8.  0.]]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Let's check the package decomposition!"
      ],
      "metadata": {
        "id": "L4wtwinxR3vd"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Calculate the LU decomposition\n",
        "from scipy.linalg import lu, inv\n",
        "P, L, U = lu(A)\n",
        "\n",
        "print(\"Permutation Matrix (P):\\n\", P)\n",
        "print(\"Lower Triangular Matrix (L):\\n\", L)\n",
        "print(\"Upper Triangular Matrix (U):\\n\", U)\n",
        "\n",
        "\n",
        "\n",
        "print(\"\\nMultiply L and U:\\n\", L@U, \"\\nwhich is correct but pivoted!\")\n",
        "\n",
        "print(\"\\nMultiply PLU:\\n\", P@L@U, \"\\nwhich is the original matrix!\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "12yze0ctR2vy",
        "outputId": "3f2b4514-efb0-468c-e9ea-87f27af12759"
      },
      "execution_count": 24,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Permutation Matrix (P):\n",
            " [[0. 0. 1.]\n",
            " [0. 1. 0.]\n",
            " [1. 0. 0.]]\n",
            "Lower Triangular Matrix (L):\n",
            " [[ 1.    0.    0.  ]\n",
            " [-0.25  1.    0.  ]\n",
            " [ 0.5   0.5   1.  ]]\n",
            "Upper Triangular Matrix (U):\n",
            " [[ 8.   8.   0. ]\n",
            " [ 0.  -2.   5. ]\n",
            " [ 0.   0.  -7.5]]\n",
            "\n",
            "Multiply L and U:\n",
            " [[ 8.  8.  0.]\n",
            " [-2. -4.  5.]\n",
            " [ 4.  3. -5.]] \n",
            "which is correct but pivoted!\n",
            "\n",
            "Multiply PLU:\n",
            " [[ 4.  3. -5.]\n",
            " [-2. -4.  5.]\n",
            " [ 8.  8.  0.]] \n",
            "which is the original matrix!\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "NB: $P$ in the above is the permutation matrix that, when multiplied by LU recovers the original matrix. It is *not* the pivoting operation that is done internally (although that matrix is easily obtained!)."
      ],
      "metadata": {
        "id": "h9Ro4bk8y-C1"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Matrix types\n",
        "\n",
        "LU decomposition has useful properties depending on the types of matricies involved."
      ],
      "metadata": {
        "id": "u8XBgRiu4mL-"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "|Symbol | Matrix Type | Example 2x2 | Interesting Properties |\n",
        "| --- |---|---|---|\n",
        "|0| Zero | $\\begin{bmatrix} 0 & 0 \\\\ 0 & 0 \\end{bmatrix}$ | A-A = 0 |\n",
        "|I | Identity | $\\begin{bmatrix} 1 & 0 \\\\ 0 & 1 \\end{bmatrix}$ | A I = A |\n",
        "|D|  Diagonal | $\\begin{bmatrix} d_1 & 0 \\\\ 0 & d_2 \\end{bmatrix}$ | $D^{-1}$ = $\\begin{bmatrix} d_1^{-1} & 0 \\\\ 0 & d_2^{-1} \\end{bmatrix}$ |\n",
        "|U| Upper Triangular| $\\begin{bmatrix} 1 & 2 \\\\ 0 & 4 \\end{bmatrix}$ | $U^{-1}$ is another upper triangular matrix|\n",
        "|L| Lower Triangular | $\\begin{bmatrix} 5 & 0 \\\\ -1 & 2 \\end{bmatrix}$ | $L^{-1}$ is another lower triangular matrix |\n",
        "| P | Permutation | $\\begin{bmatrix} 0 & 1 \\\\ 1 & 0 \\end{bmatrix}$ | Its transpose is its inverse ($P^T = P^{-1}$) |\n",
        "| S | Symmetric | $\\begin{bmatrix} 1 & 3 \\\\ 3 & -4 \\end{bmatrix}$ |  Equal to its transpose ($S = S^T$). |\n",
        "| SPD | Symmetric positive definite | $\\begin{bmatrix} 1 & 3 \\\\ 3 & 4 \\end{bmatrix}$ |  $x^T [SPD] x >0$  - think of a quadratic|\n"
      ],
      "metadata": {
        "id": "xqCj0hQK5A3_"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "##Sparse matricies and sparsity patterns\n",
        "Sparsity refers to the fraction of non-zero elements in a matrix. Sparse matricies have a significant fraction of 0 elements, which can **drastically** reduce storage and computational requirements.\n",
        "\n",
        "Storage: Typically one (conceptually) stores a matrix as a 2D array. With significant sparsity, one can instead store only non-zero entries. There is a trade-off here between access / modification and memory efficiency.\n",
        "\n",
        "Computation: (Fill-in)\n",
        "Fill-in is the phenomena where, during computation, zeros become non-zeros. This disrupts the sparsity pattern and the benefits that accompany it, and is hence something to be avoided.\n",
        "\n",
        "\n",
        "\n"
      ],
      "metadata": {
        "id": "jfwpmcCfFaZN"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Sparsity patterns and matrix\n",
        "\n",
        "Sparsity patterns are condensly viewed as monochrome matricies filled in according to 0 or non-zero.\n",
        "\n",
        "Example: A sparsity pattern from a finite element problem in 2D.\n",
        "\n",
        "![Finite_element_sparse_matrix.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAANwAAADcCAAAAAA+MZBdAAAAAmJLR0QA/4ePzL8AAAAJcEhZcwAAUvkAAFL5Acb0BsEAABa4SURBVHjaxV1NbFvHtb62LNMuHPu5aJPYlYT6T5ssChRxbMey1NqKV0FXtpFVsyhsiQ+WjAJeSdvuDKjxg8VYxlsE8sLojtvGaEwyBZhkYy28aNiNSZOFslD8AAHkFWmK78794/ycM3/3kpkdRere+e7MnfnOOd854zgHC71BtKyT7Q2wZR3Hv0F1wiEtB/7IOXD0abszgLsvOouW/9lpktbqwt92W/53i/sO7iM3qE1mDh48mHkE/tApPP10Kj8AcJWvKpb/mZ+e8drVDfjbjav+d5VnhWfkBm75WaFQeNYAf+j0mjPIoP5ULedPNOxtKRx09N4k8sNeJ7+6PFccVs/bK7eEexVvrbypLGazi8FYv1jNee1RA75A41H8XXE+O7+8hd3J+6ETvJ1DGztonuScmbb2gHBDPFGV/MAHV1xdvvkkre5vLc9HYwCtFvkHL5iBnPOe/mq+Sw0I1p78ieujN8Srj7fpWy8tbQngUl24yeKsPQbeQDq6t1b3sTo2VgXAPbnpPb5kO0IxHBJ85Bpr69vhT3Ph3cjI3dKcNOo+IiPnT/xmwm1V9eIWMuEbQnZgm7uZ9jEGl5/69Ok3rsF/djfY3WXlQ3y7DHaljdkbm8FPp2am7rhkkyq9NtoBL95R9JBse/1uxeA67advna6ZrHuze9aYRb6JT5mFvZm9C71us7Ub/rTZJL2snT5cMuMuqqdfOz2aGV0TwPV65dPvf/sfyZze3azS127dOLau26ul4xPHl8Q/18+fLgcPvLq5qw+xXhMwbr2sviSXqF8Ymxi7X98SwLmvvv3tB3X8oq1rJ8oM1pfb2tuDd+8tSTfLJ661tMHVz50sC0/v3bFj5BKd+stqbeGXSwI4798++O23r1wFuO2XwFP2kFarBKtbq1aZBws9Z6GlBi748AsIXOc/375/uqyYluvHbogdad04PjFGZml5cnx8sgxNPek6kNa0DD9sAeDIC/mWygJa23v5R3dXWF1Gghe5dGT//iMl/orNZju0ZJrdIVJwDpz7jdICahS/uBwu6cy+ECzBr0uFArO+kytemplaUVgywwDX07GA4s3YiGLJLZnhgCMWkIJFxzQqtFYWlrak70h+NfegqLBkCHdaqIiWjHeDtgUVDBA4ydwfvrViNpA44c4UREsmZ8nTshg4lKFy9gp5RHN3V+9ll7f0lzrxEiGLRkZu7u6Dx3r7KT2h8JFDnxf0PhZGjR6teAkVi9aeGOIPYXAei/5yQ1y03TsX7xUrLHHeuLJoQre9S+Rjklt2KRaNtY2P/pDXYPTRIk7zeRicx6J/dhl4mq4779yOxuCK85BgbLpGr4Mb8+vqiRNVikWj1odb+fUJ9dg9dH7vb79htyTger3nV/4IWkAr51eiXi5+2B8DM9OFtNdf5z++tqn1081rH+e/jm4gbKSBTVW8dzGYQVS3UHBdxAJqN9viGJiaLqSVDp+q7Ohxrt2dyqn4BgIFCmwqJ+u6fLdQcL4F9N0r1ZQLWLQWf6R44LZ/+fO0BSIQbsbsoG4gkNeAK7/NWFTbVZ9q4uDcV9+9r+xzwKK1mD/F4Nd9NHV6twH6TBmM1A2Ap0Bz5bBbYxPHvY45UrMWY9G7rdCmXtsz2xS+Y8gx95FMIdaCZ2YbNesBU1+bq6yNZka8jsnA4Sx681q4FjSeCTvG5o3ZDfwjefm/AihYuE701yvASUN9p9oXPBrvdcwx3HJhByFIo8CPCXyT5r5VR0GWYBZNgavc5pjt9vpaA/+oYrz/g0Qt2itz7ynAVRb50IFj9SwpcIb0y7ZpuKZJuGGiagSOsOgiYJ7cD40XYeQkb4I6HkBPmocNduT8oEIn1ZGL7Qfpm6XHbA8aRRFGC5p83vKdC60O4Ilx4IhBuvpCcSEqBqc2h9ortytDAAdelAPnvxK55O9VzrQfycEBcQR/++rHvd07UzPTiGPJdx6BA8YRbt8catCxdC6wro4VWIAD4gg+8VhwFvosuonFCojbb+8CTIFYwk0479qeBZoPM//XMTSvHL2fwSw6oKsw26M5IhwpgAn3+rtLfd7J8WH/Xi/rnZTBwSw6oKswT2fYPRgpgAn3tvfT2GLY4v/Pu9fYhboVuJgPYywadBfDFpaqdYmQZhez9bBgmnevUe87l38FYHPeAfkwzKL/eOU5aFODtrGqPf9o+tJ15G6vS2UX+6bwrOzuLPKBzpWLM1OiL4cBF3o1kDe3eflnKSipopWneHjfyKmq1hALA0z7SYJ2e8/onitNKbjQH4Wsud2NL5MrqeI9IxgDjf+Ig8091BgK4tKiuy6x39KGABtdJIFDGwb3IndvHuDDjAVkorIw5Gkc3b6vCEYYbwWYJdO3gEz0McmazDQ2AReblpgl048jMCMnhgDsjFXYXbaUzshpWDLg+4i/jTa2UVrN4b07l78oNkxZNBMC6EmuaOOaTg/crvvj5b1rpiy61xM5A3hFG9d0euA47QxMC1VqHPyKZq5pPYJqAo7RzrSuHRt7d0lk0XI1Dn5FM9c084TAniSyCsgljy3UqryuRq7GGUQLe2IgV1GCI5Ohfj9wvDPvnVSNM4gW9ATS9+iBYxgqFVMMHe+RbdHuW0DNZjtZj1tmmhsgMgExcyiExTBUKhrcoElp5LCPtTNJmhhHUDAxgBz3JJuSgzDUodCvQWzwNJ2gwFFSB+/DQ8ZVGQtphgtOEKTrMPM+EdT0oVDy5GGCSyiah8BFoRaK88Yj13+WGt57NSNekCf82IrmScaB938yeZQti06z2cqjfOEOBG7lwt0CWZY4rzLCor/cGKSEMqLpomdB5lbKT/vCHQjcm91Vh3hb/G3odqyqQVg0qMZJrXV3/3HoZBXyCYlvk//D0MHle/ogcJVSpFihNzbEWsHUOEkaNQaN0heXiQcw1M68Lj5H58nzWd5VCIGjFCsUJZmEHa9d03wEjUY5ZWPRcUA8Soc/aur7ACFwUQiAccrjq4uWGseoUQIcTuhSGD33fSNRrCAMAbASHxyclhrHqFHPlZMoFQ5kfmUZK2CoLDU3iAM6npaANzpk0Sl4o5FLxBZ86eiBUf2XwEGpLL2yXryzEy8ogDc6YNHTyb3RmCMmvrm2n1rf+8XspJj3awhhY2v6FcjANcC9eAAQohhcRcc7jPBh3PnJ9NECXABAAxz20ANwWp59Wz5sHfAPBgT2D4ODxT301VyOxBEIuDdEZSGzVrLOe3M26QJUP57cLBq/c4kaGRACziX6GLlEzRl8EmnK4DwD5XOSZNABRi40QsIfzk8nB/fJ58VhgpO9FDl2sFJI/9WwlB1kQzFvlWeF0j3PQCkUoBB5aISEPy3eM5TLiGbZlIZ3is+fmzxi68xf2JsZyXaeHj1w8OCBo8JFOrTcgGHmdo3IF9qa4CLBB7cVmOSlBkm3RJnHaPMgD7qgndFzStjGCpbfXobAmWQUb72sVrdEcFDsw7ybdo+Eq83ATUteIK7Y7NxOH1zfNS2NWrn6dJsTgpmD4xYUDc8F/ZJ7RnsMTmXBq7gysGIZ18YJwdlXDKIacbeUggXlaeemFsXyaNv9TrNpt7oQ40sabnBsnwryaIlWqlB4+un0cV1wZ34nTzGTeMYu3mn/IAs3pLuJMzw6q/tDWy7mMQFX6lt15hcraUGKGFd75U/H9cF9spaz8CoT4uzx6K4cXIrJzdGqFA6HJric5b1yCtfOwEZuLpvViyNAxaR8v5BKmuOP3Gru87/8fN+f/7cx1HcuMTnWitwFafWSuTdAcMa1cTTB9V1VfkGES0f2f1ZqiBuit24PEFwngS9aBi7mPH4pi9aPl5wMW2sv2LRLR/YPEByjxtmW0kmhZEj9Au7mDdlqRFBbN45n4lXJp+nBwlY+Mz5QcJQaZ/0dGesVir3I1Djhd1FJl90f/nU+BucXggnAubWq05NdFDJ5wEo22CViNY4mOAN7Ia5XQ2yXxzS4yMThFxTOQIFMnnWF0oW5RKzG0ZyWBmZN/znw+ivNSjaxyUOFh8BKNtD7Tl1RW41DbPR2aq+FvJJNbPJQAhywkg3uhjFT4xDvysqgwGGEiBXgWFWy0d76s4MAF7j5NcAxahxNcO/9d1Z0MkNSjRcPioMAFwDQAGdj/4CGzaD1zwI4mMpW9KT9QlBhe319u0OyoEFwSpGNrIFZ1yzhdtIZHYwrhykB2QQmqcrAkswF4Z1TBsmMwIUpAQG4M5/F1YTUcSNtA0syF2zqfplbOdnwrUtYFMr6nYuEmiK4hr7zCMuYXZma+d0ZD9w7f/vn6wTgdMSkILhIYiuCezSi7fbDMmaJ1OG+B25PJowi2IHTkgFD4CLKKIITtlVOfqPZHo+/sye2mh8rpdhcqVAtWouCi8i+BjhOfqPZtut/OxCDU5kWHme/LlaxVP8XAI4uB9u8OvKInV7UtIwEWHYGdunoAVJjX1DjkLsJwq3m7KFSKgsKXcjX58rMwtCg8vgvBhmZZSvXCPFGQ2oc724dIXmz+7z4OhVwWXEmDjOtM8VqhAC4Ir9vB3uswI4GBM67G0SPMJ9mApNHxmsHmJALgEt+N04eJWGhPQU1e3KLZHnqVPLMLc8L9g9gRiUigiy4lMQTmvIoZ7j0KzG4J55hM/0Qr+RJnS6QAriizqimB07VZ6qa6tDlUQLnFXYz5c4jL8OZn56Z7le9v2QURxBC9YW7bCEFODLsoJxX4CHNjxTZtPIynJ1+jRDi6jOKI3DuwubVfatsCQxY6YBtBW4jf4pjkFCxKvtmkNUr5FmKm8TS24tAyicGrnzi4wqncU8XnEFWr5AhK4Lbqn/2rlhb3kEsGS404X33pnX9FFfSXzu5GQwBaGf11s+eZZ8CFRzom0PvSMExlgwXVPK/4w9jwOomiQ0MAWhn9XbqnP2oG59xorWgLbVkgO/wild+tRomPIAom1LJ6lVvBUTMJLNkgO8ktcruTM3MMOEBRMOTSlavGlyqp8WZSTUGWFQFBmdXa77PeP96RkeCHIGzk6oLDS0ByYFLysW0rZVE8ijeaoKLdwoHfQ0OnHjQ13vW4Di3NTJyum6GFMCZV73Hm7KXjINIXvfLhH5EGtH4isI3NN2+ZJfVy7qxfG90eEpihQL3JiyiJ6/YZoQudN/FV+zPS75gn8+irbJ6mV4+9HOIq6dGRkdH99ymwFXCyDdXay+FikFC9b74GBimkaxem4J99DB+SHKIN69f8rbYKKzug1scCeRFXJXEBEkGMVFqCb7WEeokrhhy+ylx2Ca5WzCMuy2aHEELCpJkkEqLj14C1nFtVzvM9RmiuYWBi2juIMDFh2bh4AzCHVjUYj0oV8RMS47mynJwsRXEtoVxBHIRg0AVxvXDwrTBgvJVA3JIyHJwsbXftgVxBN8xYhBixLh+GN6wreNsy7cGfhEdN0PEjnx3vU5aJ+ZHNJBjaJxUlSa4FJ5leucVDGjkjBVR9MgtZ+/ltH3+XI199PiCNMD1V+pk8iiD4WeC1IlrEAFHb4nPq2Jb+bSnHWcLY0rFueVcfJabqnqUYkI5oBEylGqq2MtpcAqfYkJpg8NZtG2InOfYUT4VAZc5h52fyNwNKINsBQ7PYLQVNwgcO8yEI+DG/93uatyt60pPnILAxepuuuF1NzGCR85sslHj5JzM2Nnv6q7G3dxa4w2kxgm+q0LgtqSn2WpT8/WxCYMK9TS48xW8Ng7r9D99vrEDqHH870jSROsqS5sBS8nO8bY2mhmFjnBQtUcjV128xr4wR11EjeOnuyiTbvMf2rmESelIbqnZvD77XOP/Nrp4jX1hdcFioopEpfA5P5SWNTNb7k8eKkrP5qa2/mQ19jtYUTNucW58lVr6oL+fyE5Vp50tiWrs+1aakcgmnZbV90Yn1Ow6akKE0y+L2po9lcqiePNJGhaQb6UlIs6D4GnpEmfg4ekSZ7+2ZtfeHJKPXM++eKcE3E9ztrGth8MQXFS8E1jB+C0xP/XpP+RJWQlbkqxeEFy7/WAv6F0TK4h02k8PydPpErYkWb3IgrJ+7A+VH8Qet24cu1/nzzg6eb1l2e+6Fq02rLGvBrf9Kv9rUdbR292sLfySP52qumk5boK8BNn4zWrs6/hQAAWR3/ppsFoWQyp9tq6xb6D9ws0hmzMWkmi/knu/QpMXPJeWikyR8vh2BzcbqfbesmHRjtSc6CHn0lIxxY0riztWR26beF7ws3otwcl29J/knNWcLTiEAmPn0lLguPL4CXM18fdOh0XzpxPJ5VGYesngbOM0TSUVK+S75UgZalDZPicUCaJOyBUsGXzkEjJsNYtGRs6geiM2PmobNDHDti67ikdP/eqN/IrI1EQOjgJXeg8U/mEtFm2kxsHV6ZSvpUkXqWacR/58PDxKVBZqv0+32UrGsE1r7OOuPW4jwteJ51fPZQ6kuobge6CZGgd3ypaOTNa0wHXd78fTXSBx9mKmxmHc6YwRwr3+0qpAyHd2mbvypFj43G6G8gYf+FhB/dy5OgZOoyqQUSelbEvyUGBwjPw9+MDHCmonT9awaUnu2eymNb3kJ+RKr9lX41CNkb8HH/hYAaOOFN5qd/HKRloLg7cvXI3zlgyvGatxaO+ceFC3WVGzdNNdREVy+izaoJBgJ/9QUdzvFpWtKfdGF+fu9nMFrVj0vMYhr+mWgKR4mpzmDqkCd4rgip9/0u9z4LaOqWyRtVcS10OhwaFakDTB4ZVs0i0KxYNDX0H/nSsmzVmOxy7MJPGFO6QqKjZyvaTD6L3f85Eap4OOnJ61YrOcab9XOeu8AoUaZzDgfG2pdp8Dayu2qfyYuI7LSa7GGRi4Xtt9ow8usLZi/+DanllXy1koV+MQcFWPFi6LWZLJ2u4P1WrjL2azLfbsrh+7sdPQcfMGapxXVZT0OhOTZY9Gb6cLjgjsJ35uBi5m30R6/6ah4aD31ThnxwmCwe9z1IIyO5LJ7HMuWQe3dndcPTVOdcJD8HfEGz0YcCSp5es/O/910dZpoqvG6RFwe38zkx8iuODBZjLO+Pc7VmPnq3EUa5YfxKhNZvbyFVbi8MYAwTVKn+3PfGA3dhpbQRDEcL/5+294cHF4Y4Dgkpk1GgSMqtSeg77pOTnKc8x4bAV1RpucJWQ2do9ws8ZeNR26tCNBSXxSVZ8ERlITB6S5gA8ldWPV2v0M+eFgAwsD13gk1KRJF1xlcbFfEVIIFWmMHMuj4zgCeraxPOyWLri0g3wR2UNOpaYCpo3iF78PT1rtu15gB5Hs5FPFXhbcjdyAOr3cthE1DhGv7sDniVOhblKIeif4G3XuC+jaKx2atXzo4d3IDbrNxOL9DvFFZyZrwEnwnEiBSj5UhfzsRTaaN+i/a4rqK0iGvyAvodJGVfe2F9mYggMlMLrgYCc2mGRgQPXrqvMKyA20npAKXHlyYvzsdzUqJZaaluXTF+rmk0FhpEXRhzBFFplt5VMac1vVE7dWrRELaPwxsKAw22o3nRzcOPoQpshiMyr5VhCY88RIALcCBpx6cdbKwY2jD3KJewrgAkeMHjj17dJMn01tE0fBNR7d75NZNSHSqNAoVihJA5xAv3yPc3hO7GOaRePEeZgmjxE4ZU3RvsnDPhN/5CKTQXlmk8pZrD1yldsrHY27yUaOY9FFARyHXEPolGIN9DRlVWFkAgIXhdaGC04YkARlV4tzvgWUMJdnkNXr08jlwfPnNmVHDtK7C9bsU3Xdcun/dpCNVCuO4FtAssxHjTxFeTGR0pEzlikBtTNHShgF8r/TsoCc9GtZDMpPYrjyePvQ/wMATb06WoZfWQAAAFx0RVh0Y29tbWVudABGaWxlIHNvdXJjZTogaHR0cDovL2NvbW1vbnMud2lraW1lZGlhLm9yZy93aWtpL0ZpbGU6RmluaXRlX2VsZW1lbnRfc3BhcnNlX21hdHJpeC5wbmcWxZJ5AAAAJXRFWHRkYXRlOmNyZWF0ZQAyMDEzLTExLTA4VDE5OjQzOjEzKzAwOjAwPqHvwgAAACV0RVh0ZGF0ZTptb2RpZnkAMjAxMy0xMS0wOFQxOTo0MzoxMyswMDowME/8V34AAABGdEVYdHNvZnR3YXJlAEltYWdlTWFnaWNrIDYuNi45LTcgMjAxMi0wOC0xNyBRMTYgaHR0cDovL3d3dy5pbWFnZW1hZ2ljay5vcmecvblIAAAAGHRFWHRUaHVtYjo6RG9jdW1lbnQ6OlBhZ2VzADGn/7svAAAAGHRFWHRUaHVtYjo6SW1hZ2U6OmhlaWdodAA4MTbPZQcbAAAAF3RFWHRUaHVtYjo6SW1hZ2U6OldpZHRoADgxNhPJVJYAAAAZdEVYdFRodW1iOjpNaW1ldHlwZQBpbWFnZS9wbmc/slZOAAAAF3RFWHRUaHVtYjo6TVRpbWUAMTM4MzkzOTc5M/q2o+AAAAATdEVYdFRodW1iOjpTaXplADEuMzNLQkIhdFWPAAAAM3RFWHRUaHVtYjo6VVJJAGZpbGU6Ly8vdG1wL2xvY2FsY29weV81NjVmOGZiNTIzZGYtMS5wbmeYx4H2AAAAAElFTkSuQmCC)"
      ],
      "metadata": {
        "id": "7qALRcqjN_cR"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We already have examples of sparse matricies with\n",
        "- diagonal\n",
        "- upper and lower triangular\n",
        "and have seen their importance in controlling complexity and roundoff error."
      ],
      "metadata": {
        "id": "TlY-1ThVN-17"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Banded matricies"
      ],
      "metadata": {
        "id": "Yu7Mb8jaMtdv"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Banded matricies only have non-zero elements parallel to the diagonal. Sometimes they are described by the *bandwidth*, the number of parallel bands.\n",
        "\n"
      ],
      "metadata": {
        "id": "laiYdnSOGlMU"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Tridiagonal (bandwidth 3)"
      ],
      "metadata": {
        "id": "S7dFY1EtP9VT"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# prompt: give me a tridiagonal matrix of dimension 10 and view its sparsity pattern with spy\n",
        "\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Create a tridiagonal matrix of dimension 10\n",
        "n = 10\n",
        "T = np.zeros((n, n))\n",
        "np.fill_diagonal(T, 2)  # Main diagonal\n",
        "np.fill_diagonal(T[1:], -1)  # Lower diagonal\n",
        "np.fill_diagonal(T[:, 1:], -1)  # Upper diagonal\n",
        "\n",
        "print(T)\n",
        "\n",
        "# View the sparsity pattern\n",
        "plt.spy(T)\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 609
        },
        "id": "g3ZyT869PVgo",
        "outputId": "08428b34-861a-4330-bd77-2c7151c0fe3e"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[[ 2. -1.  0.  0.  0.  0.  0.  0.  0.  0.]\n",
            " [-1.  2. -1.  0.  0.  0.  0.  0.  0.  0.]\n",
            " [ 0. -1.  2. -1.  0.  0.  0.  0.  0.  0.]\n",
            " [ 0.  0. -1.  2. -1.  0.  0.  0.  0.  0.]\n",
            " [ 0.  0.  0. -1.  2. -1.  0.  0.  0.  0.]\n",
            " [ 0.  0.  0.  0. -1.  2. -1.  0.  0.  0.]\n",
            " [ 0.  0.  0.  0.  0. -1.  2. -1.  0.  0.]\n",
            " [ 0.  0.  0.  0.  0.  0. -1.  2. -1.  0.]\n",
            " [ 0.  0.  0.  0.  0.  0.  0. -1.  2. -1.]\n",
            " [ 0.  0.  0.  0.  0.  0.  0.  0. -1.  2.]]\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAGiCAYAAADa2tCeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAUOUlEQVR4nO3dfWxV9f3A8c+ljtKYtkEcKKEIM0tQwMcikSYui0RjmJnL4h6CCcP/THmSZFndgg9BqCybIRHHhCyERPAhWZjOhC2kizAmhCpiNJuyxWTrNIAmphcxqaa9vz/2k/36A/Re6Kf3Xni9khuz4z29n3xX7zvnnt5zCqVSqRQAMMLGVHsAAM5PAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQIq6DcyTTz4Z06ZNi3HjxsXcuXPjwIED1R6ppnR3d8ecOXOiubk5Jk6cGHfddVe888471R6rpj322GNRKBRixYoV1R6l5rz33ntxzz33xIQJE6KpqSlmz54dr776arXHqimDg4OxatWqmD59ejQ1NcWVV14Zq1evjgv5alx1GZjnnnsuVq5cGQ899FAcPHgwrr322rj99tvj2LFj1R6tZuzevTs6Oztj//79sWvXrvjss8/itttuixMnTlR7tJrU29sbTz31VFxzzTXVHqXmfPTRR9HR0RFf+cpXYufOnfHXv/41fvnLX8b48eOrPVpNWbduXWzcuDE2bNgQf/vb32LdunXx85//PJ544olqj1Y1hXq82OXcuXNjzpw5sWHDhoiIGBoaira2tli6dGl0dXVVebra9MEHH8TEiRNj9+7dccstt1R7nJry8ccfxw033BC/+tWv4tFHH43rrrsu1q9fX+2xakZXV1f85S9/iT//+c/VHqWmfetb34pJkybFb37zm5Pbvvvd70ZTU1M8/fTTVZyseuruCObTTz+N1157LebPn39y25gxY2L+/Pmxb9++Kk5W2/r7+yMi4pJLLqnyJLWns7MzFixYMOx3iv968cUXo729Pe6+++6YOHFiXH/99bF58+Zqj1Vz5s2bFz09PXH48OGIiHjjjTdi7969cccdd1R5suq5qNoDVOrDDz+MwcHBmDRp0rDtkyZNirfffrtKU9W2oaGhWLFiRXR0dMSsWbOqPU5NefbZZ+PgwYPR29tb7VFq1rvvvhsbN26MlStXxk9/+tPo7e2NZcuWxdixY2PRokXVHq9mdHV1RbFYjBkzZkRDQ0MMDg7GmjVrYuHChdUerWrqLjBUrrOzM956663Yu3dvtUepKX19fbF8+fLYtWtXjBs3rtrj1KyhoaFob2+PtWvXRkTE9ddfH2+99Vb8+te/Fpj/4/nnn49t27bF9u3bY+bMmXHo0KFYsWJFTJ48+YJdp7oLzKWXXhoNDQ1x9OjRYduPHj0al112WZWmql1LliyJl156Kfbs2RNTpkyp9jg15bXXXotjx47FDTfccHLb4OBg7NmzJzZs2BADAwPR0NBQxQlrw+WXXx5XX331sG1XXXVV/Pa3v63SRLXpxz/+cXR1dcUPfvCDiIiYPXt2/POf/4zu7u4LNjB1dw5m7NixceONN0ZPT8/JbUNDQ9HT0xM333xzFSerLaVSKZYsWRI7duyIP/3pTzF9+vRqj1Rzbr311njzzTfj0KFDJx/t7e2xcOHCOHTokLj8r46OjlP+xP3w4cNxxRVXVGmi2vTJJ5/EmDHD31IbGhpiaGioShNVX90dwURErFy5MhYtWhTt7e1x0003xfr16+PEiROxePHiao9WMzo7O2P79u3xwgsvRHNzcxw5ciQiIlpbW6OpqanK09WG5ubmU85JXXzxxTFhwgTnqv6P+++/P+bNmxdr166N733ve3HgwIHYtGlTbNq0qdqj1ZQ777wz1qxZE1OnTo2ZM2fG66+/Ho8//njce++91R6tekp16oknnihNnTq1NHbs2NJNN91U2r9/f7VHqikRcdrHli1bqj1aTfvGN75RWr58ebXHqDm///3vS7NmzSo1NjaWZsyYUdq0aVO1R6o5xWKxtHz58tLUqVNL48aNK33ta18r/exnPysNDAxUe7SqqcvvwQBQ++ruHAwA9UFgAEghMACkEBgAUggMACkEBoAUAgNAiroNzMDAQDz88MMxMDBQ7VFqmnUqj3Uqj3Uqj3X6j7r9omWxWIzW1tbo7++PlpaWao9Ts6xTeaxTeaxTeazTf9TtEQwAtU1gAEgx6ldTHhoaivfffz+am5ujUCic9c8pFovD/snpWafyWKfyWKfynM/rVCqV4vjx4zF58uRTbk/w/436OZh///vf0dbWNpovCcAI6+vr+9KbGI76EUxzc/Nov2Td6u/vr/YIAMMUi8Voa2sr67181ANzLh+LXWgu5L8+AWpbOe/lTvIDkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApDirwDz55JMxbdq0GDduXMydOzcOHDgw0nMBUOcqDsxzzz0XK1eujIceeigOHjwY1157bdx+++1x7NixjPkAqFMV3w9m7ty5MWfOnNiwYUNE/OcGYm1tbbF06dLo6ur60v0/v1c1X26Ub9UD8KU+fw/v7+//0iu+V3QE8+mnn8Zrr70W8+fP/+8PGDMm5s+fH/v27TvtPgMDA1EsFoc9ADj/VRSYDz/8MAYHB2PSpEnDtk+aNCmOHDly2n26u7ujtbX15MPdLAEuDOl/RfbAAw9Ef3//yUdfX1/2SwJQAyq6o+Wll14aDQ0NcfTo0WHbjx49Gpdddtlp92lsbIzGxsaznxCAulTREczYsWPjxhtvjJ6enpPbhoaGoqenJ26++eYRHw6A+lXREUxExMqVK2PRokXR3t4eN910U6xfvz5OnDgRixcvzpgPgDpVcWC+//3vxwcffBAPPvhgHDlyJK677rr4wx/+cMqJfwAubBV/D+Zc+R5M+XwPBqg1ad+DAYByCQwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUFV/sktFTKBSqPcIpXB8NKJcjGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNAiouq9cL9/f3R0tJSrZc/RaFQqPYIdaEW16lUKlV7BOA0HMEAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFBUFpru7O+bMmRPNzc0xceLEuOuuu+Kdd97Jmg2AOlZRYHbv3h2dnZ2xf//+2LVrV3z22Wdx2223xYkTJ7LmA6BOFUrncLemDz74ICZOnBi7d++OW265pax9isVitLa2uuEYI8YNx2D0VPIefk53tOzv74+IiEsuueSMzxkYGIiBgYFhwwFw/jvrk/xDQ0OxYsWK6OjoiFmzZp3xed3d3dHa2nry0dbWdrYvCUAdOeuPyO67777YuXNn7N27N6ZMmXLG553uCKatrc1HZIwYH5HB6En/iGzJkiXx0ksvxZ49e74wLhERjY2N0djYeDYvA0AdqygwpVIpli5dGjt27IiXX345pk+fnjUXAHWuosB0dnbG9u3b44UXXojm5uY4cuRIRES0trZGU1NTyoAA1KeKzsGc6TzFli1b4kc/+lFZP8OfKTPSnIOB0ZN2DsZ/yACUy7XIAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFKc0y2Tzye1eJ01F+AsTy2uUy3+PsFocwQDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEhxUbUH4MxKpVK1RzhFoVCo9gh1oRbXqRZ/nzi/OYIBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKc4pMI899lgUCoVYsWLFCI0DwPnirAPT29sbTz31VFxzzTUjOQ8A54mzCszHH38cCxcujM2bN8f48eNHeiYAzgNnFZjOzs5YsGBBzJ8//0ufOzAwEMVicdgDgPNfxbdMfvbZZ+PgwYPR29tb1vO7u7vjkUceqXgwAOpbRUcwfX19sXz58ti2bVuMGzeurH0eeOCB6O/vP/no6+s7q0EBqC+FUqlUKvfJv/vd7+I73/lONDQ0nNw2ODgYhUIhxowZEwMDA8P+3ekUi8VobW2N/v7+aGlpOfvJqYpCoVDtEThLFfynDmdUyXt4RR+R3XrrrfHmm28O27Z48eKYMWNG/OQnP/nSuABw4agoMM3NzTFr1qxh2y6++OKYMGHCKdsBuLD5Jj8AKSr+K7L/7+WXXx6BMQA43ziCASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEhxztci48JSi/cUcY+a8tTiOtXi7xMjxxEMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASDFRdUeAM5VqVSq9ginKBQK1R6hLtTiOtXi71O9cgQDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkqDgw7733Xtxzzz0xYcKEaGpqitmzZ8err76aMRsAdayiy/V/9NFH0dHREd/85jdj586d8dWvfjX+/ve/x/jx47PmA6BOVRSYdevWRVtbW2zZsuXktunTp4/4UADUv4o+InvxxRejvb097r777pg4cWJcf/31sXnz5qzZAKhjFQXm3XffjY0bN8bXv/71+OMf/xj33XdfLFu2LLZu3XrGfQYGBqJYLA57AHD+K5QquD/o2LFjo729PV555ZWT25YtWxa9vb2xb9++0+7z8MMPxyOPPHLK9v7+/mhpaTmLkaH21eKtgCmPWyZ/sWKxGK2trWW9h1d0BHP55ZfH1VdfPWzbVVddFf/617/OuM8DDzwQ/f39Jx99fX2VvCQAdaqik/wdHR3xzjvvDNt2+PDhuOKKK864T2NjYzQ2Np7ddADUrYqOYO6///7Yv39/rF27Nv7xj3/E9u3bY9OmTdHZ2Zk1HwB1qqLAzJkzJ3bs2BHPPPNMzJo1K1avXh3r16+PhQsXZs0HQJ2q6CT/SKjkBBHUKyf565eT/F8s7SQ/AJRLYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSVHQ/GKA8tXjBRBfgLE8trlMt/j6VwxEMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASDFRdUeABgdpVKp2iOcolAoVHuEulCv6+QIBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASFFRYAYHB2PVqlUxffr0aGpqiiuvvDJWr15dk1dpBaC6Krpc/7p162Ljxo2xdevWmDlzZrz66quxePHiaG1tjWXLlmXNCEAdqigwr7zySnz729+OBQsWRETEtGnT4plnnokDBw6kDAdA/aroI7J58+ZFT09PHD58OCIi3njjjdi7d2/ccccdKcMBUL8qOoLp6uqKYrEYM2bMiIaGhhgcHIw1a9bEwoULz7jPwMBADAwMnPzfxWLx7KcFoG5UdATz/PPPx7Zt22L79u1x8ODB2Lp1a/ziF7+IrVu3nnGf7u7uaG1tPfloa2s756EBqH2FUgV/AtbW1hZdXV3R2dl5ctujjz4aTz/9dLz99tun3ed0RzBtbW3R398fLS0t5zA6UO/q9V7zRFnv4RV9RPbJJ5/EmDHDD3oaGhpiaGjojPs0NjZGY2NjJS8DwHmgosDceeedsWbNmpg6dWrMnDkzXn/99Xj88cfj3nvvzZoPgDpV0Udkx48fj1WrVsWOHTvi2LFjMXny5PjhD38YDz74YIwdO7asn1EsFqO1tdVHZICPyOpYOe/hFQVmJAgM8DmBqV/lvIe7FhkAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIEVFl+sHGEmjfK3dsrgA58hxBANACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQ4qLRfsFSqRQREcVicbRfGoAR8vl7+RcZ9cAcP348IiLa2tpG+6UBGCHHjx+P1tbWL3xOoVROhkbQ0NBQvP/++9Hc3ByFQuGsf06xWIy2trbo6+uLlpaWEZzw/GKdymOdymOdynM+r1OpVIrjx4/H5MmTY8yYLz7LMupHMGPGjIkpU6aM2M9raWk57/4PzGCdymOdymOdynO+rtOXHbl8zkl+AFIIDAAp6jYwjY2N8dBDD0VjY2O1R6lp1qk81qk81qk81uk/Rv0kPwAXhro9ggGgtgkMACkEBoAUAgNACoEBIIXAAJBCYABIITAApPgf8x05gPAsLjIAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Pentadiagonal (bandwidth 5)"
      ],
      "metadata": {
        "id": "hPfI_RzmPrNG"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# prompt: give me a pentadiagonal matrix of dimension 10 and view its sparsity pattern with spy\n",
        "\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Create a pentadiagonal matrix of dimension 10\n",
        "n = 10\n",
        "A = np.zeros((n, n))\n",
        "np.fill_diagonal(A, 2)  # Main diagonal\n",
        "np.fill_diagonal(A[1:], -1)  # Upper diagonal\n",
        "np.fill_diagonal(A[:, 1:], -1)  # Lower diagonal\n",
        "np.fill_diagonal(A[2:], -1)  # Second upper diagonal\n",
        "np.fill_diagonal(A[:, 2:], -1)  # Second lower diagonal\n",
        "\n",
        "print(A)\n",
        "\n",
        "# View the sparsity pattern\n",
        "plt.spy(A)\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 609
        },
        "id": "lyo0fqEpPttO",
        "outputId": "8500298c-a733-4cba-8a31-0ce7d43f2656"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[[ 2. -1. -1.  0.  0.  0.  0.  0.  0.  0.]\n",
            " [-1.  2. -1. -1.  0.  0.  0.  0.  0.  0.]\n",
            " [-1. -1.  2. -1. -1.  0.  0.  0.  0.  0.]\n",
            " [ 0. -1. -1.  2. -1. -1.  0.  0.  0.  0.]\n",
            " [ 0.  0. -1. -1.  2. -1. -1.  0.  0.  0.]\n",
            " [ 0.  0.  0. -1. -1.  2. -1. -1.  0.  0.]\n",
            " [ 0.  0.  0.  0. -1. -1.  2. -1. -1.  0.]\n",
            " [ 0.  0.  0.  0.  0. -1. -1.  2. -1. -1.]\n",
            " [ 0.  0.  0.  0.  0.  0. -1. -1.  2. -1.]\n",
            " [ 0.  0.  0.  0.  0.  0.  0. -1. -1.  2.]]\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAGiCAYAAADa2tCeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAUM0lEQVR4nO3df2hd9d3A8c9tXNMgSdC6VEtT7WRQbf1RTVu04BgWRTqZY7gfVHD1P4m2tTAWN2qVWmPHJgXrOiujFLT+gNHphG6UDtt1WvpLRdlsN4QtU9oqSG6tECW5zx977J48bfXeNJ+ee9PXCy7i8Z7cD8fkvvnek5xTqlQqlQCAUTau6AEAGJsEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFA0bmCeeeCIuueSSmDBhQsydOzd2795d9Eh1pbe3N2bPnh2tra3R0dERt912Wxw4cKDoserao48+GqVSKZYuXVr0KHXnvffeizvuuCMmTpwYLS0tccUVV8TevXuLHquuDA4OxvLly2PatGnR0tISl156aaxcuTLO5qtxNWRgnn/++Vi2bFmsWLEi9u/fH1dddVXcfPPNceTIkaJHqxvbt2+P7u7u2LVrV2zdujU+++yzuOmmm+LYsWNFj1aX9uzZE08++WRceeWVRY9Sdz766KOYN29efOUrX4ktW7bEX//61/jlL38Z5513XtGj1ZXVq1fHunXrYu3atfG3v/0tVq9eHT//+c/j8ccfL3q0wpQa8WKXc+fOjdmzZ8fatWsjImJoaCg6Ozvj3nvvjZ6enoKnq08ffPBBdHR0xPbt2+OGG24oepy68vHHH8c111wTv/rVr+Lhhx+Oq6++OtasWVP0WHWjp6cn/vKXv8Sf//znokepa9/61rdi0qRJ8Zvf/Ob4tu9+97vR0tISTz/9dIGTFafhVjCffvpp7Nu3L+bPn39827hx42L+/Pnx2muvFThZfevv74+IiPPPP7/gSepPd3d3LFiwYNj3FP/10ksvRVdXV9x+++3R0dERs2bNiqeeeqroserO9ddfH9u2bYuDBw9GRMSbb74ZO3fujFtuuaXgyYpzTtED1OrDDz+MwcHBmDRp0rDtkyZNinfeeaegqerb0NBQLF26NObNmxczZ84sepy68txzz8X+/ftjz549RY9St959991Yt25dLFu2LH7605/Gnj17YvHixTF+/Pi48847ix6vbvT09ES5XI7p06dHU1NTDA4OxqpVq2LhwoVFj1aYhgsMtevu7o633347du7cWfQodaWvry+WLFkSW7dujQkTJhQ9Tt0aGhqKrq6ueOSRRyIiYtasWfH222/Hr3/9a4H5P1544YV45plnYtOmTTFjxox44403YunSpTF58uSz9jg1XGAuuOCCaGpqisOHDw/bfvjw4bjwwgsLmqp+3XPPPfHyyy/Hjh07YsqUKUWPU1f27dsXR44ciWuuueb4tsHBwdixY0esXbs2BgYGoqmpqcAJ68NFF10Ul19++bBtl112Wfz2t78taKL69OMf/zh6enriBz/4QUREXHHFFfHPf/4zent7z9rANNw5mPHjx8e1114b27ZtO75taGgotm3bFtddd12Bk9WXSqUS99xzT2zevDn+9Kc/xbRp04oeqe7ceOON8dZbb8Ubb7xx/NHV1RULFy6MN954Q1z+17x58074FfeDBw/GxRdfXNBE9emTTz6JceOGv6U2NTXF0NBQQRMVr+FWMBERy5YtizvvvDO6urpizpw5sWbNmjh27FgsWrSo6NHqRnd3d2zatClefPHFaG1tjUOHDkVERHt7e7S0tBQ8XX1obW094ZzUueeeGxMnTnSu6v+477774vrrr49HHnkkvve978Xu3btj/fr1sX79+qJHqyu33nprrFq1KqZOnRozZsyI119/PR577LG46667ih6tOJUG9fjjj1emTp1aGT9+fGXOnDmVXbt2FT1SXYmIkz42bNhQ9Gh17Rvf+EZlyZIlRY9Rd37/+99XZs6cWWlubq5Mnz69sn79+qJHqjvlcrmyZMmSytSpUysTJkyofO1rX6v87Gc/qwwMDBQ9WmEa8u9gAKh/DXcOBoDGIDAApBAYAFIIDAApBAaAFAIDQAqBASBFwwZmYGAgHnzwwRgYGCh6lLrmOFXHcaqO41Qdx+k/GvYPLcvlcrS3t0d/f3+0tbUVPU7dcpyq4zhVx3GqjuP0Hw27ggGgvgkMACnO+NWUh4aG4v3334/W1tYolUoj/jrlcnnYPzk5x6k6jlN1HKfqjOXjVKlU4ujRozF58uQTbk/w/53xczD//ve/o7Oz80y+JACjrK+v70tvYnjGVzCtra1n+iUZRf39/UWPABSoXC5HZ2dnVe/lZzwwp/OxGMU7m38jBvivat7LneQHIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASDGiwDzxxBNxySWXxIQJE2Lu3Lmxe/fu0Z4LgAZXc2Cef/75WLZsWaxYsSL2798fV111Vdx8881x5MiRjPkAaFA13w9m7ty5MXv27Fi7dm1E/OcGYp2dnXHvvfdGT0/Pl+7/+b2qaUxn+PZBQJ35/D28v7//S6+uXtMK5tNPP419+/bF/Pnz//sFxo2L+fPnx2uvvXbSfQYGBqJcLg97ADD21RSYDz/8MAYHB2PSpEnDtk+aNCkOHTp00n16e3ujvb39+MPdLAHODum/RXb//fdHf3//8UdfX1/2SwJQB2q6o+UFF1wQTU1Ncfjw4WHbDx8+HBdeeOFJ92lubo7m5uaRTwhAQ6ppBTN+/Pi49tprY9u2bce3DQ0NxbZt2+K6664b9eEAaFw1rWAiIpYtWxZ33nlndHV1xZw5c2LNmjVx7NixWLRoUcZ8ADSomgPz/e9/Pz744IN44IEH4tChQ3H11VfHH/7whxNO/ANwdqv572BOl7+DaWz+DgbObml/BwMA1RIYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKWq+2CVnt1KpVPQIJ3B9NKhPVjAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBTnFD0AnK5SqVT0CCeoVCpFjwCFs4IBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKWoKTG9vb8yePTtaW1ujo6Mjbrvttjhw4EDWbAA0sJoCs3379uju7o5du3bF1q1b47PPPoubbropjh07ljUfAA2qVDmNOyN98MEH0dHREdu3b48bbrihqn3K5XK0t7eP9CWhIbjhGGPV5+/h/f390dbW9oXPPa07Wvb390dExPnnn3/K5wwMDMTAwMCw4QAY+0Z8kn9oaCiWLl0a8+bNi5kzZ57yeb29vdHe3n780dnZOdKXBKCBjPgjsrvvvju2bNkSO3fujClTppzyeSdbwYgMY52PyBir0j8iu+eee+Lll1+OHTt2fGFcIiKam5ujubl5JC8DQAOrKTCVSiXuvffe2Lx5c7zyyisxbdq0rLkAaHA1Baa7uzs2bdoUL774YrS2tsahQ4ciIqK9vT1aWlpSBgSgMdV0DqZUKp10+4YNG+JHP/pRVV/DrylzNnAOhrEq7RyMHxoAquVaZACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApTuuWyaejmgulnUmnupAnjEQ9fj+5liBnmhUMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASDFOUUPUC8qlUrRI5ygVCoVPQJjSD1+P9Xjzx2jxwoGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApDitwDz66KNRKpVi6dKlozQOAGPFiAOzZ8+eePLJJ+PKK68czXkAGCNGFJiPP/44Fi5cGE899VScd955oz0TAGPAiALT3d0dCxYsiPnz53/pcwcGBqJcLg97ADD21XzL5Oeeey72798fe/bsqer5vb298dBDD9U8GACNraYVTF9fXyxZsiSeeeaZmDBhQlX73H///dHf33/80dfXN6JBAWgspUqlUqn2yb/73e/iO9/5TjQ1NR3fNjg4GKVSKcaNGxcDAwPD/tvJlMvlaG9vj/7+/mhraxv55GeBUqlU9AiQqoa3H+pELe/hNX1EduONN8Zbb701bNuiRYti+vTp8ZOf/ORL4wLA2aOmwLS2tsbMmTOHbTv33HNj4sSJJ2wH4OzmL/kBSFHzb5H9f6+88soojAHAWGMFA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJDitK9FRp56vFeGe9Qwmurx+6kef+4alRUMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASDFOUUPQGOpVCpFj3CCUqlU9AiMIfX4/VSPP3fVsIIBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABS1ByY9957L+64446YOHFitLS0xBVXXBF79+7NmA2ABlbT5fo/+uijmDdvXnzzm9+MLVu2xFe/+tX4+9//Huedd17WfAA0qJoCs3r16ujs7IwNGzYc3zZt2rRRHwqAxlfTR2QvvfRSdHV1xe233x4dHR0xa9aseOqpp7JmA6CB1RSYd999N9atWxdf//rX449//GPcfffdsXjx4ti4ceMp9xkYGIhyuTzsAcDYV6rUcC/O8ePHR1dXV7z66qvHty1evDj27NkTr7322kn3efDBB+Ohhx46YXt/f3+0tbWNYGQYrh5vcQujqZ5umVwul6O9vb2q9/CaVjAXXXRRXH755cO2XXbZZfGvf/3rlPvcf//90d/ff/zR19dXy0sC0KBqOsk/b968OHDgwLBtBw8ejIsvvviU+zQ3N0dzc/PIpgOgYdW0grnvvvti165d8cgjj8Q//vGP2LRpU6xfvz66u7uz5gOgQdUUmNmzZ8fmzZvj2WefjZkzZ8bKlStjzZo1sXDhwqz5AGhQNZ3kHw21nCCCajjJz1h3VpzkB4BqCQwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNAipruBwP1qJ4uBPg5F+BkNDXq95MVDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEgxTlFDwBjUaVSKXqEE5RKpaJH4CxjBQNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKSoKTCDg4OxfPnymDZtWrS0tMSll14aK1eurMsrxwJQrJou17969epYt25dbNy4MWbMmBF79+6NRYsWRXt7eyxevDhrRgAaUE2BefXVV+Pb3/52LFiwICIiLrnkknj22Wdj9+7dKcMB0Lhq+ojs+uuvj23btsXBgwcjIuLNN9+MnTt3xi233JIyHACNq6YVTE9PT5TL5Zg+fXo0NTXF4OBgrFq1KhYuXHjKfQYGBmJgYOD4v5fL5ZFPC0DDqGkF88ILL8QzzzwTmzZtiv3798fGjRvjF7/4RWzcuPGU+/T29kZ7e/vxR2dn52kPDUD9K1Vq+BWwzs7O6Onpie7u7uPbHn744Xj66afjnXfeOek+J1vBdHZ2Rn9/f7S1tZ3G6EAtSqVS0SMwhlTzHl7TR2SffPJJjBs3fNHT1NQUQ0NDp9ynubk5mpuba3kZAMaAmgJz6623xqpVq2Lq1KkxY8aMeP311+Oxxx6Lu+66K2s+ABpUTR+RHT16NJYvXx6bN2+OI0eOxOTJk+OHP/xhPPDAAzF+/Piqvka5XI729nYfkcEZ5iMyRlM17+E1BWY0CAwUQ2AYTdW8h7sWGQApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEgRU2X6wca1xm+rm1VXIBzbLOCASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEhxzpl+wUqlEhER5XL5TL80AKPk8/fyL3LGA3P06NGIiOjs7DzTLw3AKDl69Gi0t7d/4XNKlWoyNIqGhobi/fffj9bW1iiVSiP+OuVyOTo7O6Ovry/a2tpGccKxxXGqjuNUHcepOmP5OFUqlTh69GhMnjw5xo374rMsZ3wFM27cuJgyZcqofb22trYx9z8wg+NUHcepOo5TdcbqcfqylcvnnOQHIIXAAJCiYQPT3NwcK1asiObm5qJHqWuOU3Ucp+o4TtVxnP7jjJ/kB+Ds0LArGADqm8AAkEJgAEghMACkEBgAUggMACkEBoAUAgNAiv8B78wqgESxhh8AAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "One can imagine the benefit of storing such a matrix"
      ],
      "metadata": {
        "id": "jNwEgVVNJ-VD"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### General banded matrix"
      ],
      "metadata": {
        "id": "rDw8T5BpPeY0"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Create a pentadiagonal matrix of dimension 10\n",
        "n = 10\n",
        "A = np.zeros((n, n))\n",
        "np.fill_diagonal(A, 2)  # Main diagonal\n",
        "np.fill_diagonal(A[1:], -1)  # Upper diagonal\n",
        "np.fill_diagonal(A[:, 1:], -1)  # Lower diagonal\n",
        "np.fill_diagonal(A[4:], -2)  # Second upper diagonal\n",
        "np.fill_diagonal(A[:, 4:], -1)  # Second lower diagonal\n",
        "\n",
        "print(A)\n",
        "\n",
        "# View the sparsity pattern\n",
        "plt.spy(A)\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 609
        },
        "id": "8vZ1HBTuPhsP",
        "outputId": "55df358d-fa5a-467f-ca0e-b62b7a5c3ddd"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[[ 2. -1.  0.  0. -1.  0.  0.  0.  0.  0.]\n",
            " [-1.  2. -1.  0.  0. -1.  0.  0.  0.  0.]\n",
            " [ 0. -1.  2. -1.  0.  0. -1.  0.  0.  0.]\n",
            " [ 0.  0. -1.  2. -1.  0.  0. -1.  0.  0.]\n",
            " [-2.  0.  0. -1.  2. -1.  0.  0. -1.  0.]\n",
            " [ 0. -2.  0.  0. -1.  2. -1.  0.  0. -1.]\n",
            " [ 0.  0. -2.  0.  0. -1.  2. -1.  0.  0.]\n",
            " [ 0.  0.  0. -2.  0.  0. -1.  2. -1.  0.]\n",
            " [ 0.  0.  0.  0. -2.  0.  0. -1.  2. -1.]\n",
            " [ 0.  0.  0.  0.  0. -2.  0.  0. -1.  2.]]\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAGiCAYAAADa2tCeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAU6ElEQVR4nO3df6jV9f3A8dfxNq93ce/BbFrSNV0MTO33VSmhMZIiXKwx2g8MnMH+iFtqwphtlIXZzbGFUM0yNhHKfsBwtaCNcGRzJVpWFFu5EWx3hVoQ95jBLe493z/2nfv6Vds5dl/38znXxwMOY5/dc8+L9+fsPPmcz/XzqdTr9XoAwAgbV/QAAIxNAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQIqWDcwDDzwQ06dPjwkTJsT8+fNj165dRY9UKn19fTF37tzo7OyMyZMnx7XXXhtvv/120WOV2j333BOVSiVWrFhR9Cil8+6778b1118fkyZNio6OjjjvvPPi5ZdfLnqsUhkaGorbbrstZsyYER0dHXHOOefEmjVr4mS+GldLBuaJJ56IlStXxurVq2PPnj1xwQUXxFVXXRUHDhwoerTS2L59e/T29sbOnTvjueeei08//TSuvPLKOHToUNGjldLu3bvjoYceivPPP7/oUUrnww8/jAULFsQXvvCFePbZZ+PPf/5z/PznP4+JEycWPVqprFu3LjZs2BD3339//OUvf4l169bFT3/607jvvvuKHq0wlVa82OX8+fNj7ty5cf/990dExPDwcHR3d8fNN98cq1atKni6cnr//fdj8uTJsX379rj88suLHqdUPvroo7j44ovjF7/4Rdx1111x4YUXxvr164seqzRWrVoVf/rTn+KPf/xj0aOU2te//vWYMmVK/PKXvzy87Vvf+lZ0dHTEI488UuBkxWm5I5hPPvkkXnnllVi4cOHhbePGjYuFCxfGSy+9VOBk5TYwMBAREaeddlrBk5RPb29vLFq06Ij3FP/x9NNPR09PT1x33XUxefLkuOiii+Lhhx8ueqzSueyyy2Lbtm2xd+/eiIh4/fXXY8eOHXH11VcXPFlxTil6gGZ98MEHMTQ0FFOmTDli+5QpU+Ktt94qaKpyGx4ejhUrVsSCBQtizpw5RY9TKo8//njs2bMndu/eXfQopfXOO+/Ehg0bYuXKlfHjH/84du/eHcuWLYvx48fHkiVLih6vNFatWhW1Wi1mzpwZbW1tMTQ0FGvXro3FixcXPVphWi4wNK+3tzfefPPN2LFjR9GjlEp/f38sX748nnvuuZgwYULR45TW8PBw9PT0xN133x0RERdddFG8+eab8eCDDwrM//Hkk0/Go48+Glu2bInZs2fHa6+9FitWrIipU6eetOvUcoE5/fTTo62tLfbv33/E9v3798cZZ5xR0FTlddNNN8UzzzwTL7zwQpx11llFj1Mqr7zyShw4cCAuvvjiw9uGhobihRdeiPvvvz8GBwejra2twAnL4cwzz4xZs2Ydse3cc8+NX//61wVNVE4//OEPY9WqVfHd7343IiLOO++8+Pvf/x59fX0nbWBa7hzM+PHj45JLLolt27Yd3jY8PBzbtm2LSy+9tMDJyqVer8dNN90UW7dujT/84Q8xY8aMokcqnSuuuCLeeOONeO211w4/enp6YvHixfHaa6+Jy/9asGDBUX/ivnfv3jj77LMLmqicPv744xg37siP1La2thgeHi5oouK13BFMRMTKlStjyZIl0dPTE/PmzYv169fHoUOHYunSpUWPVhq9vb2xZcuWeOqpp6KzszP27dsXERHVajU6OjoKnq4cOjs7jzondeqpp8akSZOcq/o/brnllrjsssvi7rvvjm9/+9uxa9eu2LhxY2zcuLHo0UrlmmuuibVr18a0adNi9uzZ8eqrr8a9994bN9xwQ9GjFafeou677776tGnT6uPHj6/PmzevvnPnzqJHKpWIOOZj06ZNRY9Wal/96lfry5cvL3qM0vntb39bnzNnTr29vb0+c+bM+saNG4seqXRqtVp9+fLl9WnTptUnTJhQ//KXv1z/yU9+Uh8cHCx6tMK05L+DAaD8Wu4cDACtQWAASCEwAKQQGABSCAwAKQQGgBQCA0CKlg3M4OBg3HHHHTE4OFj0KKVmnRpjnRpjnRpjnf6lZf+hZa1Wi2q1GgMDA9HV1VX0OKVlnRpjnRpjnRpjnf6lZY9gACg3gQEgxahfTXl4eDjee++96OzsjEqlcsK/p1arHfGfHJt1aox1aox1asxYXqd6vR4HDx6MqVOnHnV7gv9v1M/B/POf/4zu7u7RfEkARlh/f/9/vYnhqB/BdHZ2jvZLtqyBgYGiR2gJ1Wq16BGOYt8xVtVqteju7m7os3zUA/N5vhY72ZzMf33S6uw7xrpGPsud5AcghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIcUKBeeCBB2L69OkxYcKEmD9/fuzatWuk5wKgxTUdmCeeeCJWrlwZq1evjj179sQFF1wQV111VRw4cCBjPgBaVNOBuffee+MHP/hBLF26NGbNmhUPPvhgfPGLX4xf/epXGfMB0KKaCswnn3wSr7zySixcuPA/v2DcuFi4cGG89NJLx3zO4OBg1Gq1Ix4AjH1NBeaDDz6IoaGhmDJlyhHbp0yZEvv27Tvmc/r6+qJarR5+uJslwMkh/a/Ibr311hgYGDj86O/vz35JAEqgqTtann766dHW1hb79+8/Yvv+/fvjjDPOOOZz2tvbo729/cQnBKAlNXUEM378+Ljkkkti27Zth7cNDw/Htm3b4tJLLx3x4QBoXU0dwURErFy5MpYsWRI9PT0xb968WL9+fRw6dCiWLl2aMR8ALarpwHznO9+J999/P26//fbYt29fXHjhhfG73/3uqBP/AJzcKvV6vT6aL1ir1aJarY7mS7asUd41LatSqRQ9wlHsO8aqf3+GDwwMRFdX12f+rGuRAZBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKRo+mKXjB7X2GpMGWey78ARDABJBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEgxSlFvfDAwEB0dXUV9fJHqVQqRY/QEsq4TvV6vegRjlLGmew7RpsjGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJCiqcD09fXF3Llzo7OzMyZPnhzXXnttvP3221mzAdDCmgrM9u3bo7e3N3bu3BnPPfdcfPrpp3HllVfGoUOHsuYDoEVV6p/jjj/vv/9+TJ48ObZv3x6XX355Q8+p1WpRrVbdcIwR46ZVjSnje9y+az3NfIZ/rjtaDgwMRETEaaeddtyfGRwcjMHBwSOGA2DsO+GT/MPDw7FixYpYsGBBzJkz57g/19fXF9Vq9fCju7v7RF8SgBZywl+R3XjjjfHss8/Gjh074qyzzjruzx3rCKa7u9tXZIwYX7M0pozvcfuu9aR/RXbTTTfFM888Ey+88MJnxiUior29Pdrb20/kZQBoYU0Fpl6vx8033xxbt26N559/PmbMmJE1FwAtrqnA9Pb2xpYtW+Kpp56Kzs7O2LdvX0REVKvV6OjoSBkQgNbU1DmY432Hu2nTpvj+97/f0O/wZ8qMNN/jN6aM73H7rvWknYPxZgCgUa5FBkAKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJDic90yeSwp43XWynhxwjIq4zqV8f1Uxpnsu7HNEQwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIMUpRb1wtVot6qWPqV6vFz3CUco4U6VSKXqEllDGdSrj+6mMM9l3I8cRDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEjxuQJzzz33RKVSiRUrVozQOACMFSccmN27d8dDDz0U559//kjOA8AYcUKB+eijj2Lx4sXx8MMPx8SJE0d6JgDGgBMKTG9vbyxatCgWLlz4X392cHAwarXaEQ8Axr6mb5n8+OOPx549e2L37t0N/XxfX1/ceeedTQ8GQGtr6gimv78/li9fHo8++mhMmDChoefceuutMTAwcPjR399/QoMC0Foq9Xq93ugP/+Y3v4lvfvOb0dbWdnjb0NBQVCqVGDduXAwODh7xvx1LrVaLarV64hMnaWIZTmqVSqXoEThB3uONKeN7vEz77t+f4QMDA9HV1fWZP9vUV2RXXHFFvPHGG0dsW7p0acycOTN+9KMf/de4AHDyaCownZ2dMWfOnCO2nXrqqTFp0qSjtgNwcvMv+QFI0fRfkf1/zz///AiMAcBY4wgGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIMXnvhbZiWrkXgKjyT0gGlPGmcq478qojOtUxvdTGWcq475rhCMYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0CKU4oeoCzq9XrRIxylUqkUPcJRyrhOZZypjPuujMq4TmV8P5VpplqtFtVqtaGfdQQDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkaDow7777blx//fUxadKk6OjoiPPOOy9efvnljNkAaGFNXa7/ww8/jAULFsTXvva1ePbZZ+NLX/pS/PWvf42JEydmzQdAi2oqMOvWrYvu7u7YtGnT4W0zZswY8aEAaH1NfUX29NNPR09PT1x33XUxefLkuOiii+Lhhx/Omg2AFtZUYN55553YsGFDfOUrX4nf//73ceONN8ayZcti8+bNx33O4OBg1Gq1Ix4AjH2VehP34hw/fnz09PTEiy++eHjbsmXLYvfu3fHSSy8d8zl33HFH3HnnnUdtHxgYiK6urhMY+eThdrKtq4z7jsZ4j3+2f98yuZHP8KaOYM4888yYNWvWEdvOPffc+Mc//nHc59x6660xMDBw+NHf39/MSwLQopo6yb9gwYJ4++23j9i2d+/eOPvss4/7nPb29mhvbz+x6QBoWU0dwdxyyy2xc+fOuPvuu+Nvf/tbbNmyJTZu3Bi9vb1Z8wHQopoKzNy5c2Pr1q3x2GOPxZw5c2LNmjWxfv36WLx4cdZ8ALSopk7yj4RmThCd7Mp4otgJ0MaUcd/RGO/xz5Z2kh8AGiUwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACmauh8Mo6uMF90r40Ucy7hOZZypjPuujMq4TmV8PzXCEQwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIMUpRQ9Aa6nX60WPcJRKpVL0CEcp4zqVcaYy7rsyatV1cgQDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkaCowQ0NDcdttt8WMGTOio6MjzjnnnFizZk0pr9IKQLGaulz/unXrYsOGDbF58+aYPXt2vPzyy7F06dKoVquxbNmyrBkBaEFNBebFF1+Mb3zjG7Fo0aKIiJg+fXo89thjsWvXrpThAGhdTX1Fdtlll8W2bdti7969ERHx+uuvx44dO+Lqq69OGQ6A1tXUEcyqVauiVqvFzJkzo62tLYaGhmLt2rWxePHi4z5ncHAwBgcHD//3Wq124tMC0DKaOoJ58skn49FHH40tW7bEnj17YvPmzfGzn/0sNm/efNzn9PX1RbVaPfzo7u7+3EMDUH6VehN/Atbd3R2rVq2K3t7ew9vuuuuueOSRR+Ktt9465nOOdQTT3d0dAwMD0dXV9TlGh38p4/3K/WVlY8q472hMI5/hTX1F9vHHH8e4cUce9LS1tcXw8PBxn9Pe3h7t7e3NvAwAY0BTgbnmmmti7dq1MW3atJg9e3a8+uqrce+998YNN9yQNR8ALaqpr8gOHjwYt912W2zdujUOHDgQU6dOje9973tx++23x/jx4xv6HbVaLarVqq/IGDFl/JrFV2SNKeO+ozGNfIY3FZiRIDCMtDJ+SAlMY8q472hMI5/hrkUGQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEjR1OX6oYzKeGHJMl7EsYzrVMaZyrjvWpUjGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUp4z2C9br9YiIqNVqo/3ScFLz/zlG0r8/yz/LqAfm4MGDERHR3d092i8NJ7VqtVr0CIwhBw8e/K/vqUq9kQyNoOHh4Xjvvfeis7MzKpXKCf+eWq0W3d3d0d/fH11dXSM44dhinRpjnRpjnRozltepXq/HwYMHY+rUqTFu3GefZRn1I5hx48bFWWedNWK/r6ura8ztwAzWqTHWqTHWqTFjdZ0aPRp2kh+AFAIDQIqWDUx7e3usXr062tvbix6l1KxTY6xTY6xTY6zTv4z6SX4ATg4tewQDQLkJDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKT4HwqIjvuxW4sDAAAAAElFTkSuQmCC\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "## Block matricies"
      ],
      "metadata": {
        "id": "S8W1_sPd9bpO"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Block matricies can be partitioned:\n",
        "\n",
        "$$ M=\n",
        "\\begin{bmatrix}\n",
        "2 & -1 & 0 & 1 & 0 & 0 \\\\\n",
        "-1 & 2 & -1 & 0 & 1 & 0 \\\\\n",
        "0 & -1 & 2 & 0 & 0 & 1 \\\\\n",
        "1 & 0 & 0 & 2 & -1 & 0 \\\\\n",
        "0 & 1 & 0 & -1 & 2 & -1 \\\\\n",
        "0 & 0 & 1 & 0 & -1 & 2\n",
        "\\end{bmatrix} = \\begin{bmatrix}\n",
        "A & B \\\\\n",
        "C & D\n",
        "\\end{bmatrix}\n",
        "$$\n",
        "\n",
        "with\n",
        "\n",
        "$$\n",
        "A=D = \\begin{bmatrix}\n",
        "2 & -1 & 0  \\\\\n",
        "-1 & 2 & -1  \\\\\n",
        "0 & -1 & 2 \\\\\n",
        "\\end{bmatrix}\n",
        "$$\n",
        "\n",
        "and\n",
        "\n",
        "$$\n",
        "B=C = \\begin{bmatrix}\n",
        "1 & 0 & 0  \\\\\n",
        "0 & 1 & 0  \\\\\n",
        "0 & 0 & 1 \\\\\n",
        "\\end{bmatrix}\n",
        "$$\n",
        "\n",
        "The inverse can be formed:\n",
        "\n",
        "$$\n",
        "M^{-1} =\n",
        "\\begin{bmatrix}\n",
        "A^{-1} + A^{-1}B(D - CA^{-1}B)^{-1}CA^{-1} & -A^{-1}B(D - CA^{-1}B)^{-1} \\\\\n",
        "-(D - CA^{-1}B)^{-1}CA^{-1} & (D - CA^{-1}B)^{-1}\n",
        "\\end{bmatrix}\n",
        "$$\n",
        "\n",
        "This is a substantial savings since the cost of each block scales exponentially!\n",
        "\n",
        "Sometimes we get lucky and one of the blocks is $0$ or $I$ which is a HUGE benefit!\n",
        "\n",
        "If the blocks are each $n \\times n$, $M$ is $2n \\times 2n$ and takes $O((2n)^3)$ to solve. But if we partition into blocks, the inverse can be formed in 2 separate $O(n^3)$ (plus matrix multiplications)."
      ],
      "metadata": {
        "id": "W8AtA9bOlWUs"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "This situation arises commonly in multiphysics applications, with two fields (e.g.: T and c). The blocks are then defined corresponding to the unknown vectors for each field.\n",
        "\n",
        "Certaint types of physics is solved more efficiently with certain solvers. This decomposition, sometimes called *operator splitting* or *segregated solving* allows each variable (block) to be treated optimally and inverse of the full system assembled.  \n",
        "\n",
        "In contrast, solving the matrix as a whole is called a *monolithic* solution."
      ],
      "metadata": {
        "id": "ZbLvSENLxIQM"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Block diagonal"
      ],
      "metadata": {
        "id": "O4Sx4JVU9ggE"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Taking $B=C=0$ in the above we can see the inverse of a *block diagonal* matrix is simply the diagonal assembly of the inverse of its blocks:\n",
        "\n",
        "$$\n",
        "\\begin{bmatrix}\n",
        "A & 0 \\\\\n",
        "0 & D\n",
        "\\end{bmatrix}^{-1} = \\begin{bmatrix}\n",
        "A^{-1} & 0 \\\\\n",
        "0 & D^{-1}\n",
        "\\end{bmatrix}\n",
        "$$"
      ],
      "metadata": {
        "id": "EI8QL5IYlxv9"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Matrix decompositions"
      ],
      "metadata": {
        "id": "nc0PKYF-77e6"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "We can now summarize some decompositions:\n",
        "\n",
        "| Matrix equation | Operation count |\n",
        "|---|---|\n",
        "| $A = LU$ |  $n^3/3$ |\n",
        "| $S = L D L^T$ | $n^3/6$|\n",
        "| $SPD = L L^T$|  $n^3/6$ + n|\n",
        "\n",
        "The later two are called Choleski decomposition and take advantage of matrix symmetry for efficiency.\n",
        "\n",
        "Since in general\n",
        "\n",
        "$A^{-1} = [LU]^{-1} = U^{-1} L^{-1}$\n",
        "\n",
        "it would be tempting to calculate and store $A^{-1}$ instead of $L$ and $U$, but this is not preferred due to the impact on sparsity, in addition to round-off error $\\propto |A|^{-1}$."
      ],
      "metadata": {
        "id": "aeonjmEk8o9F"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "##Decompositions of certain matricies"
      ],
      "metadata": {
        "id": "aW7I_9NK9nE3"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Inverse of a triangular matrix is triangular"
      ],
      "metadata": {
        "id": "6mzIifj_6uki"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# prompt: show the inverse of an upper triangular matrix with numpy\n",
        "\n",
        "import numpy as np\n",
        "# Create an upper triangular matrix\n",
        "A = np.array([[1, 2, 3],\n",
        "              [0, 4, 5],\n",
        "              [0, 0, 6]])\n",
        "\n",
        "# Calculate the inverse\n",
        "A_inv = np.linalg.inv(A)\n",
        "\n",
        "print(\"Upper triangular matrix:\\n\", A)\n",
        "print(\"\\nInverse of the upper triangular matrix:\\n\", A_inv)\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "i8TNGXEp6Q0f",
        "outputId": "cba9058e-5717-42de-e28b-795cddc351d7"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Upper triangular matrix:\n",
            " [[1 2 3]\n",
            " [0 4 5]\n",
            " [0 0 6]]\n",
            "\n",
            "Inverse of the upper triangular matrix:\n",
            " [[ 1.         -0.5        -0.08333333]\n",
            " [ 0.          0.25       -0.20833333]\n",
            " [ 0.          0.          0.16666667]]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### LU decomposition of a banded matrix is banded-triangular"
      ],
      "metadata": {
        "id": "pJ-2nwTN73gZ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# prompt: Show that the choleski decomposition of a banded matrix is triangular-banded\n",
        "\n",
        "import numpy as np\n",
        "from scipy.linalg import cholesky\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Create a pentadiagonal matrix of dimension 10\n",
        "n = 10\n",
        "A = np.zeros((n, n))\n",
        "np.fill_diagonal(A, 2)  # Main diagonal\n",
        "np.fill_diagonal(A[1:], 1)  # Upper diagonal\n",
        "np.fill_diagonal(A[:, 1:], 1)  # Lower diagonal\n",
        "np.fill_diagonal(A[2:], .1)  # Second upper diagonal\n",
        "np.fill_diagonal(A[:, 2:], .1)  # Second lower diagonal\n",
        "\n",
        "# Calculate the Cholesky decomposition\n",
        "L = cholesky(A, lower=True)\n",
        "\n",
        "print(\"Original matrix:\\n\", A)\n",
        "print(\"\\nCholesky decomposition (L):\\n\", np.around(L, decimals=2))\n",
        "\n",
        "# View the sparsity pattern\n",
        "plt.spy(L)\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 834
        },
        "id": "fZiwnljO8CLN",
        "outputId": "09e9a04c-7b93-45c4-a989-bed37c5d63ff"
      },
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Original matrix:\n",
            " [[2.  1.  0.1 0.  0.  0.  0.  0.  0.  0. ]\n",
            " [1.  2.  1.  0.1 0.  0.  0.  0.  0.  0. ]\n",
            " [0.1 1.  2.  1.  0.1 0.  0.  0.  0.  0. ]\n",
            " [0.  0.1 1.  2.  1.  0.1 0.  0.  0.  0. ]\n",
            " [0.  0.  0.1 1.  2.  1.  0.1 0.  0.  0. ]\n",
            " [0.  0.  0.  0.1 1.  2.  1.  0.1 0.  0. ]\n",
            " [0.  0.  0.  0.  0.1 1.  2.  1.  0.1 0. ]\n",
            " [0.  0.  0.  0.  0.  0.1 1.  2.  1.  0.1]\n",
            " [0.  0.  0.  0.  0.  0.  0.1 1.  2.  1. ]\n",
            " [0.  0.  0.  0.  0.  0.  0.  0.1 1.  2. ]]\n",
            "\n",
            "Cholesky decomposition (L):\n",
            " [[1.41 0.   0.   0.   0.   0.   0.   0.   0.   0.  ]\n",
            " [0.71 1.22 0.   0.   0.   0.   0.   0.   0.   0.  ]\n",
            " [0.07 0.78 1.18 0.   0.   0.   0.   0.   0.   0.  ]\n",
            " [0.   0.08 0.79 1.17 0.   0.   0.   0.   0.   0.  ]\n",
            " [0.   0.   0.08 0.8  1.16 0.   0.   0.   0.   0.  ]\n",
            " [0.   0.   0.   0.09 0.8  1.16 0.   0.   0.   0.  ]\n",
            " [0.   0.   0.   0.   0.09 0.8  1.16 0.   0.   0.  ]\n",
            " [0.   0.   0.   0.   0.   0.09 0.8  1.16 0.   0.  ]\n",
            " [0.   0.   0.   0.   0.   0.   0.09 0.8  1.16 0.  ]\n",
            " [0.   0.   0.   0.   0.   0.   0.   0.09 0.8  1.16]]\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAGiCAYAAADa2tCeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAUPklEQVR4nO3df6jV9f3A8de5t3m9xL0Xs2mJ13QxsNTKuiopNEZShIs1RvuBgbP/4vorYey6YRamN8cWQjaXMkRI+wHD1QI3xJHOpfgro9jSjWC7K/wRxD1mcIt7z/ePfXPf+1XrHPXlOUcfDziMPp3PPS/e3Z0nn/M59/MplEqlUgDARdZQ7QEAuDwJDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKeo2MM8++2yMHTs2hg4dGtOmTYu9e/dWe6Sa0t3dHVOmTImWlpYYMWJEPPDAA3H48OFqj1XTnnrqqSgUCrFo0aJqj1Jz3n///XjooYdi+PDh0dzcHJMmTYr9+/dXe6ya0t/fH0uXLo1x48ZFc3Nz3HjjjbF8+fK4kq/GVZeBeemll2Lx4sWxbNmyOHjwYNx6661x7733xvHjx6s9Ws3YsWNHdHZ2xp49e2Lbtm3x2WefxT333BOnTp2q9mg1ad++ffHcc8/FLbfcUu1Ras5HH30UM2bMiK985SuxdevW+Otf/xq//OUvY9iwYdUeraasWrUq1q5dG2vWrIm//e1vsWrVqvj5z38ezzzzTLVHq5pCPV7sctq0aTFlypRYs2ZNREQMDAxEe3t7zJ8/P7q6uqo8XW06ceJEjBgxInbs2BF33XVXtcepKR9//HHcfvvt8atf/SqefPLJuO2222L16tXVHqtmdHV1xV/+8pf485//XO1Ratq3vvWtGDlyZPzmN785ve273/1uNDc3x/PPP1/Fyaqn7o5gPv300zhw4EDMnDnz9LaGhoaYOXNm7N69u4qT1bbe3t6IiLjmmmuqPEnt6ezsjFmzZg36neK/Xn311ejo6IgHH3wwRowYEZMnT47169dXe6yaM3369Ni+fXscOXIkIiLeeuut2LVrV9x3331Vnqx6rqr2AJX68MMPo7+/P0aOHDlo+8iRI+Pdd9+t0lS1bWBgIBYtWhQzZsyIiRMnVnucmvLiiy/GwYMHY9++fdUepWa99957sXbt2li8eHH89Kc/jX379sWCBQtiyJAhMWfOnGqPVzO6urqiWCzG+PHjo7GxMfr7+2PFihUxe/bsao9WNXUXGCrX2dkZ77zzTuzatavao9SUnp6eWLhwYWzbti2GDh1a7XFq1sDAQHR0dMTKlSsjImLy5MnxzjvvxK9//WuB+T9efvnl2LRpU2zevDkmTJgQhw4dikWLFsWoUaOu2HWqu8Bce+210djYGMeOHRu0/dixY3HddddVaaraNW/evHjttddi586dMXr06GqPU1MOHDgQx48fj9tvv/30tv7+/ti5c2esWbMm+vr6orGxsYoT1obrr78+br755kHbbrrppvjtb39bpYlq049//OPo6uqKH/zgBxERMWnSpPjnP/8Z3d3dV2xg6u4czJAhQ+KOO+6I7du3n942MDAQ27dvjzvvvLOKk9WWUqkU8+bNiy1btsSf/vSnGDduXLVHqjl33313vP3223Ho0KHTj46Ojpg9e3YcOnRIXP7XjBkzzviK+5EjR+KGG26o0kS16ZNPPomGhsFvqY2NjTEwMFCliaqv7o5gIiIWL14cc+bMiY6Ojpg6dWqsXr06Tp06FXPnzq32aDWjs7MzNm/eHK+88kq0tLTE0aNHIyKira0tmpubqzxdbWhpaTnjnNTVV18dw4cPd67q/3j00Udj+vTpsXLlyvje974Xe/fujXXr1sW6deuqPVpNuf/++2PFihUxZsyYmDBhQrz55pvx9NNPx8MPP1zt0aqnVKeeeeaZ0pgxY0pDhgwpTZ06tbRnz55qj1RTIuKsjw0bNlR7tJr2jW98o7Rw4cJqj1Fzfv/735cmTpxYampqKo0fP760bt26ao9Uc4rFYmnhwoWlMWPGlIYOHVr62te+VvrZz35W6uvrq/ZoVVOXfwcDQO2ru3MwANQHgQEghcAAkEJgAEghMACkEBgAUggMACnqNjB9fX3x+OOPR19fX7VHqWnWqTzWqTzWqTzW6T/q9g8ti8VitLW1RW9vb7S2tlZ7nJplncpjncpjncpjnf6jbo9gAKhtAgNAikt+NeWBgYH44IMPoqWlJQqFwnn/nGKxOOh/OTvrVB7rVB7rVJ7LeZ1KpVKcPHkyRo0adcbtCf6/S34O5t///ne0t7dfypcE4CLr6en50psYXvIjmJaWlkv9kmXp7e2t9ggANa9YLEZ7e3tZ7+WXPDAX8rFYpiv5mx4AlSrnvdxJfgBSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAU5xWYZ599NsaOHRtDhw6NadOmxd69ey/2XADUuYoD89JLL8XixYtj2bJlcfDgwbj11lvj3nvvjePHj2fMB0Cdqvh+MNOmTYspU6bEmjVrIuI/NxBrb2+P+fPnR1dX15fu//m9qmvNJb4tDkBd+vw9vLe390uvQl/REcynn34aBw4ciJkzZ/73BzQ0xMyZM2P37t1n3aevry+KxeKgBwCXv4oC8+GHH0Z/f3+MHDly0PaRI0fG0aNHz7pPd3d3tLW1nX64myXAlSH9W2RLliyJ3t7e04+enp7slwSgBlR0R8trr702Ghsb49ixY4O2Hzt2LK677rqz7tPU1BRNTU3nPyEAdamiI5ghQ4bEHXfcEdu3bz+9bWBgILZv3x533nnnRR8OgPpV0RFMRMTixYtjzpw50dHREVOnTo3Vq1fHqVOnYu7cuRnzAVCnKg7M97///Thx4kQ89thjcfTo0bjtttviD3/4wxkn/gG4slX8dzAXyt/BANSvtL+DAYByCQwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUFV/s8nJVKBSqPcIZXB8NqGeOYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKa6q9gCcW6FQqPYIZyiVStUeAagTjmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNAiooC093dHVOmTImWlpYYMWJEPPDAA3H48OGs2QCoYxUFZseOHdHZ2Rl79uyJbdu2xWeffRb33HNPnDp1Kms+AOpUoXQBd5A6ceJEjBgxInbs2BF33XVXWfsUi8Voa2s735ekytxwDK5sn7+H9/b2Rmtr6xc+94LuaNnb2xsREddcc805n9PX1xd9fX2DhgPg8nfeJ/kHBgZi0aJFMWPGjJg4ceI5n9fd3R1tbW2nH+3t7ef7kgDUkfP+iOyRRx6JrVu3xq5du2L06NHnfN7ZjmBEpn75iAyubOkfkc2bNy9ee+212Llz5xfGJSKiqakpmpqazudlAKhjFQWmVCrF/PnzY8uWLfH666/HuHHjsuYCoM5VFJjOzs7YvHlzvPLKK9HS0hJHjx6NiIi2trZobm5OGRCA+lTROZhCoXDW7Rs2bIgf/ehHZf0MX1Oub87BwJUt7RyMNxcAyuVaZACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApLuiWyReinAulXUrnupAng9XiOrlGHtQmRzAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBRXVXuAWlEqlao9whkKhUK1R6gLtbhOtfj7BJeaIxgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQ4oIC89RTT0WhUIhFixZdpHEAuFycd2D27dsXzz33XNxyyy0Xcx4ALhPnFZiPP/44Zs+eHevXr49hw4Zd7JkAuAycV2A6Oztj1qxZMXPmzC99bl9fXxSLxUEPAC5/Fd8y+cUXX4yDBw/Gvn37ynp+d3d3PPHEExUPBkB9q+gIpqenJxYuXBibNm2KoUOHlrXPkiVLore39/Sjp6fnvAYFoL4USqVSqdwn/+53v4vvfOc70djYeHpbf39/FAqFaGhoiL6+vkH/7myKxWK0tbVFb29vtLa2nv/kV4BCoVDtEThPFfzfCupKJe/hFX1Edvfdd8fbb789aNvcuXNj/Pjx8ZOf/ORL4wLAlaOiwLS0tMTEiRMHbbv66qtj+PDhZ2wH4MrmL/kBSFHxt8j+v9dff/0ijAHA5cYRDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0CKC74WGXlq8Z4i7lFTnlpcp1r8feLy5ggGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJDiqmoPQH0plUrVHuEMhUKh2iPUhVpcp1r8feLicQQDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkqDgw77//fjz00EMxfPjwaG5ujkmTJsX+/fszZgOgjlV0uf6PPvooZsyYEd/85jdj69at8dWvfjX+/ve/x7Bhw7LmA6BOVRSYVatWRXt7e2zYsOH0tnHjxl30oQCofxV9RPbqq69GR0dHPPjggzFixIiYPHlyrF+/Pms2AOpYRYF57733Yu3atfH1r389/vjHP8YjjzwSCxYsiI0bN55zn76+vigWi4MeAFz+CqUK7lk6ZMiQ6OjoiDfeeOP0tgULFsS+ffti9+7dZ93n8ccfjyeeeOKM7b29vdHa2noeI8NgtXgrYMrjlsn1p1gsRltbW1nv4RUdwVx//fVx8803D9p20003xb/+9a9z7rNkyZLo7e09/ejp6ankJQGoUxWd5J8xY0YcPnx40LYjR47EDTfccM59mpqaoqmp6fymA6BuVXQE8+ijj8aePXti5cqV8Y9//CM2b94c69ati87Ozqz5AKhTFQVmypQpsWXLlnjhhRdi4sSJsXz58li9enXMnj07az4A6lRFJ/kvhkpOEEE5nOSvX07y15+0k/wAUC6BASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEhR0f1goBbV4gUTXYCzPLW4TrX4+1SvHMEAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFJcVe0B4HJUKpWqPcIZCoVCtUeoC7W4TrX4+1QORzAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0CKigLT398fS5cujXHjxkVzc3PceOONsXz58rq90icAeSq6XP+qVati7dq1sXHjxpgwYULs378/5s6dG21tbbFgwYKsGQGoQxUF5o033ohvf/vbMWvWrIiIGDt2bLzwwguxd+/elOEAqF8VfUQ2ffr02L59exw5ciQiIt56663YtWtX3HfffSnDAVC/KjqC6erqimKxGOPHj4/Gxsbo7++PFStWxOzZs8+5T19fX/T19Z3+52KxeP7TAlA3KjqCefnll2PTpk2xefPmOHjwYGzcuDF+8YtfxMaNG8+5T3d3d7S1tZ1+tLe3X/DQANS+QqmCr4C1t7dHV1dXdHZ2nt725JNPxvPPPx/vvvvuWfc52xFMe3t79Pb2Rmtr6wWMDlSiFu81T3lq6Zu6xWIx2traynoPr+gjsk8++SQaGgYf9DQ2NsbAwMA592lqaoqmpqZKXgaAy0BFgbn//vtjxYoVMWbMmJgwYUK8+eab8fTTT8fDDz+cNR8Adaqij8hOnjwZS5cujS1btsTx48dj1KhR8cMf/jAee+yxGDJkSFk/o5LDK+Di8RFZ/arXj8gqCszFIDBQHQJTv+o1MK5FBkAKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIUdHl+oH6VUsXTPycC3CWp17XyREMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQIqrLvULlkqliIgoFouX+qUBuEg+fy//Ipc8MCdPnoyIiPb29kv90gBcJCdPnoy2trYvfE6hVE6GLqKBgYH44IMPoqWlJQqFwnn/nGKxGO3t7dHT0xOtra0XccLLi3Uqj3Uqj3Uqz+W8TqVSKU6ePBmjRo2KhoYvPstyyY9gGhoaYvTo0Rft57W2tl52/wEzWKfyWKfyWKfyXK7r9GVHLp9zkh+AFAIDQIq6DUxTU1MsW7Ysmpqaqj1KTbNO5bFO5bFO5bFO/3HJT/IDcGWo2yMYAGqbwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0CK/wHuN0KAoVjmIgAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Inverse of a banded matrix is **not** banded (except diagonal)"
      ],
      "metadata": {
        "id": "SOMVhoTm6zqZ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# prompt: show the inverse of a banded matrix in numpy\n",
        "\n",
        "import numpy as np\n",
        "\n",
        "# Create a tridiagonal matrix of dimension 10\n",
        "n = 10\n",
        "A = np.zeros((n, n))\n",
        "np.fill_diagonal(A, 2)  # Main diagonal\n",
        "np.fill_diagonal(A[1:], -1)  # Upper diagonal\n",
        "np.fill_diagonal(A[:, 1:], -1)  # Lower diagonal\n",
        "\n",
        "# Calculate the inverse\n",
        "A_inv = np.linalg.inv(A)\n",
        "\n",
        "print(\"Banded matrix:\\n\")\n",
        "print(A)\n",
        "plt.spy(A)\n",
        "plt.show()\n",
        "print(\"\\nInverse of the banded matrix:\\n\")\n",
        "print(A_inv)\n",
        "plt.spy(A_inv)\n",
        "plt.show()\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "id": "snlLxPZF6S5T",
        "outputId": "0c881edf-cdc9-4fe7-a895-eb0206708aac"
      },
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Banded matrix:\n",
            "\n",
            "[[ 2. -1.  0.  0.  0.  0.  0.  0.  0.  0.]\n",
            " [-1.  2. -1.  0.  0.  0.  0.  0.  0.  0.]\n",
            " [ 0. -1.  2. -1.  0.  0.  0.  0.  0.  0.]\n",
            " [ 0.  0. -1.  2. -1.  0.  0.  0.  0.  0.]\n",
            " [ 0.  0.  0. -1.  2. -1.  0.  0.  0.  0.]\n",
            " [ 0.  0.  0.  0. -1.  2. -1.  0.  0.  0.]\n",
            " [ 0.  0.  0.  0.  0. -1.  2. -1.  0.  0.]\n",
            " [ 0.  0.  0.  0.  0.  0. -1.  2. -1.  0.]\n",
            " [ 0.  0.  0.  0.  0.  0.  0. -1.  2. -1.]\n",
            " [ 0.  0.  0.  0.  0.  0.  0.  0. -1.  2.]]\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAGiCAYAAADa2tCeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAUOUlEQVR4nO3dfWxV9f3A8c+ljtKYtkEcKKEIM0tQwMcikSYui0RjmJnL4h6CCcP/THmSZFndgg9BqCybIRHHhCyERPAhWZjOhC2kizAmhCpiNJuyxWTrNIAmphcxqaa9vz/2k/36A/Re6Kf3Xni9khuz4z29n3xX7zvnnt5zCqVSqRQAMMLGVHsAAM5PAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQIq6DcyTTz4Z06ZNi3HjxsXcuXPjwIED1R6ppnR3d8ecOXOiubk5Jk6cGHfddVe888471R6rpj322GNRKBRixYoV1R6l5rz33ntxzz33xIQJE6KpqSlmz54dr776arXHqimDg4OxatWqmD59ejQ1NcWVV14Zq1evjgv5alx1GZjnnnsuVq5cGQ899FAcPHgwrr322rj99tvj2LFj1R6tZuzevTs6Oztj//79sWvXrvjss8/itttuixMnTlR7tJrU29sbTz31VFxzzTXVHqXmfPTRR9HR0RFf+cpXYufOnfHXv/41fvnLX8b48eOrPVpNWbduXWzcuDE2bNgQf/vb32LdunXx85//PJ544olqj1Y1hXq82OXcuXNjzpw5sWHDhoiIGBoaira2tli6dGl0dXVVebra9MEHH8TEiRNj9+7dccstt1R7nJry8ccfxw033BC/+tWv4tFHH43rrrsu1q9fX+2xakZXV1f85S9/iT//+c/VHqWmfetb34pJkybFb37zm5Pbvvvd70ZTU1M8/fTTVZyseuruCObTTz+N1157LebPn39y25gxY2L+/Pmxb9++Kk5W2/r7+yMi4pJLLqnyJLWns7MzFixYMOx3iv968cUXo729Pe6+++6YOHFiXH/99bF58+Zqj1Vz5s2bFz09PXH48OGIiHjjjTdi7969cccdd1R5suq5qNoDVOrDDz+MwcHBmDRp0rDtkyZNirfffrtKU9W2oaGhWLFiRXR0dMSsWbOqPU5NefbZZ+PgwYPR29tb7VFq1rvvvhsbN26MlStXxk9/+tPo7e2NZcuWxdixY2PRokXVHq9mdHV1RbFYjBkzZkRDQ0MMDg7GmjVrYuHChdUerWrqLjBUrrOzM956663Yu3dvtUepKX19fbF8+fLYtWtXjBs3rtrj1KyhoaFob2+PtWvXRkTE9ddfH2+99Vb8+te/Fpj/4/nnn49t27bF9u3bY+bMmXHo0KFYsWJFTJ48+YJdp7oLzKWXXhoNDQ1x9OjRYduPHj0al112WZWmql1LliyJl156Kfbs2RNTpkyp9jg15bXXXotjx47FDTfccHLb4OBg7NmzJzZs2BADAwPR0NBQxQlrw+WXXx5XX331sG1XXXVV/Pa3v63SRLXpxz/+cXR1dcUPfvCDiIiYPXt2/POf/4zu7u4LNjB1dw5m7NixceONN0ZPT8/JbUNDQ9HT0xM333xzFSerLaVSKZYsWRI7duyIP/3pTzF9+vRqj1Rzbr311njzzTfj0KFDJx/t7e2xcOHCOHTokLj8r46OjlP+xP3w4cNxxRVXVGmi2vTJJ5/EmDHD31IbGhpiaGioShNVX90dwURErFy5MhYtWhTt7e1x0003xfr16+PEiROxePHiao9WMzo7O2P79u3xwgsvRHNzcxw5ciQiIlpbW6OpqanK09WG5ubmU85JXXzxxTFhwgTnqv6P+++/P+bNmxdr166N733ve3HgwIHYtGlTbNq0qdqj1ZQ777wz1qxZE1OnTo2ZM2fG66+/Ho8//njce++91R6tekp16oknnihNnTq1NHbs2NJNN91U2r9/f7VHqikRcdrHli1bqj1aTfvGN75RWr58ebXHqDm///3vS7NmzSo1NjaWZsyYUdq0aVO1R6o5xWKxtHz58tLUqVNL48aNK33ta18r/exnPysNDAxUe7SqqcvvwQBQ++ruHAwA9UFgAEghMACkEBgAUggMACkEBoAUAgNAiroNzMDAQDz88MMxMDBQ7VFqmnUqj3Uqj3Uqj3X6j7r9omWxWIzW1tbo7++PlpaWao9Ts6xTeaxTeaxTeazTf9TtEQwAtU1gAEgx6ldTHhoaivfffz+am5ujUCic9c8pFovD/snpWafyWKfyWKfynM/rVCqV4vjx4zF58uRTbk/w/436OZh///vf0dbWNpovCcAI6+vr+9KbGI76EUxzc/Nov2Td6u/vr/YIAMMUi8Voa2sr67181ANzLh+LXWgu5L8+AWpbOe/lTvIDkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApDirwDz55JMxbdq0GDduXMydOzcOHDgw0nMBUOcqDsxzzz0XK1eujIceeigOHjwY1157bdx+++1x7NixjPkAqFMV3w9m7ty5MWfOnNiwYUNE/OcGYm1tbbF06dLo6ur60v0/v1c1X26Ub9UD8KU+fw/v7+//0iu+V3QE8+mnn8Zrr70W8+fP/+8PGDMm5s+fH/v27TvtPgMDA1EsFoc9ADj/VRSYDz/8MAYHB2PSpEnDtk+aNCmOHDly2n26u7ujtbX15MPdLAEuDOl/RfbAAw9Ef3//yUdfX1/2SwJQAyq6o+Wll14aDQ0NcfTo0WHbjx49Gpdddtlp92lsbIzGxsaznxCAulTREczYsWPjxhtvjJ6enpPbhoaGoqenJ26++eYRHw6A+lXREUxExMqVK2PRokXR3t4eN910U6xfvz5OnDgRixcvzpgPgDpVcWC+//3vxwcffBAPPvhgHDlyJK677rr4wx/+cMqJfwAubBV/D+Zc+R5M+XwPBqg1ad+DAYByCQwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUFV/sktFTKBSqPcIpXB8NKJcjGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNAiouq9cL9/f3R0tJSrZc/RaFQqPYIdaEW16lUKlV7BOA0HMEAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFBUFpru7O+bMmRPNzc0xceLEuOuuu+Kdd97Jmg2AOlZRYHbv3h2dnZ2xf//+2LVrV3z22Wdx2223xYkTJ7LmA6BOFUrncLemDz74ICZOnBi7d++OW265pax9isVitLa2uuEYI8YNx2D0VPIefk53tOzv74+IiEsuueSMzxkYGIiBgYFhwwFw/jvrk/xDQ0OxYsWK6OjoiFmzZp3xed3d3dHa2nry0dbWdrYvCUAdOeuPyO67777YuXNn7N27N6ZMmXLG553uCKatrc1HZIwYH5HB6En/iGzJkiXx0ksvxZ49e74wLhERjY2N0djYeDYvA0AdqygwpVIpli5dGjt27IiXX345pk+fnjUXAHWuosB0dnbG9u3b44UXXojm5uY4cuRIRES0trZGU1NTyoAA1KeKzsGc6TzFli1b4kc/+lFZP8OfKTPSnIOB0ZN2DsZ/yACUy7XIAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFKc0y2Tzye1eJ01F+AsTy2uUy3+PsFocwQDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEhxUbUH4MxKpVK1RzhFoVCo9gh1oRbXqRZ/nzi/OYIBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKc4pMI899lgUCoVYsWLFCI0DwPnirAPT29sbTz31VFxzzTUjOQ8A54mzCszHH38cCxcujM2bN8f48eNHeiYAzgNnFZjOzs5YsGBBzJ8//0ufOzAwEMVicdgDgPNfxbdMfvbZZ+PgwYPR29tb1vO7u7vjkUceqXgwAOpbRUcwfX19sXz58ti2bVuMGzeurH0eeOCB6O/vP/no6+s7q0EBqC+FUqlUKvfJv/vd7+I73/lONDQ0nNw2ODgYhUIhxowZEwMDA8P+3ekUi8VobW2N/v7+aGlpOfvJqYpCoVDtEThLFfynDmdUyXt4RR+R3XrrrfHmm28O27Z48eKYMWNG/OQnP/nSuABw4agoMM3NzTFr1qxh2y6++OKYMGHCKdsBuLD5Jj8AKSr+K7L/7+WXXx6BMQA43ziCASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEhxztci48JSi/cUcY+a8tTiOtXi7xMjxxEMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASDFRdUeAM5VqVSq9ginKBQK1R6hLtTiOtXi71O9cgQDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkqDgw7733Xtxzzz0xYcKEaGpqitmzZ8err76aMRsAdayiy/V/9NFH0dHREd/85jdj586d8dWvfjX+/ve/x/jx47PmA6BOVRSYdevWRVtbW2zZsuXktunTp4/4UADUv4o+InvxxRejvb097r777pg4cWJcf/31sXnz5qzZAKhjFQXm3XffjY0bN8bXv/71+OMf/xj33XdfLFu2LLZu3XrGfQYGBqJYLA57AHD+K5QquD/o2LFjo729PV555ZWT25YtWxa9vb2xb9++0+7z8MMPxyOPPHLK9v7+/mhpaTmLkaH21eKtgCmPWyZ/sWKxGK2trWW9h1d0BHP55ZfH1VdfPWzbVVddFf/617/OuM8DDzwQ/f39Jx99fX2VvCQAdaqik/wdHR3xzjvvDNt2+PDhuOKKK864T2NjYzQ2Np7ddADUrYqOYO6///7Yv39/rF27Nv7xj3/E9u3bY9OmTdHZ2Zk1HwB1qqLAzJkzJ3bs2BHPPPNMzJo1K1avXh3r16+PhQsXZs0HQJ2q6CT/SKjkBBHUKyf565eT/F8s7SQ/AJRLYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSVHQ/GKA8tXjBRBfgLE8trlMt/j6VwxEMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASDFRdUeABgdpVKp2iOcolAoVHuEulCv6+QIBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASFFRYAYHB2PVqlUxffr0aGpqiiuvvDJWr15dk1dpBaC6Krpc/7p162Ljxo2xdevWmDlzZrz66quxePHiaG1tjWXLlmXNCEAdqigwr7zySnz729+OBQsWRETEtGnT4plnnokDBw6kDAdA/aroI7J58+ZFT09PHD58OCIi3njjjdi7d2/ccccdKcMBUL8qOoLp6uqKYrEYM2bMiIaGhhgcHIw1a9bEwoULz7jPwMBADAwMnPzfxWLx7KcFoG5UdATz/PPPx7Zt22L79u1x8ODB2Lp1a/ziF7+IrVu3nnGf7u7uaG1tPfloa2s756EBqH2FUgV/AtbW1hZdXV3R2dl5ctujjz4aTz/9dLz99tun3ed0RzBtbW3R398fLS0t5zA6UO/q9V7zRFnv4RV9RPbJJ5/EmDHDD3oaGhpiaGjojPs0NjZGY2NjJS8DwHmgosDceeedsWbNmpg6dWrMnDkzXn/99Xj88cfj3nvvzZoPgDpV0Udkx48fj1WrVsWOHTvi2LFjMXny5PjhD38YDz74YIwdO7asn1EsFqO1tdVHZICPyOpYOe/hFQVmJAgM8DmBqV/lvIe7FhkAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIEVFl+sHGEmjfK3dsrgA58hxBANACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQ4qLRfsFSqRQREcVicbRfGoAR8vl7+RcZ9cAcP348IiLa2tpG+6UBGCHHjx+P1tbWL3xOoVROhkbQ0NBQvP/++9Hc3ByFQuGsf06xWIy2trbo6+uLlpaWEZzw/GKdymOdymOdynM+r1OpVIrjx4/H5MmTY8yYLz7LMupHMGPGjIkpU6aM2M9raWk57/4PzGCdymOdymOdynO+rtOXHbl8zkl+AFIIDAAp6jYwjY2N8dBDD0VjY2O1R6lp1qk81qk81qk81uk/Rv0kPwAXhro9ggGgtgkMACkEBoAUAgNACoEBIIXAAJBCYABIITAApPgf8x05gPAsLjIAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        },
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "\n",
            "Inverse of the banded matrix:\n",
            "\n",
            "[[0.90909091 0.81818182 0.72727273 0.63636364 0.54545455 0.45454545\n",
            "  0.36363636 0.27272727 0.18181818 0.09090909]\n",
            " [0.81818182 1.63636364 1.45454545 1.27272727 1.09090909 0.90909091\n",
            "  0.72727273 0.54545455 0.36363636 0.18181818]\n",
            " [0.72727273 1.45454545 2.18181818 1.90909091 1.63636364 1.36363636\n",
            "  1.09090909 0.81818182 0.54545455 0.27272727]\n",
            " [0.63636364 1.27272727 1.90909091 2.54545455 2.18181818 1.81818182\n",
            "  1.45454545 1.09090909 0.72727273 0.36363636]\n",
            " [0.54545455 1.09090909 1.63636364 2.18181818 2.72727273 2.27272727\n",
            "  1.81818182 1.36363636 0.90909091 0.45454545]\n",
            " [0.45454545 0.90909091 1.36363636 1.81818182 2.27272727 2.72727273\n",
            "  2.18181818 1.63636364 1.09090909 0.54545455]\n",
            " [0.36363636 0.72727273 1.09090909 1.45454545 1.81818182 2.18181818\n",
            "  2.54545455 1.90909091 1.27272727 0.63636364]\n",
            " [0.27272727 0.54545455 0.81818182 1.09090909 1.36363636 1.63636364\n",
            "  1.90909091 2.18181818 1.45454545 0.72727273]\n",
            " [0.18181818 0.36363636 0.54545455 0.72727273 0.90909091 1.09090909\n",
            "  1.27272727 1.45454545 1.63636364 0.81818182]\n",
            " [0.09090909 0.18181818 0.27272727 0.36363636 0.45454545 0.54545455\n",
            "  0.63636364 0.72727273 0.81818182 0.90909091]]\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAGiCAYAAADa2tCeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAATTElEQVR4nO3df6jV9f3A8de917zexb0XtWmJmhZjllpZ/pgKjZEU4WKN0TYwcPZfXH9NGNNFc8PZzcFCyOZShmuU1mC0foAb4pjOLfFXtmKbBsF2V/gjiHPM2Cnu/Xz/2Jf73f2WdY/5uuccfTzgTdxPn8/5vPgU98nnfO49t6koiiIA4AJrrvUAAFycBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQNG5jHHnssJk2aFCNGjIg5c+bEgQMHaj1SXenu7o5Zs2ZFe3t7jBkzJu6+++44duxYrceqaw8//HA0NTXFypUraz1K3XnzzTfj3nvvjdGjR0dbW1tMnz49Dh06VOux6kpvb288+OCDMXny5Ghra4trr7021q1bF5fyp3E1ZGCeeeaZWLVqVaxduzaOHDkSN954Y9xxxx1x6tSpWo9WN/bs2RNdXV2xf//+2LVrV3zwwQdx++23x9mzZ2s9Wl06ePBgPP7443HDDTfUepS6884778T8+fPjsssui507d8Zf//rX+MlPfhIjR46s9Wh1ZcOGDbF58+bYtGlT/O1vf4sNGzbEj3/843j00UdrPVrNNDXih13OmTMnZs2aFZs2bYqIiL6+vpgwYUIsW7YsVq9eXePp6tPp06djzJgxsWfPnrj11ltrPU5deffdd+Pmm2+On/70p/GjH/0obrrppti4cWOtx6obq1evjj/96U/xxz/+sdaj1LUvf/nLMXbs2Pj5z3/ev+1rX/tatLW1xZNPPlnDyWqn4e5g3n///Th8+HAsWLCgf1tzc3MsWLAgXnrppRpOVt9KpVJERIwaNarGk9Sfrq6uWLhw4YD/p/g/zz//fMycOTPuueeeGDNmTMyYMSO2bt1a67Hqzrx582L37t1x/PjxiIh45ZVXYt++fXHnnXfWeLLaGVbrAar19ttvR29vb4wdO3bA9rFjx8bf//73Gk1V3/r6+mLlypUxf/78mDZtWq3HqStPP/10HDlyJA4ePFjrUerWG2+8EZs3b45Vq1bF9773vTh48GAsX748hg8fHosXL671eHVj9erVUS6XY8qUKdHS0hK9vb2xfv36WLRoUa1Hq5mGCwzV6+rqitdeey327dtX61HqSk9PT6xYsSJ27doVI0aMqPU4dauvry9mzpwZDz30UEREzJgxI1577bX42c9+JjD/5Ve/+lU89dRTsX379pg6dWocPXo0Vq5cGePGjbtkr1PDBeaKK66IlpaWOHny5IDtJ0+ejCuvvLJGU9WvpUuXxosvvhh79+6N8ePH13qcunL48OE4depU3Hzzzf3bent7Y+/evbFp06aoVCrR0tJSwwnrw1VXXRXXX3/9gG3XXXdd/PrXv67RRPXpO9/5TqxevTq++c1vRkTE9OnT4x//+Ed0d3dfsoFpuGcww4cPj1tuuSV2797dv62vry92794dc+fOreFk9aUoili6dGk8++yz8fvf/z4mT55c65Hqzm233RavvvpqHD16tH/NnDkzFi1aFEePHhWX/zV//vwP/Yj78ePH4+qrr67RRPXpvffei+bmgd9SW1paoq+vr0YT1V7D3cFERKxatSoWL14cM2fOjNmzZ8fGjRvj7NmzsWTJklqPVje6urpi+/bt8dxzz0V7e3ucOHEiIiI6Ozujra2txtPVh/b29g89k7r88stj9OjRnlX9l29/+9sxb968eOihh+LrX/96HDhwILZs2RJbtmyp9Wh15a677or169fHxIkTY+rUqfHyyy/HI488Evfdd1+tR6udokE9+uijxcSJE4vhw4cXs2fPLvbv31/rkepKRHzk2rZtW61Hq2tf/OIXixUrVtR6jLrzwgsvFNOmTStaW1uLKVOmFFu2bKn1SHWnXC4XK1asKCZOnFiMGDGiuOaaa4oHHnigqFQqtR6tZhry92AAqH8N9wwGgMYgMACkEBgAUggMACkEBoAUAgNACoEBIEXDBqZSqcQPfvCDqFQqtR6lrrlOg+M6DY7rNDiu03807C9alsvl6OzsjFKpFB0dHbUep265ToPjOg2O6zQ4rtN/NOwdDAD1TWAASDHkn6bc19cXb731VrS3t0dTU9N5v065XB7wTz6a6zQ4rtPguE6DczFfp6Io4syZMzFu3LgP/XmC/2/In8H861//igkTJgzlKQG4wHp6ej7xjxgO+Vtk7e3tQ31KAC6wwXwvH/LAfJq3xQCoD4P5Xu4hPwApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0CK8wrMY489FpMmTYoRI0bEnDlz4sCBAxd6LgAaXNWBeeaZZ2LVqlWxdu3aOHLkSNx4441xxx13xKlTpzLmA6BRFVWaPXt20dXV1f91b29vMW7cuKK7u3tQx5dKpSIiLMuyrAZepVLpE7/fV3UH8/7778fhw4djwYIF/duam5tjwYIF8dJLL33kMZVKJcrl8oAFwMWvqsC8/fbb0dvbG2PHjh2wfezYsXHixImPPKa7uzs6Ozv7l79mCXBpSP8psjVr1kSpVOpfPT092acEoA4Mq2bnK664IlpaWuLkyZMDtp88eTKuvPLKjzymtbU1Wltbz39CABpSVXcww4cPj1tuuSV2797dv62vry92794dc+fOveDDAdC4qrqDiYhYtWpVLF68OGbOnBmzZ8+OjRs3xtmzZ2PJkiUZ8wHQoKoOzDe+8Y04ffp0fP/7348TJ07ETTfdFL/97W8/9OAfgEtbU1EUxVCesFwuR2dn51CeEoALrFQqRUdHx8fu47PIAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASBFVYHp7u6OWbNmRXt7e4wZMybuvvvuOHbsWNZsADSwqgKzZ8+e6Orqiv3798euXbvigw8+iNtvvz3Onj2bNR8ADaqpKIrifA8+ffp0jBkzJvbs2RO33nrroI4pl8vR2dl5vqcEoA6USqXo6Oj42H2GfdoTRESMGjXqnPtUKpWoVCr9X5fL5U9zSgAaxHk/5O/r64uVK1fG/PnzY9q0aefcr7u7Ozo7O/vXhAkTzveUADSQ836L7P7774+dO3fGvn37Yvz48efc76PuYEQGoLGlvUW2dOnSePHFF2Pv3r0fG5eIiNbW1mhtbT2f0wDQwKoKTFEUsWzZsnj22WfjD3/4Q0yePDlrLgAaXFWB6erqiu3bt8dzzz0X7e3tceLEiYiI6OzsjLa2tpQBAWhMVT2DaWpq+sjt27Zti29961uDeg0/pgzQ+C74M5hP8SszAFxifBYZACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKT4VIF5+OGHo6mpKVauXHmBxgHgYnHegTl48GA8/vjjccMNN1zIeQC4SJxXYN59991YtGhRbN26NUaOHHmhZwLgInBegenq6oqFCxfGggULPnHfSqUS5XJ5wALg4jes2gOefvrpOHLkSBw8eHBQ+3d3d8cPf/jDqgcDoLFVdQfT09MTK1asiKeeeipGjBgxqGPWrFkTpVKpf/X09JzXoAA0lqaiKIrB7vyb3/wmvvrVr0ZLS0v/tt7e3mhqaorm5uaoVCoD/t1HKZfL0dnZef4TA1BzpVIpOjo6Pnafqt4iu+222+LVV18dsG3JkiUxZcqU+O53v/uJcQHg0lFVYNrb22PatGkDtl1++eUxevToD20H4NLmN/kBSFHVM5gLwTMYgMY3mGcw7mAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNAiqoD8+abb8a9994bo0ePjra2tpg+fXocOnQoYzYAGtiwanZ+5513Yv78+fGlL30pdu7cGZ/97Gfj9ddfj5EjR2bNB0CDqiowGzZsiAkTJsS2bdv6t02ePPmCDwVA46vqLbLnn38+Zs6cGffcc0+MGTMmZsyYEVu3bs2aDYAGVlVg3njjjdi8eXN87nOfi9/97ndx//33x/Lly+OJJ5445zGVSiXK5fKABcAloKjCZZddVsydO3fAtmXLlhVf+MIXznnM2rVri4iwLMuyLqJVKpU+sRlV3cFcddVVcf311w/Ydt1118U///nPcx6zZs2aKJVK/aunp6eaUwLQoKp6yD9//vw4duzYgG3Hjx+Pq6+++pzHtLa2Rmtr6/lNB0DjquYtsgMHDhTDhg0r1q9fX7z++uvFU089VXzmM58pnnzyyUG/RqlUqvmtnWVZlvXp1mDeIqsqMEVRFC+88EIxbdq0orW1tZgyZUqxZcuWqo4XGMuyrMZfgwlMU1EURQyhcrkcnZ2dQ3lKAC6wUqkUHR0dH7uPzyIDIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBRVBaa3tzcefPDBmDx5crS1tcW1114b69ati6IosuYDoEENq2bnDRs2xObNm+OJJ56IqVOnxqFDh2LJkiXR2dkZy5cvz5oRgAZUVWD+/Oc/x1e+8pVYuHBhRERMmjQpduzYEQcOHEgZDoDGVdVbZPPmzYvdu3fH8ePHIyLilVdeiX379sWdd96ZMhwAjauqO5jVq1dHuVyOKVOmREtLS/T29sb69etj0aJF5zymUqlEpVLp/7pcLp//tAA0jqIKO3bsKMaPH1/s2LGj+Mtf/lL88pe/LEaNGlX84he/OOcxa9euLSLCsizLuohWqVT6xGZUFZjx48cXmzZtGrBt3bp1xec///lzHvPvf/+7KJVK/aunp6fmF8ayLMv6dGswganqLbL33nsvmpsHPrZpaWmJvr6+cx7T2toara2t1ZwGgItAVYG56667Yv369TFx4sSYOnVqvPzyy/HII4/EfffdlzUfAI2qmrfIyuVysWLFimLixInFiBEjimuuuaZ44IEHikqlMujXKJVKNb+1syzLsj7dGsxbZE1FMbS/hl8ul6Ozs3MoTwnABVYqlaKjo+Nj9/FZZACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEghMACkEBgAUggMACkEBoAUAgNACoEBIIXAAJBCYABIITAApBAYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASCEwAKQQGABSCAwAKQQGgBQCA0AKgQEghcAAkEJgAEgx5IEpimKoTwnABTaY7+VDHpgzZ84M9SkBuMAG8728qRjiW4q+vr546623or29PZqams77dcrlckyYMCF6enqio6PjAk54cXGdBsd1GhzXaXAu5utUFEWcOXMmxo0bF83NH3+PMmyIZurX3Nwc48ePv2Cv19HRcdH9B8zgOg2O6zQ4rtPgXKzXqbOzc1D7ecgPQAqBASBFwwamtbU11q5dG62trbUepa65ToPjOg2O6zQ4rtN/DPlDfgAuDQ17BwNAfRMYAFIIDAApBAaAFAIDQAqBASCFwACQQmAASPE/NKMbgtv7MpEAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Yikes!!! :-(\n"
      ],
      "metadata": {
        "id": "eTYiBukM9A1D"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#Package implementations\n",
        "\n",
        "Numpy and Scipy linear solvers both do (P)LU decomposition and then solution. In general they can be accessed through:\n",
        "\n",
        ">numpy.linalg.solve(A,b)\n",
        "\n",
        "and\n",
        "\n",
        ">scipy.linalg.solve(A,b)\n",
        "\n",
        "Examining the options of scipy show that you can choose to:\n",
        "1. specify the matrix type\n",
        "2. overwrite the original matricies (which may boost efficiency)\n",
        "\n",
        "The advent of distributed computing (HPCs) motivated new algorithms that are better suited to large, sparse, systems (You will run in to these names in the future!)\n",
        "\n",
        "- PARADISO (PARallel Direct SOlver)\n",
        "- SuperLU (Supernodal LU)\n",
        "- UMFPACK (Unsymmetric-pattern MultiFrontal method)\n",
        "- MUMPS (Multifrontal Massively Parallel Sparse Direct Solver)\n"
      ],
      "metadata": {
        "id": "ypxXu0OlhhFV"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# prompt: use numpy to solve a linear system then repeat with scipy\n",
        "\n",
        "import numpy as np\n",
        "from scipy.linalg import solve\n",
        "\n",
        "# Define the coefficient matrix A\n",
        "A = np.array([[2, -1, 0],\n",
        "              [-1, 2, -1],\n",
        "              [0, -1, 2]])\n",
        "\n",
        "# Define the right-hand side vector b\n",
        "b = np.array([1, 2, 3])\n",
        "\n",
        "# Solve the linear system using numpy.linalg.solve\n",
        "x_numpy = np.linalg.solve(A, b)\n",
        "print(\"Solution using numpy.linalg.solve:\\n\", x_numpy)\n",
        "\n",
        "# Solve the linear system using scipy.linalg.solve\n",
        "x_scipy = solve(A, b)\n",
        "print(\"\\nSolution using scipy.linalg.solve:\\n\", x_scipy)\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "7eG4Jxj1iqNm",
        "outputId": "2cf66944-2d51-432e-93c7-933063436195"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Solution using numpy.linalg.solve:\n",
            " [2.5 4.  3.5]\n",
            "\n",
            "Solution using scipy.linalg.solve:\n",
            " [2.5 4.  3.5]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# prompt: solve a 200x200 sparse system with sparse solvers then again  with a dense solver\n",
        "\n",
        "import numpy as np\n",
        "from scipy.sparse import diags\n",
        "from scipy.sparse.linalg import spsolve\n",
        "from scipy.linalg import solve\n",
        "\n",
        "# Generate a random sparse matrix\n",
        "n = 200\n",
        "\n",
        "main_diag = np.full(n, 2)\n",
        "upper_diag = np.full(n - 1, -1)\n",
        "lower_diag = np.full(n - 1, -1)\n",
        "A_sparse = diags([lower_diag, main_diag, upper_diag], offsets=[-1, 0, 1], format='csr')\n",
        "\n",
        "\n",
        "# row_ind = np.arange(n)\n",
        "# col_ind = np.arange(n)\n",
        "# data = np.random.rand(n)\n",
        "# A_sparse = csr_matrix((data, (row_ind, col_ind)), shape=(n, n))\n",
        "\n",
        "# Generate a random right-hand side vector\n",
        "b = np.random.rand(n)\n",
        "\n",
        "# Solve the sparse system\n",
        "%timeit x_sparse = spsolve(A_sparse, b)\n",
        "#print(\"Solution using sparse solver:\\n\", x_sparse, \"\\n\")\n",
        "\n",
        "# Convert the sparse matrix to a dense matrix\n",
        "A_dense = A_sparse.toarray()\n",
        "\n",
        "# Solve the dense system\n",
        "\n",
        "%timeit x_dense = solve(A_dense, b)\n",
        "#print(\"Solution using dense solver:\\n\", x_dense)\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "1tRNkmH3iuaQ",
        "outputId": "d7e93eff-304d-4e70-ec26-417a81bf775d"
      },
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "307 µs ± 142 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
            "1.29 ms ± 101 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Dr. Mike's Tips!\n",
        "\n",
        "- Direct solver are your 'black box' for most of your needs.\n",
        "- They are the most robust for ill-conditioned systems.\n",
        "- They scale *terribly* (both in system size and parallelization)\n",
        "- If you use them, start with a small system and work upwards.\n",
        "- Generally speaking you won't see a speedup with parallelization until you get a large # of nodes\n",
        "- Warning: Some implementations (numpy) are sophisticated enough to handle singular matricies as well as non-singular (be careful with the answer!)\n",
        "- Sparse matricies are your saving grace! Do your best to protect them (hence store the LU factors, not the inverse!)"
      ],
      "metadata": {
        "id": "2P5RjZcxmXRA"
      }
    }
  ]
}