02_LU_decomposition.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LU Decomposition\n",
    "\n",
    "Any matrix $A$ has decomposition of the form $A = PLU$ where\n",
    "\n",
    "* $P$ is a permutation matrix\n",
    "* $L$ is unit lower triangular\n",
    "* $U$ is upper triangular\n",
    "\n",
    "This is called the LU decomposition with partial pivoting."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1\n",
    "\n",
    "Find a decomposition $A = LU$ for\n",
    "$$\n",
    "A = \\begin{bmatrix} 2 & 1 & 1 \\\\ 2 & 0 & 2 \\\\ 4 & 3 & 4 \\end{bmatrix}\n",
    "$$\n",
    "using Gaussian elimination without pivoting and compare to the LU decomposition with partial pivoting computed by `scipy.linalg.lu`.\n",
    "\n",
    "First, let us comute by hand:\n",
    "\\begin{align*}\n",
    "\\begin{bmatrix} \\phantom{-}1 & 0 & 0 \\\\ -1 & 1 & 0 \\\\ -2 & 0 & 1 \\end{bmatrix}\n",
    "\\begin{bmatrix} 2 & 1 & 1 \\\\ 2 & 0 & 2 \\\\ 4 & 3 & 4 \\end{bmatrix}\n",
    "&=\n",
    "\\left[ \\begin{array}{rrr} 2 & 1 & 1 \\\\ 0 & -1 & \\phantom{+}1 \\\\ 0 & 1 & 2 \\end{array} \\right] \\\\\n",
    "\\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 1 & 0 \\\\ 0 & 1 & 1 \\end{bmatrix}\n",
    "\\left[ \\begin{array}{rrr} 2 & 1 & 1 \\\\ 0 & -1 & \\phantom{+}1 \\\\ 0 & 1 & 2 \\end{array} \\right]\n",
    "&=\n",
    "\\left[ \\begin{array}{rrr} 2 & 1 & 1 \\\\ 0 & -1 & \\phantom{+}1 \\\\ 0 & 0 & 3 \\end{array} \\right]\n",
    "\\end{align*}\n",
    "Therefore $A = LU$ with\n",
    "$$\n",
    "L = \\left[ \\begin{array}{rrr} 1 & 0 & 0 \\\\ 1 & 1 & 0 \\\\ 2 & -1 & \\phantom{+}1 \\end{array} \\right]\n",
    "\\hspace{10mm}\n",
    "U = \\left[ \\begin{array}{rrr} 2 & 1 & 1 \\\\ 0 & -1 & \\phantom{+}1 \\\\ 0 & 0 & 3 \\end{array} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align*}\n",
    "L_2 L_1 A &= U \\\\\n",
    "A &= L_1^{-1} L_2^{-1} U \\\\\n",
    "A &= LU\n",
    "\\end{align*}\n",
    "where $L = L_1^{-1} L_2^{-1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's verify the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1  0  0]\n",
      " [ 1  1  0]\n",
      " [ 2 -1  1]]\n"
     ]
    }
   ],
   "source": [
    "L = np.array([[1,0,0],[1,1,0],[2,-1,1]])\n",
    "print(L)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2  1  1]\n",
      " [ 0 -1  1]\n",
      " [ 0  0  3]]\n"
     ]
    }
   ],
   "source": [
    "U = np.array([[2,1,1],[0,-1,1],[0,0,3]])\n",
    "print(U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2, 1, 1],\n",
       "       [2, 0, 2],\n",
       "       [4, 3, 4]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L @ U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now compute with `scipy.linalg.lu`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "la.lu?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[2,1,1],[2,0,2],[4,3,4]])\n",
    "p,l,u = la.lu(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 0. 1.]\n",
      " [0. 1. 0.]\n",
      " [1. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "print(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.         0.         0.        ]\n",
      " [0.5        1.         0.        ]\n",
      " [0.5        0.33333333 1.        ]]\n"
     ]
    }
   ],
   "source": [
    "print(l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 4.   3.   4. ]\n",
      " [ 0.  -1.5  0. ]\n",
      " [ 0.   0.  -1. ]]\n"
     ]
    }
   ],
   "source": [
    "print(u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2., 1., 1.],\n",
       "       [2., 0., 2.],\n",
       "       [4., 3., 4.]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p@l@u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With partial pivoting the matrix $L$ always has entries with absolute value less than or equal to 1 in the lower left. That is, $|\\ell_{i,j}| \\leq 1$ for $i > j$ where $\\ell_{i,j}$ is the entry of $L$ in row $i$ and column $j$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2\n",
    "\n",
    "The $N$ by $N$ matrix\n",
    "$$\n",
    "A = \n",
    "\\left[\n",
    "\\begin{array}{cccccccc}\n",
    "2 & -1 & 0 & 0 & \\cdots & 0 & 0 & 0 \\\\\n",
    "-1 & 2 & -1 & 0 & & 0 & 0  & 0 \\\\\n",
    "0 & -1 & 2 & -1 & \\cdots & 0 & 0 & 0 \\\\\n",
    "\\vdots &  & \\vdots &  & \\ddots & & \\vdots & \\\\\n",
    "0 & 0 & 0 & 0 & \\cdots & -1 & 2 & -1 \\\\\n",
    "0 & 0 & 0 & 0 & \\cdots & 0 & -1 & 2 \\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "$$\n",
    "is symmetric and positive definite. Compute the solution of the system $A \\mathbf{x} = \\mathbf{b}$ for\n",
    "$$\n",
    "\\mathbf{b} = \\left[ \\begin{array}{r} 1 \\\\ \\vdots \\\\ 1 \\end{array} \\right]\n",
    "$$\n",
    "by computing the Cholesky decomposition and using the function `scipy.linalg.solve_triangular`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100\n",
    "A1 = 2*np.eye(N)\n",
    "A2 = np.diag(-np.ones(N-1),1)\n",
    "A = A1 + A2 + A2.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "evals,evecs = la.eig(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAVtElEQVR4nO3df6zdd13H8ee7dy2gQIYUpbY9K8bGXxAZOdmsGHPDMG5lsSbwxzA4XAw3I1vclMU4TSRizDQxZJDBmhuYMCEshs3RLEUlV242km6urTDXVbQ63UrrxmbYmCOUsrd/fL+Xnp6ec+/39n7PPef7Pc9HcnLP+X6/957PJ9tefHh/Pt/vJzITSVLzbRh3AyRJ9TDQJaklDHRJagkDXZJawkCXpJa4YFxfvHnz5tyxY8e4vl6SGunQoUPPZObrBp0bW6Dv2LGDgwcPjuvrJamRIuK/h52z5CJJLWGgS1JLGOiS1BIGuiS1hIEuSS1ROdAjYiYi/jki7htwLiLioxFxLCIeiYi31NtMSdJKVrNs8QbgKPDqAeeuAHaWr0uB28ufkjT1DhyAxUWYnS0+L73ftave76kU6BGxDXgH8GfA7w24ZA9wZxbP4n0wIi6MiC2ZebK+pkrSZOkN6l27Bgf3a18LN94Ip07BzAxEwOnTsGkTLCzUG+pVR+i3Ar8PvGrI+a3Akz2fj5fHzgr0iJgD5gA6nc6qGipJk2AptHuDetMmuPXWwcEdAS+9dOYFkFlct7i4zoEeEVcCT2fmoYiYHXbZgGPn7JyRmfPAPEC323VnDUkTa6XRdm9QnzoFd99d/Pz+988O7g0bzgR8/wh96W/XpcoI/a3Ar0XEbuDlwKsj4jOZ+Z6ea44D23s+bwNO1NdMSRq9QaPvYaPt3qDetAne+U544IHBpZVbb4Vnn52AGnpm3gzcDFCO0G/qC3OAfcD1EXEXxWToc9bPJTXBoBAfVibpD/HeoN61C970pmqTn3UH+ZLzfjhXRFwLkJl7gf3AbuAY8CJwTS2tk6QRWCnElyuT9Id4r127zj42quAeZlWBnpmLwGL5fm/P8QSuq7NhklSn1YT4epZJ6jS2x+dK0qitJcTXq0xSJwNdUqvUGeJNY6BLao0DB+Cyy6YrxHsZ6JIab2lU/sQTZ9aCT0uI9zLQJTXSsDXjF5SpNi0h3stAl9Q4w0orAO97H3Q60xPivQx0SY1RpbRy9dXTF+RLDHRJE83SSnUGuqSJZWlldQx0SROl9ymHi4uWVlbDQJc0MXpH5EullE2bzv5saWU4A13S2A2a7Dx1qgjvhYVmPEdlEhjoksaqd1TeP9m5FOIGeTUGuqSxGDQqByc718JAl7RuqixBdLLz/BnoktaFSxBHr8om0S8H7gdeVl7/+cz8YN81s8AXgMfLQ/dk5ofqbaqkJvLuzvVTZYT+XeBtmflCRGwEvhIRX8zMB/uueyAzr6y/iZKaarkJT5cg1q/KJtEJvFB+3Fi+cpSNktRsTniOR6UaekTMAIeAnwQ+lpkPDbhsV0R8DTgB3JSZRwb8nTlgDqDT6Zx3oyVNruVG5ZZWRqtSoGfm94E3R8SFwN9GxBsz89GeSw4DF5Vlmd3AvcDOAX9nHpgH6Ha7jvKlFuq9XR8cla+nVa1yycxvRcQicDnwaM/x53ve74+Ij0fE5sx8praWSppovUsSe2/Xd1S+fqqscnkd8L0yzF8BvB34i75rXg88lZkZEZcAG4BnR9FgSZNn0DNYnPBcf1VG6FuAT5d19A3A32TmfRFxLUBm7gXeBbw/Ik4D3wGuKidTJbXYcs9gufnmcbdu+lRZ5fIIcPGA43t73t8G3FZv0yRNspWewaL1552iklbFJYmTy0CXVJlLEiebgS5pRY7Km8FAl7QsR+XNYaBLWpY3CjWHgS5pIG8Uah4DXdI5vFGomQx0ST/gjULNZqBLArxRqA0MdEmAk59tYKBLU87Jz/Yw0KUp5uRnuxjo0pRZGpHPzp5dZnHys/kMdGmKDBqR95ZZnPxsNgNdmiKDRuQLC2dG7JZZmq3KjkUvB+4HXlZe//nM/GDfNQF8BNgNvAj8VmYerr+5ks7HsInPpRA3yNuhygj9u8Dbyg2gNwJfiYgvZuaDPddcQbEp9E7gUuD28qekMXPic3pU2bEogRfKjxvLV//2cnuAO8trH4yICyNiS2aerLW1klbNic/pUamGXu4negj4SeBjmflQ3yVbgSd7Ph8vjxno0pgsV2ZRO1UK9Mz8PvDmiLgQ+NuIeGNmPtpzSQz6tf4DETEHzAF0Op3zaK6kKiyzTKdVrXLJzG9FxCJwOdAb6MeB7T2ftwEnBvz+PDAP0O12zwl8SWvjw7WmW5VVLq8DvleG+SuAtwN/0XfZPuD6iLiLYjL0Oevn0vry4VqqMkLfAny6rKNvAP4mM++LiGsBMnMvsJ9iyeIximWL14yovZKG8OFaqrLK5RHg4gHH9/a8T+C6epsmqQofrqUl3ikqNZiTn+ploEsN5hpz9TLQpQZyjbkGMdClhrHMomEMdKlhLLNoGANdagjLLFqJgS41gGUWVWGgSxPMW/m1Gga6NKG8lV+rZaBLE8pb+bVaBro0YbyVX+fLQJcmiJOfWgsDXZogrjHXWhjo0gRwjbnqYKBLY2aZRXUx0KUxs8yiumxY6YKI2B4RX46IoxFxJCJuGHDNbEQ8FxFfLV9/PJrmSu1x4ADccsuZMsvMjGUWrU2VEfpp4AOZeTgiXgUciogvZeZjfdc9kJlX1t9EqX0ss2gUqmxBdxI4Wb7/dkQcBbYC/YEuqSLLLBqFVdXQI2IHxf6iDw04vSsivgacAG7KzCMDfn8OmAPodDqrbavUeK5m0ShVDvSIeCVwN3BjZj7fd/owcFFmvhARu4F7gZ39fyMz54F5gG63m+fdaqmBLLNo1CoFekRspAjzz2bmPf3newM+M/dHxMcjYnNmPlNfU6Vms8yiUVsx0CMigE8CRzPzw0OueT3wVGZmRFxCsXrm2VpbKjWUZRatlyoj9LcCvwn8S0R8tTz2h0AHIDP3Au8C3h8Rp4HvAFdlpiUVTT3LLFpPVVa5fAWIFa65DbitrkZJbWGZRevJO0WlEbDMonEw0KWaWWbRuBjoUs0ss2hcDHSpBkslltnZ4mWZReNgoEtr1F9iWVgoXksBb5lF68VAl9aov8SyuFiUWAxyrTcDXTpPrmTRpDHQpfPgShZNIgNdOg+uZNEkMtClVbDMoklmoEsVWWbRpDPQpYoss2jSGejSCiyzqCkMdGkZllnUJAa6tAzLLGqSDeNugDTJlp7LMjNjmUWTb8VAj4jtEfHliDgaEUci4oYB10REfDQijkXEIxHxltE0V1ofBw7ALbcU7xcW4E//tPhpmUWTrErJ5TTwgcw8HBGvAg5FxJcy87Gea64AdpavS4Hby59S4wx62JZlFjXBiiP0zDyZmYfL998GjgJb+y7bA9yZhQeBCyNiS+2tldbBoIdtSU2wqknRiNgBXAw81HdqK/Bkz+fj5bGTfb8/B8wBdDqd1bVUGjGXJ6rpKgd6RLwSuBu4MTOf7z894FfynAOZ88A8QLfbPee8NC4uT1QbVAr0iNhIEeafzcx7BlxyHNje83kbcGLtzZPWh8sT1QYrBnpEBPBJ4GhmfnjIZfuA6yPiLorJ0Ocy8+SQa6WJYZlFbVJlhP5W4DeBf4mIr5bH/hDoAGTmXmA/sBs4BrwIXFN/U6V6WWZR26wY6Jn5FQbXyHuvSeC6uholrQfLLGobb/3X1LHMorYy0DVVLLOozQx0TRXLLGozA11TwTKLpoGBrtazzKJpYaCr9SyzaFoY6GotyyyaNga6Wskyi6aRga5WssyiaWSgq1Uss2iaGehqDcssmnYGulrDMoumnYGuxrPMIhUMdDWaZRbpDANdjWaZRTrDQFcjWWaRzlVlC7o7gCuBpzPzjQPOzwJfAB4vD92TmR+qs5FSL8ss0mBVRuifAm4D7lzmmgcy88paWiStwDKLNFiVLejuj4gdo2+KNNxSiWV2tnhZZpHOVVcNfVdEfA04AdyUmUcGXRQRc8AcQKfTqemr1Xb9JZaFheK1FPCWWaRCHYF+GLgoM1+IiN3AvcDOQRdm5jwwD9DtdrOG79YU6C+xLC4WJRaDXDrbhrX+gcx8PjNfKN/vBzZGxOY1t0xT78ABuOWWMytZZmYssUjLWfMIPSJeDzyVmRkRl1D8j8Sza26ZpporWaTVq7Js8XPALLA5Io4DHwQ2AmTmXuBdwPsj4jTwHeCqzLScojVxJYu0elVWubx7hfO3USxrlNbMG4ak8+edopoYllmktTHQNTEss0hrY6Br7CyzSPUw0DVWllmk+hjoGivLLFJ9DHSNhWUWqX4GutadZRZpNAx0rZulUfkTT1hmkUbBQNe66B2Vz8zABeW/eZZZpPoY6FoXvZOfAO97H3Q6llmkOhnoGqlhk59XX22QS3Uz0DUyTn5K68tA18i4xlxaXwa6aucac2k8DHTVyjKLND4rbkEXEXdExNMR8eiQ8xERH42IYxHxSES8pf5mqimGlVkMc2n0quwp+ing8mXOX0GxKfROYA64fe3NUtO4/6c0flV2LLo/InYsc8ke4M5y27kHI+LCiNiSmSdraqMmnGUWaTLUUUPfCjzZ8/l4ecxAbzlv5ZcmSx2BHgOODdwkOiLmKMoydDqdGr5a4+Kt/NLkqSPQjwPbez5vA04MujAz54F5gG63OzD01Qzeyi9NnjoCfR9wfUTcBVwKPGf9vL28lV+aXCsGekR8DpgFNkfEceCDwEaAzNwL7Ad2A8eAF4FrRtVYjZeTn9Jkq7LK5d0rnE/gutpapInj5KfUDN4pqmU5+Sk1h4GuZTn5KTWHga6BnPyUmsdA1zmc/JSayUDXDzj5KTWbgS7AyU+pDQx0AU5+Sm1goE+xpRLL7GzxcvJTajYDfUr1T3wuLBSvpYA3zKXmMdCnzLCJz8VFdxaSms5AnyJOfErtZqBPgUGjcnDiU2obA73llhuVO/EptYuB3nIuR5Smh4HeUj6LRZo+BnoL+SwWaTpVCvSIuBz4CDADfCIz/7zv/CzwBeDx8tA9mfmhGtupCnwWizTdqmxBNwN8DPgVig2hH46IfZn5WN+lD2TmlSNooypwSaKkKiP0S4BjmfmfAOVm0HuA/kDXGLgkUdKSKoG+FXiy5/Nx4NIB1+2KiK8BJ4CbMvNI/wURMQfMAXQ6ndW3VmdxSaKkXlUCPQYcy77Ph4GLMvOFiNgN3AvsPOeXMueBeYBut9v/N1SRo3JJg1QJ9OPA9p7P2yhG4T+Qmc/3vN8fER+PiM2Z+Uw9zdQSR+WShqkS6A8DOyPiDcA3gKuA3+i9ICJeDzyVmRkRlwAbgGfrbuw0c1QuaSUrBnpmno6I64G/p1i2eEdmHomIa8vze4F3Ae+PiNPAd4CrMtOSyhr13hx0442OyiUtL8aVu91uNw8ePDiW726C3tJKBLz0UvGamXFULk2ziDiUmd1B57xTdMIMKq1s2FAEeYSjcknDGegTZLkJT2/fl7QSA30COOEpqQ4G+pg44Smpbgb6GAyb8ARH5ZLOn4G+jpzwlDRKBvoILQX40tMOnfCUNEoG+oj0bzLx3vc64SlptAz0mg3bZALcCk7SaBnoNai6YuXqq8+UYAxzSXUz0M/ToBCvsmLFIJc0Kgb6KqwU4q5YkTROBvoKVhvirliRNC4G+gCGuKQmMtBLhrikppu6QO+/2ccQl9QWlQI9Ii4HPkKxY9EnMvPP+85HeX438CLwW5l5uOa2rspKwb0U1KdPG+KS2mHFQI+IGeBjwK9QbBj9cETsy8zHei67AthZvi4Fbi9/1m5QUC+F7bD14IOCe2lpYaYhLqkdqozQLwGOZeZ/AkTEXcAeoDfQ9wB3lvuIPhgRF0bElsw8WWdj+zeAWArqpRAeVDYZFtyDft8Ql9RkVQJ9K/Bkz+fjnDv6HnTNVuCsQI+IOWAOoNPprLatLC6euZ2+N6hPnYK77x78BMPlgnvpbxriktqgSqDHgGP9O0tXuYbMnAfmodgkusJ3n2V29szzUPqD+p3vhAceOPOslKrBbZBLaosqgX4c2N7zeRtw4jyuWbNdu2BhYXgN/U1vMrglTa8oyt7LXBBxAfBvwGXAN4CHgd/IzCM917wDuJ5ilculwEcz85Ll/m63282DBw+urfWSNGUi4lBmdgedW3GEnpmnI+J64O8pli3ekZlHIuLa8vxeYD9FmB+jWLZ4TV2NlyRVU2kdembupwjt3mN7e94ncF29TZMkrcaGcTdAklQPA12SWsJAl6SWMNAlqSVWXLY4si+O+Cbw3+f565uBZ2psTlNMY7+nsc8wnf2exj7D6vt9UWa+btCJsQX6WkTEwWHrMNtsGvs9jX2G6ez3NPYZ6u23JRdJagkDXZJaoqmBPj/uBozJNPZ7GvsM09nvaewz1NjvRtbQJUnnauoIXZLUx0CXpJZoXKBHxOUR8fWIOBYRfzDu9oxCRGyPiC9HxNGIOBIRN5THfyQivhQR/17+fM2421q3iJiJiH+OiPvKz9PQ5wsj4vMR8a/lP/NdU9Lv3y3//X40Ij4XES9vW78j4o6IeDoiHu05NrSPEXFzmW1fj4hfXe33NSrQezasvgL4WeDdEfGz423VSJwGPpCZPwP8AnBd2c8/ABYycyewUH5umxuAoz2fp6HPHwH+LjN/Gvh5iv63ut8RsRX4HaCbmW+keDT3VbSv358CLu87NrCP5X/jVwE/V/7Ox8vMq6xRgU7PhtWZeQpY2rC6VTLzZGYeLt9/m+I/8K0Uff10edmngV8fTwtHIyK2Ae8APtFzuO19fjXwy8AnATLzVGZ+i5b3u3QB8IpyE50fotjlrFX9zsz7gf/tOzysj3uAuzLzu5n5OMX+EstuFNSvaYE+bDPq1oqIHcDFwEPAj2XmSShCH/jR8bVsJG4Ffh94qedY2/v8E8A3gb8qS02fiIgfpuX9zsxvAH8JPEGxmfxzmfkPtLzfpWF9XHO+NS3QK21G3RYR8UrgbuDGzHx+3O0ZpYi4Eng6Mw+Nuy3r7ALgLcDtmXkx8H80v8yworJuvAd4A/DjwA9HxHvG26qxW3O+NS3Q12Uz6kkQERspwvyzmXlPefipiNhSnt8CPD2u9o3AW4Ffi4j/oiilvS0iPkO7+wzFv9PHM/Oh8vPnKQK+7f1+O/B4Zn4zM78H3AP8Iu3vNwzv45rzrWmB/jCwMyLeEBGbKCYQ9o25TbWLiKCoqR7NzA/3nNoHvLd8/17gC+vdtlHJzJszc1tm7qD45/qPmfkeWtxngMz8H+DJiPip8tBlwGO0vN8UpZZfiIgfKv99v4xirqjt/YbhfdwHXBURL4uINwA7gX9a1V/OzEa9KDaj/jfgP4A/Gnd7RtTHX6L4v1qPAF8tX7uB11LMiv97+fNHxt3WEfV/FrivfN/6PgNvBg6W/7zvBV4zJf3+E+BfgUeBvwZe1rZ+A5+jmCP4HsUI/LeX6yPwR2W2fR24YrXf563/ktQSTSu5SJKGMNAlqSUMdElqCQNdklrCQJekljDQJaklDHRJaon/B5On9hBF00kUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.sort(evals.real),'b.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "U = la.cholesky(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAD7CAYAAADgvbh3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXCcd5ng8e/Trfu0bFmyLclX4pyEHJgQSCAXGULIjKmFhGMDGY4JWSCEMwkztcwOW1PF1E5NDdRQsC6OCTsskArskGIzQLAT8BE7tuM4vuLEpyxbtiQfkmxd3epn/3i614rj45W71f2+3c+nqkv9dre6f3GkR8/vfERVcc65YhIrdAOccy7XPLA554qOBzbnXNHxwOacKzoe2JxzRccDm3Ou6GQV2ETkThHZISI7ReSxXDXKOeeyIRe6jk1E4sCrwB1AF7AO+Iiqbstd85xzbvLKsvje64GdqrobQER+DiwBzhrYpLlGaZ8O+xvhRAVUJ6BiPIsmOOdeZ7AbHT4u2byF3Hmx0jcU7MUbun+nqndm83lTIZvA1gbsn3DdBbzt9BeJyAPAAwB0TIOnvgL//WZ4cTa0D8CMIShL2c05l50nP579e/QNwfrPBHut/Lfm7D8w97IZYzvTX4U39GtVdamqLlbVxdAMv7kE3rMTvvFHmD4Ma9stg3POhYcGvIVUNoGtC+iYcN0OHDznd4yUwZ4mmDMI13VbV/RYFQyVw3gMUlll0M65nBD7XQxyC6lsuqLrgEUisgA4AHwY+Oh5v2ssDi+0wfaZcNVhuKwPXm6F9XOsW7rwGMRC/KfAuWKngIY3aAVxwYFNVZMi8nngd0Ac+JGqbj3vNyZi8OoMG1P7s11wzSHoq4Enr4BkDOYf98DmXKFF/Fcwm4wNVX0aePqCvjklsGu6dUOnD8MnN8LRatjRbF3UjgEo9xlT5woi4hlb4XYepAR2zIDn5sOMYXjoBbj+AGxvtscTvinCuYKJ+ORBVhlbzvTVwCszLNjdvgdG49BTC3GFmSehKlnoFjpXWjxjy4HdTfC7i+0f89GVcM82G4fbMBtOlhe6dc6VFgVSAW8hFY6MbVxgPG7jbcer7PqSIzBYAcPl0BuD+lHP3JzLF8/YcuhgPfzHxfb1S2vgi2utS7pyrnVXnXN5IBbYgtxCKhwZW0YiBolKqElCZRJqBBpHoXHEloIMVtrjvr/UuakV4omBIMIV2DKOVcHyBVCdhPtfssf+91U2g3pFL1x0tKDNc66olfIC3Sk1lp4VbRy1JSC1Y5alDZfZjOlY3GZM4yEevXQuykK8XSqIcAa2jKFyWNtmkwa37bHbqrmwYh60DcClR0AinjM7F0YR/7UKd2BLxKCzESrH4fbd1gXd3GqP1STsr0oMD27O5ZJ3RfMkGYNXmuFwnQW3Lz1vwW3DbGgYhYuP+nluzuVSxHOFcC33OJtxsUW8G2fB3H74y5dg0VHYNhP2TbPA55zLHV/ukUcpsTVucbWu6L1b4WQF7J1mkwvtA74UxLlc8Iwtz3Y3wZ/mWRf0K8/b2NuWFjvTbci3XzmXNS3tgyYLI/OPOVAJXfV2/bYDtgzkeJVlcDOGfPuVc9kIcTcziOhlbBn7G+CZiyxLe2SVjbvtaYLVHRbgnHMXzo8tKpBkzG4jZbYsJKa2tq0qaQt4j1bbOJxnbs5NXsQztugGtozeWvjDQqgfg8+us6D23bfCujZ4y0E7atw5NzkhzsaCiH5gS8Ss6xlXaBqxpSE1SVvXlojbsUdlKT9m3LmglFBPDAQR3TG20w1WwIq5sGEOfHAbfPNZ22O6fAHsmVbo1jkXLRFfx1Y8gW0sDofq7Ny2uf1wbbcdcdRXAycqbDwu4n+FnMsbnzwImdE4vDTLsrW3HoRrD9kat1VzofWEnczr5f2cO4dwZ2NBFF9gS8ZgX6ONq92yFy7vhe56q6GQwvaVemBz7uxCno0FUXyBLSMltq7tRAXMHoQH11u3dHML1CZgwXGfUHDubCKesRXPGNvpUgK7mmBNO8w6Af9lnVWd39xqRZnH4oVuoXPh5VuqIuBIjQWzshS8d6eNw3XX2fWsE76I17nTRbwrWrwZ20R7plm9hPJx2371/lfsyKN1bV631LnTZQ6a9OUeIZcSW8g7XG6FYlICl/fBxUds0/yhOtua5ZwzOVruISI/EpEeEdlyludFRL4jIjtF5GURuS4XzS+NwJZxoN42zh+pga+uhofXwoEGOwbJ65Y6d0ru1rH9K3DnOZ5/L7AofXsA+N4Ftvh1SiuwJWPW9RyNW7e0KglNw7YVKxmzo5B8UsG5nHVFVfVPwLnqZS4BfqJmDTBNRGZn2/zS7H8drbYxt6okfHKjPfaTq2HZQrjqsK11c65UaV5nPNuA/ROuu9KPdWfzpqUZ2Mbi1vVsHD01K1oxfqpm6UiZzZh6gRhXqoJPDDSLyPoJ10tVdekkPulMH5T1nGxpBraMoXJb51adhPfsgjt3Wib3x/kw77jXLXWlK/iPfZ+qLs7ik7qAjgnX7cDBLN4PKLUxttMlYlYcpqsBFhyDd+2DujE7nfdYtaXjIZ7Sdm7K5G+5x1PAx9OzozcA/aqaVTcUSj1jy0jGYPtMC3KX9dmM6d5p8EKbTS543VJXanLUURGRnwG3YF3WLuBvgXIAVf0+8DRwF7ATGAI+kYvPPW9gE5EO4CfALGwb+VJV/baITAd+AcwH9gL3quqxXDQq78bFFvHGFG7qhKsPwf+6Gp68AjoG7BReD2yuVOTwoElV/ch5nlfgczn5sAmCdEWTwFdU9XLgBuBzInIF8BiwTFUXAcvS19GWElvXtmmWlff76GY7123XdNg53ZeCuNJR7DsPVLVbVV9M3x8EtmPTsUuAx9Mvexx4/1Q1Mq92NcHKuba27cvPw8377ESQl1t9d4IrHaV00KSIzAeuBdYCrZlBPlXtFpGWs3zPA9iKYqiblUVT8ygltli3s9H+593QZfUTjtTY481et9QVs3BnY0EEnhUVkTrgl8AXVXUg6Pep6lJVXayqi6mOUO2BzkZbsDsWh6+vhI9vgtem20m8XrfUFbOg2VrUMzYRKceC2k9V9Vfphw+LyOx0tjYb6JmqRhbEuMB4erHuaNzWs3UMWJm/kTI4Um0HVnrm5opRsWdsIiLAD4HtqvpPE556Crg/ff9+4Ne5b14I9NZa5tbZCF9YC19bBf1VsHyhnQriXDEqgYMmbwQ+BmwWkZfSj/018C3gCRH5FNAJ3DM1TSywRAz6K225R8MoVCeswnxl0sbdTlbYdiw/ZtwVkxB3M4M4b2BT1ZWceT8XwO25bU6IDVbYbGllEu7damve/v0yWLbAFvVecqTQLXQuNzIHTUZYaW+pmoyxOByutZnRjn47BaQyaVuvTnrdUldkSmHywE2QqVtak7BlIIsP2vWKeXZSyKV9Xt7PRV/EMzYPbJOVjNlEQnkKbt5rXdDORlsKgsIiL8jsikDEf4Q9sF2ocYHdTTZD2j4An10HPbW2Hat2DC465hMKLprye9DklPAxtguVSge2dXOgbRAe2ABvPmxbr16bYbOpzkVVxPeKesaWC73puqXlKbj7VRuH62qw69mDvojXRU/Eu6KeVuTC3mnwx3kWwB5dCX/+KmxtgfVzvG6piybP2Nz/X4U9VG61FFJiy0GGymGwEsbK7MBKz9xcVEQ8Y/PAlksH6i2wNQ/ZKbw9tfAPN8LBBrix0yYZnAu7HB40WSge2HIpGbNbfRziahla87BlbMmYnQpSlfTMzYVcuLuZQXhgmwpHaqzaVXUCPvWi/ZD8+BrbTH/NIbjI65a6kPOuqHuDRMyONZomthuhPAWV45a1jXrdUhcBHtjcWZ2ssLqlleNw12vwvldh+QJ4doGV+/O6pS6MfBO8O6dM3dLuOgtkN3Xa4ZQH6uFYldctdeHlm+DdeSVjsG2mFWK+sgcePW5Vr9a0w/Rh228a926pC5FUtHOeaLc+KpIxW8S7rQUWHIePbLajjzLBbtyzNhcynrG5wBTbapWIWXm/j22yBbw7mm0JyLzjdhqvc4XkY2xuUlJidUvXtMOMIfjiGnjHfts4v6XF65a68PCMzU1aSuy4o33TbFb0nfvshN6+Gnt85klfxOsKKPqTWp6xFUpnoy39SAk8thI+uhl2zIDVHV631BWeZ2zugkysWzpcbuVy5vfD8RF7rK8G6sY8c3P5VwR7RT1jK7TDdVbp6kC9jbl9bbVtyVq+wOuWusLxY4tcVhIxSFTa7oS6MTtOvHYMqpO2TOREhT3nx4y7fApxNzMID2xh0V8JK+ZaEPvPm21S4ckr4A8L4Yper1vq8ivE2VgQHtjCYixu57fVjcF13dA4YmvaBiptDC4Rt+pXvkPB5YNnbC6nRstg4yybNHhnp61zWz8H/jQP5gx63VI39Ypg8sADW9gkJtQtvXUPXHwUdk23vaVlKa9b6vIg3BMDQXhgC6uUWDA7Wm1brR5aa7OkG2dBvdctdVMs4n87fblHWGUKMq+fAx0D8OkX4cpe2NxqAS/p/+vcFIr4cg//7YiCnlp4pRkqk7DkFbj2EOxrtMDn+0tdrgXddRAwqxORO0Vkh4jsFJHHzvD8LSLSLyIvpW/fyPY/wX8romDPNDv26NI+eGS1ZXH/9Vb7i9m023cnuNzLUTYmInHgu8AdQBewTkSeUtVtp710harenZMPxQNbNGRmqE5WQE8NpLCiMMNltv5tuMwOrPQA53Ild7Oi1wM7VXU3gIj8HFgCnB7Ycsq7olFyoB6euQgGquCRVfCZDXZCyJ/m2TYs53Jhcl3RZhFZP+H2wGnv1gbsn3DdlX7sdG8XkU0i8h8icmW2/wmesUVJpm7paNyuK8btiKNxsWUix6qt5J9nbi5bwbuifaq6+BzPn+mNTh+dexGYp6onROQu4N+BRUEbcCaBMzYRiYvIRhH5Tfp6uog8IyKvpb82ZdMQNwmZuqU7p1vW9jcrbHfCsgVwoKHQrXPFIHeTB11Ax4TrduDg6z5KdUBVT6TvPw2Ui0hzNs2fTFf0YWD7hOvHgGWqughYlr52+ZCI2fq2ExWWsbUN2IxpSmxr1kiZLwdxWQi41CNYVrcOWCQiC0SkAvgw8NTrPk1klohI+v71WFzKanN0oJ9+EWkH3gf8YMLDS4DH0/cfB96fTUPcBcjULV0/B+5+Fb75LEwbsYrzO6cXunUuqjJbqoLczvdWqkng88DvsMToCVXdKiIPisiD6Zd9ENgiIpuA7wAfVtWslggHHWP7Z+ARoH7CY62q2p1ufLeItJzpG9ODiTagWDfrwlvq3ihTt7QqCbftsb2kK+faDoWZJ2E8ZqeE+BYsN1k5/JFJdy+fPu2x70+4/y/Av+TuEwMENhG5G+hR1Q0icstkP0BVlwJLAaTlcv8NmwqZuqWdjfDmw7DoKLyaPma8ecjrlrrJi/hvapCu6I3AX4jIXuDnwG0i8m/AYRGZDZD+2jNlrXTnlqlbun0mLDwG92yF2YMW7LoaIn9SgyuAYt9SpapfV9V2VZ2PDfwtV9X7sAHA+9Mvux/49ZS10gWjWNb20izL1D6x0Rbybm+2DG4sXugWuqiIeDGXbKbOvgXcISKvYdslvpWbJrkLlkpvnF/dAS0n4Qtr4YYu2DQLtrbYWW/OnU+mYHKEM7ZJ/aSr6nPAc+n7R4Dbc98klxPHqmFPk00e3Lrn1Am9R6ph1glfxOvOLeLDF77YqVh1NsKz8+3+YyvhQ1vthJC17XbcuHNnldN1bAXhfZNiNbFu6WCldS8uOmr3h8ote2sY9czNnVmIx8+C8Iyt2B2ug+Xz7euX1sDDa6C3Fp71uqXuLEptjM1FUKZuadW4bZAHqB+1aljjMcvgKpO2od65jIhnbB7YSkV/JayYZ13P+162Mxd+caXVLb2yx+uWutcLcTYWhAe2UjEWh94ay9bectAytopx20g/UuZ1S93rRXxW1ANbqRkpt03zVUm4ZS/cvM820j83H9oHvG6pC/3i2yA8sJWaRMy2WVWMw617bQvW9mZb2FuVhEuEyP9Uu+x5V9RF0njMtln11lhwe3iNnRSyfg40jlih5jLvlpasiP9t8+UepWpcrPrVi7Nhbj98ciNcegS2tFjleT+osoT5Al1XDA7V2Ukg1Qn4wDY7ZnxPk3VX2wZ8EW+pyRw0GWH+Z9nZkUcr5tlM6VdXwx27YHOLZXPD5YVunSuEiJ/u4RmbO/XXeajcxtlSAm/ptmUg/ektWDOGPHMrJSHuZgbhGZs7pasBli+wQPboSvirDTZbumKuFY9xpcMzNlc0MnVLx+KWtZWl7IgjUXvM65aWjohnbB7Y3Bv11lrmVpOAz66zoPY/3wLr2mzXwsJjhW6hm0ohz8aC8MDm3igRg+NVdn/6sH2tGj+VuQ2XWzZX7hvni5bPirqidbLCjhnfMAf+0zarW1o/ZhXndzUVunVuykR/HZsHNnd2iZitcTtcC/P6YfFBqErY2W6DlTYeF/G/7O4sfPLAFb1kDDa3Qu0YXHcIruiz/aWr5lrRGK9bWlwyB01GmGds7vySMdjXCDuabeLgA9ug9aTVUOiq96ytGHnG5kpGpm7pcBnMPAmfetEmGbbMhNoELDjuEwrFwjM2VzJSYpMGa9phziA89AK89SC8PMv2mo56QeaikQp4CynP2NyFOVptQS6mtrd0LG4TDZlFvb6IN7p8jM2VrH3TrNJVTK1u6Qe3Wda2tt1mTF20RXy5h2ds7sKMi91GyqC/yrqplxyxtW8ny2G8DqaNeOYWVSGeGAjCMzaXnUN1tmC3txa+/Dw8tBa66+CP860os4ug6C/Q9YzNZScRg0SF7SutTkJyDBpHrcp8MgYDlZa1ed3S6CiCgyY9sLncOF4Ff0zXLf3Ll+yxn11ldUuv6oFFXrc0UiLeFfXA5nJjLA59NZapLT5oGVzFuG2YHymz5+NetzQyQtzNDMIDm8ut4XTd0spxuG2P3VZ12Azq3H64rM9OCXHhFvH/RT554HIrU7d0f4Ntv7ptj2Vx+xqhr9rGbiKeDZSEHE4eiMidIrJDRHaKyGNneF5E5Dvp518Wkeuybb5nbG5qZOqWHq61GqVfft62Y70wB5q8bmmo5XAfqIjEge8CdwBdwDoReUpVt0142XuBRenb24Dvpb9eMM/Y3NTI1C19udW6oB/fZMFsS6vVUfC6peGWkmC387se2Kmqu1V1DPg5sOS01ywBfqJmDTBNRGZn03zP2NzUSomtddvSYhvlP7TFisXsTtctbfe6peEzqeGCZhFZP+F6qaounXDdBuyfcN3FG7OxM72mDegO2ojTBQpsIjIN+AHwJixJ/SSwA/gFMB/YC9yrqn4Yvnuj3U22BeuKHvjK8zaZ8I1brSva7GX9Qil4V7RPVRef4/kzRcjT3z3IayYlaH/g28BvVfUy4GpgO/AYsExVFwHL0tfOndm42B7SA/WWDbztALypxypfdTXYkhAXDplN8LmZPOgCOiZctwMHL+A1k3LewCYiDcC7gB8CqOqYqh7H+sWPp1/2OPD+bBriSkCmbulIGTyyyhby7myClXMtwLnwyN1Bk+uARSKyQEQqgA8DT532mqeAj6dnR28A+lX1gruhEKwruhDoBX4sIlcDG4CHgdbMh6tqt4i0nOmbReQB4AEA6mZl01YXdZm6paNx+xpTaBu0ClijcThSbeNw3jUtvBwtyVHVpIh8HvgdEAd+pKpbReTB9PPfB54G7gJ2AkPAJ7L93CCBrQy4DnhIVdeKyLeZRLczPZC4FEBaLo/4sj+XE7218Ox8qBuDz62zoPa9t8ILbXZw5QIfqi24HO4VVdWnseA18bHvT7ivwOdy9oEEG2PrArpUdW36+kks0B3OTMmmv/bksmGuiGXqlg6VQ9MwzByyCvNlKXtuqBwSfhpvwQTthoY4TTlvYFPVQ8B+Ebk0/dDtwDasX3x/+rH7gV9PSQtd8RqstLqlG2fZQZXffNZOCPnDQq9bWmglcmzRQ8BP04N/u7E+cAx4QkQ+BXQC90xNE13RytQtrUnA1YegeRh+f5Ftpm8bPDUOFwtxalCsIv5PHiiwqepLwJnWqtye2+a4kjQWh02zLMBdfwCuOQSbW2DFPKufcGmfB7e8Cnc2FoQvHnKFl4zZPtKyFNy8zwLZ4Trba6rYWW4e2PLHD5p0LodS6f2lgxXQegI+s96WgLzcatncRce8bmm+RPzviO9EduGREtt+9UKbjbF9dh1cc9gC245m67K6/CiRyQPn8utItXVFy1Jw504Lagfr7Xr2oC/inWqesTk3BfZOs0pXFePw6CpYsgO2zrRs7kRFoVtX/Dxjc24KZM77Gi6DY1WWQVzWZ4t3T1TYhIPXLZ0aIV98G4QHNhduB+vhSA3MGIKvrbL7/+NGO3r8xv3Q0V/oFhYnnxV1bgplNs7XlUF5ymZHm4ZPZW3Hq05VxHI5Eu5uZhAe2Fw0HK22MbfqBHxqo3WVfnI1LFtouxYuPlrgBhYZ74o6lweJmM2UNsRsN0Jlurr8aBxGy+yMt7KUF4jJhcxBkxHmgc1Fy3A5rG2zSYP37LKlIM/Nt7ql84973dJcifg/oS/3cNGSiNmEwoF03dJ37bNz3fY32Oyp1y3NjdxVqSoIz9hcNCVj8EqzBbnL+uCrq23t29o2q1u66Ih3Sy9UEXRFPWNz0ZSpW7q5BRYch/tetq9bWizAed3S7ET8oEnP2Fy0pcSyNoCGUfjoZlsKsmu6TTDM7felIBfCMzbnCmxXk9UqbRqGLz8Pt+yFTa12xttQeaFbF02esTkXAimBgUrY32jZxjv228b5Y9V2BPnMk779KrDoT8B4xuaKx/5Gq1uaiMFjK+Hjm+yEkFUdtkPBBZM5aNJnRZ0LgXGB8bgt1h0tA8H2ktaN2WN9NXbfM7fzC3E3MwjP2Fzx6a21zG1/A3xhrW2eP1Ztj3XXF7p10eDHFjkXMokY9FfaOraGUdtfWjtms6PJGJyssPt+zPjZRTxj88DmitdAJayaa8s+7t1qXdNfX2p1Sy/ts4W97o18ga5zIZaIweFaO8Otox/efNgytWPpKvTJWKgHwAvKl3s4F3Kj6bql1Ul4exe89SC81Ap/mgezvW7pGUU84Htgc8UvU7e0PAU374VFR21i4dUZ1j31uqWnCffEQBAe2FzpGE+X9+uvgjnp8n69NZbN1Y553dKMkHczg/AxNlc6MnVL182B9gF4cD1c1WPbr16dYWNyzvhyD+ciqK/Gjj2qGIe7X7UFvV0NtkRkjtct9YzNuSjaO80mD6qS8MgqC25bWmD9HDjpG+d9S5VzUZT5xRwqt1oKCrypx44eH6yEsTI7LaQUM7ciWMfmgc2VtoP11i1tHoKvrbbJhH+4CQ7Uw02dNhZXiiLeFfXA5kpbpm7pWBziKagct+LMI3FIxO1UkKpk6WVunrE5VwSO1Fi1q5oEfPpF66b+6zW2/eqaEqxbmoeMTUSmA78A5gN7gXtV9dgZXrcXGATGgaSqLj7fe/vkgXNgSz2OVtux4q0noWPg1Kb5sfRRSKVURyE/W6oeA5ap6iJgWfr6bG5V1WuCBDUIGNhE5EsislVEtojIz0SkSkSmi8gzIvJa+mtTkPdyLtROVlilq/Vz4K7X4O+eg+nDduTRazMi30ULRAVSsWC37CwBHk/ffxx4f7ZvmHHelolIG/AFYLGqvgmIAx9mctHWuWjI1C3trrO6pTd1Wvf0QAMcryyduqXBM7ZmEVk/4fbAJD6lVVW7AdJfW87Rmt+LyIag7x90jK0MqBaRBFADHAS+DtySfv5x4Dng0YDv51y4JWOwbabtKb2yFx5dCTunw+p2aB6GS47YZEOxCh68+87VPRSRPwCzzvDU30yiNTeq6kERaQGeEZFXVPVP5/qG8wY2VT0gIv8IdALDwO9V9fci8rpom/7QN0hHWIuydWf673MuhJIxW8QbV7ix09a4PX4NPHkFzOuHi45a36VY5WjyQFXffbbnROSwiMxOx4/ZQM9Z3uNg+muPiPwf4HrgnIEtSFe0CesLLwDmALUict/5vm9Co5aq6mJVXUz1tKDf5lw4jItttXppli3Y/djLcG037Gi2MbexIoxumQW6U79X9Cng/vT9+4Ffn/4CEakVkfrMfeDPgC3ne+Mgo3/vBvaoaq+qJoBfAe8ADqejLOeKts5F3u4mWNNui3i/uMZK+21qtSr0I0W6Yio/s6LfAu4QkdeAO9LXiMgcEXk6/ZpWYKWIbAJeAP6vqv72fG8c5P9KJ3CDiNRgXdHbgfXASSzKfouzRFvnikZK7LijfY12httNnbaAt6/GFvG2FFnd0jzsA1XVI1g8Of3xg8Bd6fu7gasn+95BxtjWisiTwItAEtgILAXqgCdE5FNY8Ltnsh/uXKR0NtpWq44B+Hp6MuGbN1vhmFv2wqwThW5hjkR/5jdQHq2qfwv87WkPj3KGaOtc0crULR2Nw3AZiMK845bJjZRZ2b/60ehnbn7QpHMl6HCdLdg9WG9jbl9Jb55fvgAO1RW6dbnhB006V2ISMUhU2parujE7nLJ+zOqXJmO2Lasy4nVLI56xeWBz7kL1V8HKuRbEPrrZJhV+ebltnL+i1xbxRlWID5EMwruizl2oRAx6aq1OaceABbPylAW8k+U2azoewV+x/K1jmzKesTmXrdEyeHG2dUXfuc/Wua2fA8/Ng/bBaNYtjVhzT+eBzblsJWK2p7Q8Zcs+Lj4Ku5pg93Qbf4ti3dIQZ2NBeGBzLldSYgHtaDXM7YfPvwCHa2HjbFsGEqW6pRGLw6eL4ACAcyGVKci8YbYFtr/aYONuL7favtLIHFQZcHwtxFmdZ2zOTYWeWtg+02ZMl7xi43D7Gm2JSNjrlio+K+qcO4NM3dKahNUtfe9rsLnVJhWGIlC3ND+b4KeMZ2zOTYVMxjNUbtlbSqwozEgZDFTa1+khrlsa4m5mEB7YnJtKB+ptAqH1pJ3Ce6gO/v6dti3rnZ3QFtK6pSHOxoLwwObcVMrULR1NH0hZnoKWIcuIEjE4Vm3r38KUuXkleOdcIJm6pVVJeHC9BbsfXmtjbtcesqPGw8QzNufceWXqlk4bsZN4Y2ozpik5Vbe0LF7Lkd4AAAXWSURBVGW3gpPIz4p6YHMunzJ1SyvG4c932O0PC2HZQsvaLusrdAtNxLuivtzDuXzK1C09VAfzj8Pbu6A6aXVM+6ts03yhs6WgSz1C3F31jM25QkjGYGuLrXe7+pDtJ311BqzusK5qoeuWesbmnJu0ZMx2IuxohgXH4d6ttiMhU6Q5KllbSHnG5lwhKRbIEjHL1D6x0bqk22baMpD5x208Lt8KHViz5Bmbc4WUSm+cX90BM4fgC2vhbV22cX5ri+0xzTc/aNI5lzPHqyzIxRRu3WPLQHpq4Ui1lfbL5yLeEHczg/CMzbmw6Gy0RbxgG+fv3QqvNMPadttfmjd+bJFzLlcydUtHymy9G9jatoHKU5vpG/JUt9QzNudcTh2ug+Xzba3bF9fYracWns1j3VLP2JxzOZWpW1o1bue5gdUtrRuzBbyDlVCZnLrZ0iI4aNIDm3Nh1V8JK+ZZEPvYJnvsiSttC9aVPVNbtzTiXVEPbM6F1VgcemusEMxbuqF2zLK0wQoYTtctjenU7FAIcTczCA9szoXdSLkViKlKWnm/d+2zmdJn51uh5qmoW+oZm3NuSiVi0NVg2dote2HhMVsGsqfJgt0lQk4jkR806ZzLm/GYlfHrrbXg9vAaOylk/RxoHLFCzbkS8YzNl3s4FxXjAnumwcZZ1gX95EabQNjcArum57BuafqgySC3kPKMzbkoOlxre0lrEvCB7TBcZl3T4RyV9ot4V9QzNueiaO80WDnX1rZ9dTW8e7dlbicqsn/vIjho0gObc1GUEuuaDpXb6bsA16WXhORCxHceeGBzLsq6GmD5Agtwj66C9hzVKY14xuZjbM5FWaZuqYrNjMZzFG1CPDEQhKjmL+yKSC9wEghJKZ5AmolOe6PUVohWe6PS1nmqOjObNxCR32L/vUH0qeqd2XzeVMhrYAMQkfWqujivH5qFKLU3Sm2FaLU3Sm11PsbmnCtCHticc0WnEIFtaQE+MxtRam+U2grRam+U2lry8j7G5pxzU827os65ouOBzTlXdPIW2ETkThHZISI7ReSxfH1uUCLSISLPish2EdkqIg+nH58uIs+IyGvpr02FbmuGiMRFZKOI/CZ9Hea2ThORJ0XklfS/8dvD2l4R+VL6Z2CLiPxMRKrC2lZ3ZnkJbCISB74LvBe4AviIiFyRj8+ehCTwFVW9HLgB+Fy6jY8By1R1EbAsfR0WDwPbJ1yHua3fBn6rqpcBV2PtDl17RaQN+AKwWFXfBMSBDxPCtrqzy1fGdj2wU1V3q+oY8HNgSZ4+OxBV7VbVF9P3B7FfvDasnY+nX/Y48P7CtPD1RKQdeB/wgwkPh7WtDcC7gB8CqOqYqh4npO3FthpWi0gZUAMcJLxtdWeQr8DWBuyfcN2VfiyURGQ+cC2wFmhV1W6w4Ae0FK5lr/PPwCPAxEoeYW3rQqAX+HG66/wDEaklhO1V1QPAPwKdQDfQr6q/J4RtdWeXr8B2ph21oVxnIiJ1wC+BL6pqjo5KyC0RuRvoUdUNhW5LQGXAdcD3VPVabL9wKLty6bGzJcACYA5QKyL3FbZVbrLyFdi6gI4J1+1Yeh8qIlKOBbWfquqv0g8fFpHZ6ednAz2Fat8ENwJ/ISJ7sW79bSLyb4SzrWD//7tUdW36+kks0IWxve8G9qhqr6omgF8B7yCcbXVnka/Atg5YJCILRKQCG4x9Kk+fHYiICDYGtF1V/2nCU08B96fv3w/8Ot9tO52qfl1V21V1PvZvuVxV7yOEbQVQ1UPAfhG5NP3Q7cA2wtneTuAGEalJ/0zcjo23hrGt7izytvNARO7CxoXiwI9U9e/z8sEBichNwApgM6fGrf4aG2d7ApiL/dDfo6o5LAeUHRG5Bfiqqt4tIjMIaVtF5BpsoqMC2A18AvvDGrr2isjfAR/CZso3Ap8G6ghhW92Z+ZYq51zR8Z0Hzrmi44HNOVd0PLA554qOBzbnXNHxwOacKzoe2JxzRccDm3Ou6Pw/MJUOkEJTXy8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(U,cmap='winter')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAXcElEQVR4nO3dfaxdVZnH8e+PYhF0iCCFlLZ3WifVsWBm1JuGqnGa6WRAZSyJYuqE0GjHZiYoYEy0lT/MxDQwGUPQyWBSebGMSO0AGRojAlYbneRCvVUzApWhWiyXVlqLjMxLWlqe+WPvmx5uz307+7zsvdbvkzT3nHX3uXuv3r2eve6znrOPIgIzM8vDaYM+ADMz6x8HfTOzjDjom5llxEHfzCwjDvpmZhk5fdAHMJ3zzjsvFi9ePOjDMDNrlN27d/82IuZNbK990F+8eDGjo6ODPgwzs0aR9Ot27U7vmJllxEHfzCwjDvpmZhlx0Dczy4iDvplZRhz0zcwyUvuSTbNBGxmBnTth5UpYseLVz6H94xUr+n+cZjPhoG9WahfM3/hGuP56OHYM5s6FW245+XzOHJDg+PFXPx7f7sgRXwysfhz0LWvjgb41uLcGcAleeaX4d+wY3Hdf8fXEiaINIOLVj48ehU9+smib6mLgC4ANgoO+ZaddoG8N7q0B/LTTTgbuuXPhQx+CH/1o6pm+dPKiMNnFwBcAGxQHfcvCdIG+NbhPNzt/29umzulP9VfD+MXAFwAbFNX94xKHh4fD996xTsw00LcLvFAtDz/d+sBMj8MXAOuUpN0RMXxKu4O+paRKoO9HgK378Vk6Jgv6Tu9YMkZGYNWqqVM3gw6kK1ac3Od4mmjiBWCyFNCOHQ78Vp2DvjXe+Ox5//6TlTV1CvSTmc0F4NgxuOsul35adQ761kiTlVqeXp7RdQ30k5nuAjBnDtx5p0s/rToHfWucydI4AJ/4BAwNNTsYtrsA7N8PX/taMfN32seqmDboS7oDuBw4FBEXl23/CPwVcAz4JfCxiHix/N5GYB1wArg2Ih4q298JfB04E/gOcF3UfRXZamUmaZyrr04rAI5fAEZGYMsWp32summrdyS9F/hv4K6WoP+XwPcj4rikfwCIiM9JWgbcAywHLgS+B7w5Ik5I2gVcBzxKEfS/EhEPTneArt4xePXsPtd3uU737mHP+q1Vx9U7EfFDSYsntD3c8vRR4MPl49XA1og4CuyTtBdYLukZ4OyIGCkP5i7gCmDaoG95aze7hzTSOLM1Xdrn2LGiDTzzt8l1I6f/ceBb5eMFFBeBcWNl28vl44ntbUlaD6wHGBoa6sIhWhNNnN23LtKmlsaZrXZpn7lzi78Cxv/PPPO3dioFfUk3AMeBu8eb2mwWU7S3FRGbgc1QpHeqHKM1j2f3M7diRRHYx2f2O3ee/D9zvt/a6TjoS1pLscC7qmVBdgxY1LLZQuBA2b6wTbvZq3h2P3utaR8o/q/alXl61m/QYdCXdBnwOeDPIuJ/W761HfimpJspFnKXArvKhdyXJF0CPAZcDfxTtUO3lHh23x2tM/+J+X7P+g1mVr1zD7ASOA94HvgCsBE4AzhSbvZoRPxtuf0NFHn+48D14xU6koY5WbL5IPCpmZRsunonfVNV5nh22jn/v+atSvXOR9s03z7F9puATW3aR4GLp9uf5ac1Dw2e3XfLVLP+nTv9f5srvyPXBqL11sMrV57MQzt3311TVfnceKMvrDly0Le+a007jKcaWitQHIS6r3XWP/Fzf53qyYuDvvVNu8Xa8VTDxo0OPL02Puu/8UaXdebMQd/6YqpSzPFPl7L+aE2nuawzPw761hderK0PL/DmzUHfeqr1JmFerK0PL/Dmy0Hfembigm0ud8NsEi/w5ue0QR+ApWvifWCOHPGCbR2tWFH8Xo4cOXWB3dLjoG9dNzJSpAjGUzpz5njBtgnGF3jHf1/jqZ6RkUEfmXWT0zvWVU7pNJdTPXnwTN+6yimdZnOqJ30O+tYVTumkxamedDm9Y5U5pZMep3rS5Zm+VeaUTpqc6kmTg751zCmdPDjVkxand6wjTunkw6metHimbx1xSicvTvWkw0HfOjLxT36ndPLg33vzOb1js9L6iVf+4JP8tKZ6xgO+b9DWLA76NmPtPvFq48ZBH5X1W+sdOieeDw789ef0js3YxDy+87l58/nQTA76Ni2XZlo7LuVsJqd3bEouzbTJuJSzmaad6Uu6Q9IhSY+3tJ0r6RFJT5dfz2n53kZJeyU9JenSlvZ3Svp5+b2vSFL3u2Pd5tJMm4pLOZtnJumdrwOXTWjbAOyIiKXAjvI5kpYBa4CLytfcKmlO+ZqvAuuBpeW/iT/TasglejYTPk+aY9qgHxE/BF6Y0Lwa2FI+3gJc0dK+NSKORsQ+YC+wXNJ84OyIGImIAO5qeY3V0HgeH4o/1b/4Rf/JbpMbT/WMnyfg/H5ddZrTvyAiDgJExEFJ55ftC4BHW7YbK9teLh9PbLcacmmmdcKlnM3Q7eqddnn6mKK9/Q+R1ksalTR6+PDhrh2czYxL8awKnz/11mnQf75M2VB+PVS2jwGLWrZbCBwo2xe2aW8rIjZHxHBEDM+bN6/DQ7ROOT9rVfj8qbdOg/52YG35eC3wQEv7GklnSFpCsWC7q0wFvSTpkrJq5+qW11hNOI9v3eD8fr1Nm9OXdA+wEjhP0hjwBeAmYJukdcB+4EqAiHhC0jbgSeA4cE1EnCh/1N9RVAKdCTxY/rOacB7fusn5/fqaNuhHxEcn+daqSbbfBGxq0z4KXDyro7O+aZeH9eC0qnxe1Y9vw2CA87DWGz6v6se3Ycicb5VsveRbMdePg37GnMe3fnB+v16c3smY66mtn3y+1YODfsacb7V+8vlWD07vZMh5fBuEifn98ZSPz7/+ctDPjPP4Nkjj+X1wjn9QnN7JjPOqVhc+FwfDQT8zzqtaXfhcHAyndzLhPL7VjWv4B8NBPwPO41tduYa//5zeyYBzp1Z3Pkf7x0E/A86dWt35HO0fp3cy0K4+2qxOfI72j4rPKa+v4eHhGB0dHfRhNI7f9GJN5vO3Okm7I2J4Yrtn+gnyopg1mc/f3nJOP0FeFLMm8/nbWw76CfKimDWZz9/ecnonIX4DlqXAb9rqLQf9RPgNWJYSv2mrd5zeSYTzoJYin9fd56CfCOdBLUU+r7vP6Z1E+M0tliKf193noN9wE9/E4kFhqZn4wSu+AFRTKehL+jTwN0AAPwc+BpwFfAtYDDwDfCQiflduvxFYB5wAro2Ih6rsP3de5LKc+Hzvjo5z+pIWANcCwxFxMTAHWANsAHZExFJgR/kcScvK718EXAbcKmlOtcPPmxe5LCc+37uj6kLu6cCZkk6nmOEfAFYDW8rvbwGuKB+vBrZGxNGI2AfsBZZX3H/WvMhlOfH53h0dp3ci4jlJXwL2A/8HPBwRD0u6ICIOltsclHR++ZIFwKMtP2KsbDuFpPXAeoChoaFODzF5XuSynPh8746Og76kcyhm70uAF4F/lXTVVC9p09b2Fp8RsRnYDMVdNjs9xlR58dZy5UXd6qos5P4FsC8iDgNIuh94F/C8pPnlLH8+cKjcfgxY1PL6hRTpIJsFL2aZeRxUUSWnvx+4RNJZkgSsAvYA24G15TZrgQfKx9uBNZLOkLQEWArsqrD/LHkxy8zjoIoqOf3HJN0L/AQ4DvyUIiXzemCbpHUUF4Yry+2fkLQNeLLc/pqIOFHx+LMzvpg1PsPxYpblyOOgc/7krAZyLtPM42A6/uSshvPirdmreVG3Mw76DeBFK7PJeXzMju+y2QBetDKbnMfH7DjoN4DfiWg2OY+P2XF6pwH8TkSzyXl8zI6rd8zMEjRZ9Y7TOzU2MlJ8IPTIyKCPxKw5PG6m5vROTbkiwWz2PG6m55l+TbkiwWz2PG6m56BfU65IMJs9j5vpOb1TU65IMJs9j5vpOejXjG+3YFaNb88wNQf9GvEilFn3eDy155x+jXgRyqx7PJ7ac9CvES9CmXWPx1N7Tu/UiBehzLrH46k934bBzCxBvg1Djflt42a953FWcHpnwFxhYNZ7HmcneaY/YK4wMOs9j7OTHPQHzBUGZr3ncXaS0zsD5goDs97zODvJ1TtmZgnqSfWOpDdIulfSLyTtkbRC0rmSHpH0dPn1nJbtN0raK+kpSZdW2XfTuZLAbHByHn9V0ztfBr4bER+WNBc4C/g8sCMibpK0AdgAfE7SMmANcBFwIfA9SW+OiBMVj6FxXElgNji5j7+OZ/qSzgbeC9wOEBHHIuJFYDWwpdxsC3BF+Xg1sDUijkbEPmAvsLzT/TeZKwnMBif38VclvfMm4DBwp6SfSrpN0uuACyLiIED59fxy+wXAsy2vHyvbTiFpvaRRSaOHDx+ucIj15EoCs8HJffxVSe+cDrwD+FREPCbpyxSpnMmoTVvbVeSI2AxshmIht8Ix1pIrCcwGJ/fxVyXojwFjEfFY+fxeiqD/vKT5EXFQ0nzgUMv2i1pevxA4UGH/jeYPSDEbnJzHX8fpnYj4DfCspLeUTauAJ4HtwNqybS3wQPl4O7BG0hmSlgBLgV2d7r+Jcq4YMKur3MZl1eqdTwF3l5U7vwI+RnEh2SZpHbAfuBIgIp6QtI3iwnAcuCanyp3cKwbM6ijHcVkp6EfEz4BTiv8pZv3ttt8EbKqyz6ZqVzGQ+sllVnc5jkvfe6dPcq8YMKujHMel773TJ7lXDJjVUY7j0vfeMTNLkD85a0Byqwwwa7IcxqvTOz2UY2WAWVPlMl490++h3O/xYdYkuYxXB/0eyrEywKypchmvTu/0UI6VAWZNlct4dfWOmVmCXL3TRzlUAJjlIMWx7PROl+VSAWCWulTHsmf6XZZLBYBZ6lIdyw76XZZLBYBZ6lIdy07vdFkuFQBmqUt1LLt6x8wsQa7eMTMzB/1uSbG0y8xOSmWMO6ffBamWdplZIaUx7pl+F6Ra2mVmhZTGuIN+F6Ra2mVmhZTGuNM7XZBqaZeZFVIa4y7ZNDNLkEs2eyCV1Xwzm50mj/3K6R1Jc4BR4LmIuFzSucC3gMXAM8BHIuJ35bYbgXXACeDaiHio6v4HJaXVfDObuaaP/W7M9K8D9rQ83wDsiIilwI7yOZKWAWuAi4DLgFvLC0YjpbSab2Yz1/SxXynoS1oIfAC4raV5NbClfLwFuKKlfWtEHI2IfcBeYHmV/Q9SSqv5ZjZzTR/7VdM7twCfBf6gpe2CiDgIEBEHJZ1fti8AHm3ZbqxsO4Wk9cB6gKGhoYqH2Bspreab2cw1fex3HPQlXQ4ciojdklbO5CVt2tqWDkXEZmAzFNU7nR5jr61Y0bxfuJlV1+SxX2Wm/27gg5LeD7wWOFvSN4DnJc0vZ/nzgUPl9mPAopbXLwQOVNi/mZnNUsc5/YjYGBELI2IxxQLt9yPiKmA7sLbcbC3wQPl4O7BG0hmSlgBLgV0dH/kANLlMy8x6o2lxoRfvyL0J2CZpHbAfuBIgIp6QtA14EjgOXBMRJ3qw/55oepmWmXVfE+NCV96cFRE7I+Ly8vGRiFgVEUvLry+0bLcpIv4oIt4SEQ92Y9/90vQyLTPrvibGBb8jd4aaXqZlZt3XxLjgG67NUNPLtMys+5oYF3zDNTOzBPmGa2Zm5qA/naaVY5nZ4DQhXjinP4UmlmOZ2WA0JV54pj+FJpZjmdlgNCVeOOhPoYnlWGY2GE2JF07vTKGJ5VhmNhhNiRcu2TQzS5BLNs3MzEG/nSaUXZlZvdU1jjinP0FTyq7MrL7qHEc805+gKWVXZlZfdY4jDvoTNKXsyszqq85xxOmdCZpSdmVm9VXnOOKSTTOzBLlk08zMHPTH1bW8ysyar07xxTl96l1eZWbNVrf44pk+9S6vMrNmq1t8cdCn3uVVZtZsdYsvTu9Q7/IqM2u2usUXl2yamSWo6yWbkhZJ+oGkPZKekHRd2X6upEckPV1+PaflNRsl7ZX0lKRLO923mZl1pkpO/zjwmYh4K3AJcI2kZcAGYEdELAV2lM8pv7cGuAi4DLhV0pwqB29mZrPTcdCPiIMR8ZPy8UvAHmABsBrYUm62BbiifLwa2BoRRyNiH7AXWN7p/ruhTrWzZpaHQcedrizkSloMvB14DLggIg5CcWGQdH652QLg0ZaXjZVt7X7eemA9wNDQUDcO8RR1q501s/TVIe5ULtmU9HrgPuD6iPj9VJu2aWu7ihwRmyNiOCKG582bV/UQ26pb7ayZpa8OcadS0Jf0GoqAf3dE3F82Py9pfvn9+cChsn0MWNTy8oXAgSr7r6JutbNmlr46xJ2O0zuSBNwO7ImIm1u+tR1YC9xUfn2gpf2bkm4GLgSWArs63X9VdaudNbP01SHudFynL+k9wI+AnwOvlM2fp8jrbwOGgP3AlRHxQvmaG4CPU1T+XB8RD063H9fpm5nN3mR1+h3P9CPi32mfpwdYNclrNgGbOt2nmZlVk929dwZdLmVmNm4Q8Sire+/UoVzKzAwGF4+ymunXoVzKzAwGF4+yCvp1KJcyM4PBxaOs0jt1KJcyM4PBxSPfWtnMLEFdv7WymZk1j4O+mVlGsgj6rs03s7rrV5xKfiHXtflmVnf9jFPJz/Rdm29mddfPOJV80HdtvpnVXT/jVPLpHdfmm1nd9TNOuU7fzCxBrtM3MzMHfTOznCQb9F2bb2ZN1cv4leRCrmvzzaypeh2/kpzpuzbfzJqq1/EryaDv2nwza6pex68k0zuuzTezpup1/HKdvplZglynb2Zm/Q/6ki6T9JSkvZI29Hv/ZmY562vQlzQH+GfgfcAy4KOSlvXzGMzMctbvmf5yYG9E/CoijgFbgdV9PgYzs2z1O+gvAJ5teT5Wtr2KpPWSRiWNHj58uG8HZ2aWun4HfbVpO6V8KCI2R8RwRAzPmzevD4dlZpaHftfpjwGLWp4vBA5M9YLdu3f/VtKvO9zfecBvO3xtU+XYZ8iz3zn2GfLsdyd9/sN2jX2t05d0OvCfwCrgOeDHwF9HxBM92t9ouzrVlOXYZ8iz3zn2GfLsdzf73NeZfkQcl/RJ4CFgDnBHrwK+mZmdqu+3YYiI7wDf6fd+zcws/Xfkbh70AQxAjn2GPPudY58hz353rc+1v/eOmZl1T+ozfTMza+Ggb2aWkSSDfi43dZO0SNIPJO2R9ISk68r2cyU9Iunp8us5gz7WbpM0R9JPJX27fJ5Dn98g6V5Jvyh/5ytS77ekT5fn9uOS7pH02hT7LOkOSYckPd7SNmk/JW0s49tTki6dzb6SC/qZ3dTtOPCZiHgrcAlwTdnXDcCOiFgK7Cifp+Y6YE/L8xz6/GXguxHxx8CfUPQ/2X5LWgBcCwxHxMUUZd5rSLPPXwcum9DWtp/lGF8DXFS+5tYy7s1IckGfjG7qFhEHI+In5eOXKILAAor+bik32wJcMZgj7A1JC4EPALe1NKfe57OB9wK3A0TEsYh4kcT7TVFWfmb5xs6zKN7Bn1yfI+KHwAsTmifr52pga0QcjYh9wF6KuDcjKQb9Gd3ULTWSFgNvBx4DLoiIg1BcGIDzB3dkPXEL8FnglZa21Pv8JuAwcGeZ1rpN0utIuN8R8RzwJWA/cBD4r4h4mIT7PMFk/awU41IM+jO6qVtKJL0euA+4PiJ+P+jj6SVJlwOHImL3oI+lz04H3gF8NSLeDvwPaaQ1JlXmsFcDS4ALgddJumqwR1ULlWJcikF/1jd1azJJr6EI+HdHxP1l8/OS5pffnw8cGtTx9cC7gQ9KeoYidffnkr5B2n2G4rwei4jHyuf3UlwEUu73XwD7IuJwRLwM3A+8i7T73GqyflaKcSkG/R8DSyUtkTSXYsFj+4CPqSckiSLHuycibm751nZgbfl4LfBAv4+tVyJiY0QsjIjFFL/b70fEVSTcZ4CI+A3wrKS3lE2rgCdJu9/7gUsknVWe66so1q1S7nOryfq5HVgj6QxJS4ClwK4Z/9SISO4f8H6Ku3n+Erhh0MfTw36+h+LPuv8Aflb+ez/wRorV/qfLr+cO+lh71P+VwLfLx8n3GfhTYLT8ff8bcE7q/Qb+HvgF8DjwL8AZKfYZuIdi3eJlipn8uqn6CdxQxrengPfNZl++DYOZWUZSTO+YmdkkHPTNzDLioG9mlhEHfTOzjDjom5llxEHfzCwjDvpmZhn5fxDYDO75+hNtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "b = np.ones([N,1])\n",
    "y = la.solve_triangular(U.T,b,lower=True)\n",
    "x = la.solve_triangular(U,y,lower=False)\n",
    "plt.plot(x,'b.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATYAAAD8CAYAAAD9uIjPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO2de3jU9ZX/Xye3yYUEEsIlXINIFfGKCFpa13qpQi34225b3addt+rSuvIrVrvVdd3aPnW32qrbunZ1+alb+7Rb67ZutbtaARVvrQgioogKgkIg3C8Juc/k/P4435EUA0wyk8n3Ozmv55knc/lm5iMmJ+f+FlXFcRwnl8jr7wM4juNkGjdsjuPkHG7YHMfJOdywOY6Tc7hhcxwn53DD5jhOzuGGzXGcPkNExorIsyKyVkTWiMiCbq4REblbRNaLyGoRmZru56Zl2ETkIhF5JzjQjekexnGcnCMOXK+qk4EzgWtE5IRDrpkFTApu84B70/3QXhs2EckHfhIc6gTgsm4O7DjOAEZV61V1ZXC/EVgLjD7ksrnAz9R4GRgiIjXpfG5BGt87HVivqhsAROTh4IBvHe4bpLpUGVMFGyqhMQalHRCLp3EEx3H+hKataNteSect5KJjlV3NqV38av0aoLXLMwtVdWG37ytSC5wGLDvkpdHA5i6P64Ln6lM7xEdJx7B1d5gZh14kIvMw9xLGDIGHb4Qbz4dlY6B2HwxrgoJOuzmOkx5LvpD+e+xqhhVfTe1a+U6rqk476mUig4DfANeqasOhL3fzLWnNeqaTY0vpMKq6UFWnqeo09o8xo3bxu/Dgb2H8PnhlNGyuAE3rj4zjOJlEU7ylgIgUYkbtF6r6aDeX1AFjuzweA2zt5cmB9Axbzw/TGDNPbXQDzNwMFW2wrxhaCiEh0OnGzXH6n+B3MZXb0d5JRIAHgLWqetdhLnsc+KugOnomsF9Vex2GQnqh6HJgkohMALYAlwJ/ecTvKO0wL+2hU+Gx4y0U/dcn4Q9j4fnxUNkC4/Z7WOo4/YmSyQhqJvBl4A0RWRU8dxMwDkBV7wOeAGYD64Fm4CvpfmivDZuqxkVkPvAUkA88qKprjvhNsTgMb4Klteap3fMEXLUStpfBz06BeJ4ZNsdx+pcMbTNT1RfpPm3V9RoFrsnMJxrpeGyo6hOYte3BJ3bC2AaobjZPbWcZJPLghhdhSwWsHgFFCRjVCIWJdI7nOE5viXjOO/uTBwWdFoJ+bDc8Vws3nwsdefC9Z+HC9bA+aAfp8KEIx+k3Mlg86A/S8th6jaiZ1MoWCz/rKuBXU2BrOVzyNjQVwo4yC1erm82Dcxwne7jH1ttPVsunnb4VVo2Er34WXh8Jdz4F1yw3Y/fGCGgq6rcjOs6ARIHOFG8hpX88tg8/PfiXKUxAcdw8tTdGQP0gOHEH7C82j25XKZS3+5SC42SLiHts/WvYkoxphBFNFn5eNQem7IDvLIX2fPjGRfB2NZy8HWoa+/ukjjMAkMgbtnBk6AsT1uMGsKfEPLX2fGsAHNpsLSIoHCiCjvx+ParjDAi8eJBBqpth+hYLP6+7EKpa4IrXYGgL/MuZ1h5y7B7vdXOcviSzDbr9QrgMW1ECilosp7Z2mHlqw5ph4h4o6zCD15FvN1GfUHCcviLi443hMmxJytstp4bCXWdZmDqjDv58LfzuY/DMBKg5YP1wEmJ/2HGiSsR/rcKRYzuUWNwKBeXt8OI4ePJYGN0Ic982L27bINgXs8H5iLvMjhM6kqFoKreQEk6PLUksYTm1eJ55astG20aQ25fAmmHwx7EwqB3G+uC842QU99j6kMKEbQMZt9/Czx+dabm26/8AU+vhvUrYUm6zpo7jZA732LKAqOXUSjrMU7tnuo1bffVV2F0C66qgQGHEAR+/cpxM4B5bFkgOzp+0w8LP6y6E7YPg9sXwxTXwfiW8MxTaomGnHSfUaOYWTfYX0bEEorb1bVC79bvtLoHFEy0U/bP3obnQ5kqbC2FIq49fOU46hDjMTIVoeGxJRK1QcMZWCz+vmmO5t9uXwI0vWv/bqpHQEOvvkzpOtPHJgyyTVLQqUGv3aC40r21fsRUZiuP2D76v2PrfPOfmOD0n4h5b9AxbkhEHLORsKoKrLzaj9s0/WAj67U+Z+tWUnSYc4zhOzwixN5YK0TVsRQm7NReafF9x3KqmFW2WhyvrsH1RrQWuW+o4PUEJdWEgFaJr2JIMaYXT6+1/xs3n2rTCJW/D3y6HB0+z9ePj9kHtfh+/cpxUiXgoGq3iQXfE4qYmX5SAlTWwYpSNX83YYp5bQ8y8NtctdZzU8eJBSCjtsJxaJ/DAafBfJ5hgzL/9r2mWLq213W61+2wtueM4hyHcUwWpEH2PLUlRwgoF1S0Wfv7nSbaV969XwYS9sGmwtYO41+Y4RyZVby3E/kHueGxJChOWU2stME9tc4X99bn5efhgMLxWA0VxW0fuuqWO0z3usYWMgk4rFEzaY+Hnd88xL+2WpXDuRtMt/WCIbQxxHKd7fKQqhCR1S4c2m7v8wWD4+cm2x+0La6CxCOrLgZKDhQfHcQ4S4jAzFXLXbclTKxRMrbfw8+qLbd34HYtg3qsWoq4ZZn1wjuMcJIOLJkXkQRHZISJvHub1c0Rkv4isCm7fzsR/Qm56bEny1G5FcShrN5WrlTXmuZ2y3VpBOvJdt9RxDiVzHttPgXuAnx3hmhdU9eKMfSK5btiSjGmEkU2wtRyunAOnboPvPWPe2oJZtvLolG22881xnIwZNlV9XkRqM/NuqZO7oWhXChM2bgXmpTXEzKjF82wFUtKguW6p4xjZ3aB7loi8LiJPisiUTLzhwPDYkgxrguntZrwWzDKjNu9VKzL8YKbpln5st+uWOgMb7VHFs1pEVnR5vFBVF/bg01YC41X1gIjMBn4LTOrB93fLwDJsycH5XaUWfjbEzKiN339QiT6RZyr0ea5b6gxgUvfGdqnqtF5/jGpDl/tPiMi/iUi1qu7q7XvCQDNsScrbLacG5qmVdsAnN8Flb8JvJsOzE0z+b4LrljoDlCz92IvISGC7qqqITMfSY7vTfd+BadhiccurHSiy8BPMqM1aB8+Nhx1lNkCfEFtH7sbNGWhkKH8mIr8EzsFC1jrgFqAQQFXvA/4CuFpE4kALcKmqpv0LNzANW5JYwnJqiTx4dLKNYA1tgR8ugtUjzOiVu26pMwDJXFX0sqO8fg/WDpJRjloVFZGxIvKsiKwVkTUisiB4vkpEFovIuuBrZaYP1+cUJqxQMLrBtBPumQ6DW+Hal60l5L0qaxEJ8eiI42Sc5KLJCI9UpdLuEQeuV9XJwJnANSJyAnAj8LSqTgKeDh5Hkzy1nNoxe00M5u4ZsL8YrnkFPrURtlTYdhBvBXEGChEXTD6qYVPVelVdGdxvBNYCo4G5wEPBZQ8Bl/TVIfucgk4rFEzZAS+Ng+svNHm/25bA59aa4vy7Q6HNDZszQBhIa4uCDuLTgGXACFWtBzN+IjL8MN8zD5gHQGlNGkftY5K6peVtMLzJWkKenAR1FXD+Bmvo3R+zAfoq1y11cplwe2OpkPLkgYgMAn4DXNu19+RoqOpCVZ2mqtOIhTwNJ2o5tzO2wDvVNn71wjj44WL41kuwYxC8PtKMm+PkKgNl0aSIFGJG7Req+mjw9HYRqQm8tRpgR18dMqskq5/5wdfmQlt7tLvUcnAVbfbXzHVLnVwm1z02ERHgAWCtqt7V5aXHgcuD+5cDj2X+eP3IyANwZp2Fn/M+C/9vqqnN/8vv7S/VK6NhZ1l/n9Jx+oaIV0VT8dhmAl8G3hCRVcFzNwG3AY+IyJXAJuDzfXPEfiI5ftVYZC0f5e3WHlIaaJdWtNl1LYXm5fmacSeXCHGYmQpHNWyq+iJwONN8XmaPE0KqWmHaVnPNbzrPDNrn1sL8V2Dh6dbUO36fVVUdJxdILpqMMANjbVE6xOK2BaQwYT1uq0Zaz9vp9ebFNRVBe4GtQAqxa+44PWIgFA8cDuqWgnlq5e0weSfc+z82tbC01gyg65Y6uYB7bAOEpG5pVYuFn49MsQLDl1bb2qPNQeW0M/o9QI7jHttAo6DTcmrtBeapvVdlHtotS2FjJayogeI4jG3wwXknmvRs0WQocY+tpxQmrFAwMdAtvfVsa+z9x+fh7A9g/VDYNMTnSp1oE/FZUffYekueWk5NME/tp6fC9jK47A3b87a5AvJct9SJKCEOM1PBPbbektQtPa3ews9rZpvK/J2L4IrXLOfmuqVOVHGPbQCTp+axFcetStoQs4mEreXWDtJYZPoJO8us/80H552oEHGPzQ1buohaoWBkk4WfV8y1ht5bn7FxrAWzzJM7dRuMcN1SJwIkF01GGDdsmaCg0255JdBUaJ7a/hi0FdgKpNbgn7kxZl6b59ycUBPuMDMV3LBlkqRuaXugWzq8Ca5eDkNa4bZPwMtj4LjdpqHgOGHGQ1HnQ5KD8zvLYF2VeWqVrdbYW9phObm4uG6pE37csDkfoaINTgt0S7//CTNq526EL6+GX02x/rfRjVZVdWk/J2z4ELzTLbG4FQpKOyz8fK7WDNkF75m8365S63VLRD+X4eQoPlLlHJZY3HJqcTFPbckxZvDuegpW1tja8SFtrlvqhI/OaPs80T592ClKmNEadcDCz38/3QoJ81+xTSEbK6F+kHttTvhwj805KvmdFooObjNP7c6PW1vItS/DjjJYM9xmUGsO+CZep//xHJuTEgWdViiYvNPCzxvOtymFf34aPvuu6Zaur3LdUic8uMfmpERSt3Rwm+1x21EGvzvOxq8uWm+bePe5bqkTBqJf1HKPLZt8qFu61cLPK+dY1fSORXDdH2FbOax23VInBGTIYxORB0Vkh4i8eZjXRUTuFpH1IrJaRKZm4vhu2LJNUtGqoBPy1bZ/vFcFu0tg0m44fpf9tdxbYo28jpNtkrOimZHf+ylw0RFenwVMCm7zgHvTPT54KNp/jGqEoc2wpxiummNG7aYXbCLhhgusyHDSdrvOcbJNhkJRVX1eRGqPcMlc4GeqqsDLIjIkKcSezue6x9ZfFCZgUDsglm/bU2JGLZaAwa1Q2WLXNRf6Nl4n+6QeilaLyIout3k9/KTRwOYuj+uC59LCPbb+pqoFTg90S2+4wIzapW9a68dPzoAXxsMxe01nwXGyReoe2y5VnZbGJ3X3QWnXW92w9TexuN32llj4WdliRu2k7VDWYYP07fmmWwo+oeBkh+y1ctQBY7s8HgNsTfdN3bCFhdIOM2Zgnlp5O5yyDS5cD4sn2uTCsCbXLXX6nuwumnwcmC8iDwMzgP3p5tfADVt4iMWtUNBcaOFna4EZtS+ugbeGQV2FeWvj9pvz7ltBnD4jc31sIvJL4BwsF1cH3AIUAqjqfcATwGxgPdAMfCUTn+uGLWwUdkLtXisYLJpoRi2WgO8+a20hy0dBSYfrljp9S4b+bqrqZUd5XYFrMvNpB/GqaNhI6pYesxeenQDf/6QN09/8PMzcZKNXmwcfzLk5Tl/gKlVOnzGsybyydVVw/1TbzPtXr9uc6QdDAIURrlvqZJiQz4Gmgv/ZDyvJwflTtpmk3/zZsHGIjV9dvgreHwJvD4MW1y11+gD32Jw+I6lbWhK3PW4NMXhpLGypgBl1toW3zXVLnT7A5fecPkXUllXWNFr4ecUlZtRuW2Irxq/5jHlyU+thuBs2JwPkQCjqhi0KJHVLUWsDOVBkRq250AxeZxAWNMRMld5zbk66hDjMTIWUDZuI5AMrgC2qerGIVAG/AmqB94EvqOrevjikEzCiCSraLfycP9smFL6+zJp5/+mTsGyMbQdx3VInXSLusfWkeLAAWNvl8Y3A06o6CXg6eOz0JUXBgLwAGyqtabei7aAiVlHCvLe2Am8HcdIgxcJBiL26lH76RWQM8Bng/i5PzwUeCu4/BFyS2aM5h6WizXJqVS1w69lw4/nw6ffg/sdh4h6rom4eHOofPCfEZHYfW7+Q6p/1HwHfArq2uo9IznQFX4d3940iMi+50oQ2j1QzQiwOw5tsAmHZGHhpnInFfGqjVU/3lFj+LRHuHz4nxERc8+Cohk1ELgZ2qOqrvfkAVV2oqtNUdRqxyt68hXM4iuOWU6vdB788Eb51gc2b/vhJq5yuHmH9bh6WOj0l1w0bMBOYIyLvAw8D54rIz4HtIlIDEHzd0WendLonqVs68oCNXz0wFaqb4WsrzOBtrITtrlvq9IJcz7Gp6t+r6hhVrQUuBZ5R1S9h60YuDy67HHisz07pHJn8ThjTAMcG+bUfzLSJhG++BOduhPcHWw+cb+J1UmUAeGyH4zbgAhFZB1wQPHb6g+T41XG7bOXRTedBcwHc+izMXgcbqky71HVLnVRICiZH2GPrUYOuqi4Flgb3dwPnZf5ITq8QtT9Tg1utv21bOTw6GbaUw2fetWLC7lJojFk11cevnCMR8aKTTx7kEnmBbunYBlg9HH57vBm1OxfB2mr4mzk2sXDGFjdszhEItzeWCm7Yco3k8smiQL+0udCM2o4yC1VHFkMiz1pCyjrcwDndE+L8WSq4YctVkrqlu0vNUztuF3z7OQsx/u7T8GYNnLzdZk0dpyvJHFuE8QanXOVD3VIs/NxfbEYtTy0PV91srx0o8mqp81EiXhV1jy3XqWqxnFoizzy1wa22hXd4E9w9A14cBxNdt9Q5hIh7bG7Ycp2kbumeEgs/q5vNqE3eZYPz7fnQnmdem6gLxDiGV0WdSFDWYTk1ME+ttAOmbYW578ATk+CZCbYlxHVLnZCHmanghm2gEItboeBAkYWfHflm1D73FqwcCVvL7ZpOcd1Sx0NRJ2LEEpZTa88zT23lSCiNw/eXwNvVti2ktMNmUD0sHbhE/O+aV0UHGoUJKxTU7rfw846Ztv7oWy/B9C2mW1pXYcUGZ4AS/UWT7rENVEQtpxaLm6d23zRrC7lyJewrNtX5vOAa11AYWCQXTUYY/7M8UEkOzp+83cLPBbMsz/aDxXDZm4FuabWJxzgDj4j3sblhG8jkKeSr5dSqWsxTW1prBu7jmy00bS2wcaw2N3ADigyGoiJykYi8IyLrReQj2igico6I7BeRVcHt2+ke339aBzpJ3dJRjRZ+XjkXZm6C2xebQbv6YvhgMJxeD8N8rnTAkCFvLFC3+wm22qwOWC4ij6vqW4dc+oKqXpyZT3XD5sBB3dI8tYbd5kIzao1FMLrBnldsLKukw3NuA4HMFQamA+tVdQOAiDyMCUEdatgyihs25yAjDtjIVUuheWqjG+C6P1qo+p1zbDvv8TttLZKTu/Qsf1YtIiu6PF6oqgu7PB4NbO7yuA6Y0c37nCUirwNbgW+q6prUD/xR3LA5BylK2K21wMLPvCD/Vt1sA/XFcUDs9Q/V6Z2cJPWq6C5VnXaE17t7o0PN5kpgvKoeEJHZwG+BSakeoDu8eOB8lMFtllOrbDFP7cbz4eJ34cHfWg/cK6Nhc0Wo+5icdMhoH1sdMLbL4zGYV/YhqtqgqgeC+08AhSJSnc5/gRs256PE4jCsyaYUXh1l7SCjG2DmZhNr3lds4arrluYumWv3WA5MEpEJIlKECUI93vUCERkpIhLcn47Zpd3pHN9DUefwlHRYTk0FHjoVHjveet/+9Un4w1h4frx5deN8/CqnyOCiSVWNi8h84CkgH3hQVdeIyNeC1+8D/gK4WkTiQAtwqaqmVZd1w+YcnqKEFQpaC6y/bV8x3PMEXLUStpfBz04xMeZx+/v7pE6myWDzbRBePnHIc/d1uX8PcE/mPtENm5MKBZ1m4KqbzVPbWWazpDe8CFsqTHG+KGG9cIXeCpITRDx/6jk25+gkx68+thueq4Wbz4WOPPjes3Dhehuc31Bpzzm5QWeKt5DiHpuTGknd0soWCz/rKuBXU2z86pK3oSlo6t1XbJ6dN/FGFxdzcQYUSd3S07fCqpHw1c/C6yPhzqfgmuVm7N4YAU1F/X1SJ118bZEzoEhWPwsT1rDbVGjGrH4QnLjDxq7iebYCqbzddUujSog3d6SCGzand4xphBFNFn5eNQem7IDvLLVZ029cZCuPXLc0ooTbG0sFD0Wd3lGYsHErMAWs/cVm1DrFhJqHNwHquqVRJLloMpVbSHGPzUmP6mbb2xbPg+sutL1uV7wGQ1vgX8609pBj93ivW9TwUNQZ0BQloKjFcmprh5mnNqwZJu4xyb94oFnquqXRIuKhqBs2JzOUtwe6pQp3nWVh6ow6+PO18LuPmXBMTaBb6tJ+4Sfi/4s8x+ZkhqRuaXm76ZY+eSyMboS5b5sXt20Q7IvZ4HzEvYEBgbd7OE4XYgnLqcXzzFNbNto2gty+BNYMgz+Otd1urlsaXkIu1JIK7rE5mSWpWzou0C390ZmWa7v+DzC1Ht6rhC3lrlsadrwq6jjdIGo5tZIO89TumW7jVl99FXaXwLoqKHDd0nAS7jAzFVL6sykiQ0Tk1yLytoisFZGzRKRKRBaLyLrga2VfH9aJEMnB+ZN2WPh53YWwfZCpX31xDbxfCe8MdVm/sDJAdEV/DPxeVY8HTgHWAjcCT6vqJODp4LHjHEQC3dJB7dbvtrsEFk+0UPTP3ocz60wRa/sgN3BhIjkEH+HiwVENm4hUAGcDDwCoaruq7sMktB4KLnsIuKSvDulEmKRu6RlbLfy8ao7l3m5fAje+aP1vq0ZCQ6y/T+p0JeIeWyp/Jo8BdgL/ISKnAK8CC4ARqloPoKr1IjK8u28WkXnAPABKazJxZidqJBWtCtTaPZoLzWvbV2xFhuK4/ZLsK7b+N8+59T8h9sZSIZVQtACYCtyrqqcBTfQg7FTVhao6TVWnEfM03IBmxAGYscXWGl19Mdx7BnzzD/DjJ20l0iujbTuv0/9EvCqaimGrA+pUdVnw+NeYodsuIjUAwdcdfXNEJ2coSkB5m93fXGGbQUo6YEir5eHKOmwra2uB9cE5/UOqYWiIQ9Gj/vSo6jZgs4gcFzx1HiZP/zhwefDc5cBjfXJCJ/cY0mq6peVttmb8H86zLbwPPGbaCsvGuG5pfxPx4kGqpaj/C/wi0AXcAHwFM4qPiMiVwCbg831zRCfniMVhWNxyaq+MNk/tb5fblpBHplghobXA8nF5WJjqZJeI/5OnZNhUdRXQnYz9eZk9jjOgKO2AKTst/HzgNPivE0ww5t/+1zRLl9babrfafW7cskq4vbFU8ESG038UJUxhvrrF1K/+8yTbyvvXq2DCXtg02NpBQpykzkl80aTjZIDCBIzbZ+Hn8+MP5tdufh4+GAyv1UBR3NaRu25pdoi4g+wem9P/FHRC7X6YtMfCz++eY97ALUvh3I2mW/rBEK+UZpMBUjxwnL4lqVs6tNm8hQ8Gw89Ptj1uX1gDjUVQXw6UwLAmb+Lta9xjc5wMkadWKJhab+Hn1RfbuvE7FsG8Vy1EXTPMJhecviWDHpuIXCQi74jIehH5SHO/GHcHr68WkanpHt89Nidc5KndiuJQ1m4qVytrzHM7Zbu1gnTku25pX5LB5lsRyQd+AlyANfsvF5HHVfWtLpfNAiYFtxnAvcHXXuOGzQknYxphZBNsLYcr58Cp2+B7z5i3tmCWrTw6ZZvtfHMyT+YqntOB9aq6AUBEHsYWaHQ1bHOBn6mqAi8Ha9JqkrPovcFDUSecFCZs3ArMS2uImVGL59kKpKRBc93SPiDFMNRC0WoRWdHlNu+QNxsNbO7yuC54rqfX9Aj32JxwM6wJpreb8Vowy4zavFetyPCDmaZb+rHdrluaaVIPRXepanfN+0m6c/0OffdUrukRbticcFOUsNuuUgs/G2Jm1MbvP6hEn8gzFfo81y3NCMlFk5mhDhjb5fEYYGsvrukRbticaFDebjk1ME+ttAM+uQkuexN+MxmenWDyfxNctzQjZO6fcDkwSUQmAFuAS4G/POSax4H5Qf5tBrA/nfwauGFzokIsbnm1A0UWfoIZtVnr4LnxtgJpULsNzufjxi1dMlQ8UNW4iMwHnsL+zzyoqmtE5GvB6/cBTwCzgfVAM7ZkIy3csDnRIpawnFoiDx6dbCNYQ1vgh4tg9QgzeuWuW5oWmQ1FUdUnMOPV9bn7utxX4JqMfSBeFXWiRmHCCgWjG0w74Z7pMLgVrn3ZWkLeq7IWkRAPaEeCiC+adI/NiSZ5ajm1Qe0mBnP3DNhfDNe8EhQaqiG/08JXH5zvOSGeA00F99icaFLQaYWCKTvgpXFw/YUm73fbEvjcWlOcf3cotHmPW69wj81x+glRS0eXt8HwJvPUnpwEdRVw/gZr6N0fswH6qlYfv0qZcG/uSAU3bE60EbWc25gGCz+vnGNG7YeLbUPI38yB+kEwbasbtlRJLpqMMG7YnOiTrH7mB1+bC82o7S6FY/ZCRZt5IK5bmjohDjNTwQ2bkzuMPABVLRZ+zvusGbUbXzSDd/O5JhwzZadVVJ0j46Go44SE5PhVY5G1fJS3W0W0tMO8topA07Sl0Lw8r5YeHvfYHCdkVLVaTk0FbjrPDNrn1sL8V2Dh6dbUO36fVVWdj5LhBt3+wNs9nNwjFrctIIUJ63FbNdJ63k6vNy+uqQjaA7X5iCfJ+wxv93CckJLULQXz1MrbYfJOuPd/bGphaa0ZQNct/SgRN/jusTm5S1K3tKrFws9HpliB4Uurbe3R5qBy2hn9vq3M0qNFk6HEPTYn9ynotJxae4F5au9VmYd2y1LYWAkraqA4DmMbfHAeQh9mpoJ7bE7uU5iwQsHEQLf01rOtsfcfn4ezP4D1Q2HTEF8x3hX32BwnIuSp5dQE89R+eipsL4PL3rA9b5srIM91SwH32BwnMiR1S0+rt/DzmtmmMn/nIrjiNcu5uW6p0Smp3UKKe2zOwCJPzWMrjluVtCFmEwlby60dpLHI9BN2lln/20CcL82BPjY3bM7AQ9QKBSObLPy8Yq419N76jI1jLZhlntyp22DEANUtjXgo6obNGZgUdNotrwSaCs1T2x+DtgJbgdQa/Go0xsxrG2g5N/fYHCfCJHVL2wPd0uFNcPVyGNIKt30CXh4Dx+02DYWBhHtsjhNhkoPzO8tgXZV5apWt1thb2mE5ubgMPN3SgWDYROQbwEPDxJ8AAAvkSURBVFXYf+4bmDxWKfAroBZ4H/iCqu7tk1M6Tl9T0QanBbql3/+EGbVzN8KXV8Ovplj/2+hGq6rmurSfCnRGu2HiqKcXkdHA14Fpqnoitoz5UuBG4GlVnQQ8HTx2nGgSi1uhoLTDws/nas2QXfCeyfvtKrVet0S4G1MzxgAZgi8ASkSkA/PUtgJ/D5wTvP4QsBS4IcPnc5zsEotbTi0u5qktOcYM3l1PwcoaeGEcDGnLfd3SiBvvo3psqroFuAPYBNRj8vOLgBFJGfrg6/Duvl9E5onIChFZQZtHqk7IKUqY0Rp1wMLPfz/dCgnzX7FNIRsrTUMh4r/4RyXiHlsqoWglMBeYAIwCykTkS6l+gKouVNVpqjqNWGXvT+o42SS/00LRiXvNU7vz49YWcu3LlnvbNNhuuThfmmzQjfCsaCoZwvOBjaq6U1U7gEeBjwPbRaQGIPi6o++O6ThZpqDTCgWTd1r4ecP5NqXwz0/DZ9813dL1VbmrWxpxjy2VHNsm4EwRKQVagPOAFUATcDlwW/D1sb46pOP0C0nd0sFttsdtRxn87jgbv7povW3i3ZejuqVZmAMVkSpS6KwQkfeBRiABxFV12tHeO5Uc2zLg18BKrNUjD1iIGbQLRGQdcEHw2HFyi6Ru6RlbYc1w0y19eQzcsQiu+yNsK4fVI8245QxZWzTZk86KT6nqqakYNUixKqqqtwC3HPJ0G+a9OU5uk6x+FnRCvtr2j/eqYHcJTNptG3pVYG8JlLVHf/wqe2HmXPqos8InDxwnVUY1wtBm2FMMV80xo3bTCzaRcMMFVmQ4abtdF3WyUxj4k84KEem2swIzs4tERIF/V9WFR3tjN2yOkyqFCbs1xizfNrTZjFosAYNbobLFrmsuhMKI65am7rFVi8iKLo8XdjU8IrIEGNnN9/1DD04zU1W3BoZvsYi8rarPH+kb3LA5Tk+paoHTA93SGy4wo3bpm1BzAH5yBrww3lTox0dYtzT14sGuI+W9VPX8w70mIttFpCbw1g7bWaGqW4OvO0Tkv4HpwBENW7QHwhynP0jqlhZ0wuoR8OZwM2onbYeyDhukb8833dJ4BH/FstfH9jjWUQGH6awQkTIRKU/eBz4NvHm0N3aPzXF6S2mHGTMwT628HU7ZBheuh8UTbXJhWFM0dUuzc9zbgEdE5EqsrezzACIyCrhfVWcDI4D/FhEwe/Wfqvr7o72xGzbH6S2xuBUKmgst/GwtMKP2xTXw1jCoqzCvbtx+W0cepa0gWSgeqOpuuumsCELP2cH9DcApPX1vN2yOky6FnVC718arFk00oxZLwHeftbaQ5aOgpCNauqURssHdEcEEgOOEjKRu6TF74dkJ8P1PWi/bzc/DzE02erV5cITyba4E7zhOV4Y1mVe2rgrun2qbef/qdZsz/WAIoDAi5LqlSqil9VIhKn9CHCf8JAfnT9lmkn7zZ8PGITZ+dfkqeH8IvD0MWiKgWzoAhuAdx0mVpG5pSdz2uDXE4KWxsKUCZtTZFt62COiWhjjMTAU3bI6TaURtWWVNo4WfV1xiRu22JbZi/JrPmCc3tR6Gh9Ww9fcB0sMNm+P0BUndUtTaQA4UmVFrLjSD1xkk3xtipkofppybK8E7jnNERjRBRbuFn/Nn24TC15dZM+8/fRKWjYHjd4VPt9Q9NsdxDktX3dINlZDIs9za8CabXChKmPfWVmDryEPR5yaRr4q6YXOcbFDRZjk1Fbj1bDNqn34PvrIKfnmi9b+NaQiPbmnEQ1Fv93CcbBCLm5dW0mHh50vjTCzmUxuterqnxPJviRB4S6m2eoTA/h4O99gcJ5sUxy2n1inmqT15rM2b/vhJ6317frwZunH9rFvqHpvjOCmT1C0decDCzwem2gqkr60wg7exEraHQLfUPTbHcXpMfqfl1KpazFNriNlEwjdfCgRihkNRp3lz/bGJt7/D4TRxj81x+oPk+NVxu2zl0U3nQXMB3PoszF4HG6pMu7Q/dEtzQDDZPTbH6S9EzbUY3Gr9bdvK4dHJsKUcPvOuFRN2l5rGQlVLdsevQhxmpoIbNsfpT/IC3dKxDRZ+/vZ4M2p3LoK11fA3c2xi4YwtWTRs4fbGUsENm+P0N8nqZ1GgbNVcaEZtR5mFqiOLrbF3T4lpKmTDwLnH5jhORkjqlu4uNU/tuF3w7ecskf93n4Y3a+Dk7TZr2tdE3GPz4oHjhIXCBAxqt/u7SmF/sRm1PLU8XHWzvXagyNaQ9xXJRZOp3EKKe2yOEzaqWiynlsgzT21wq23hHd4Ed8+AF8fBxD7WLfVQ1HGcjBKL221PiYWf1c1m1CbvshnT9nxozzOvTbRvJhQiHoq6YXOcsFLWYTk1ME+ttAOmbYW578ATk+CZCTDiQN/olrrH5jhOnxCLW6HgQJGFnx35ZtQ+9xasHAlby+2aTsmsbqkvmnQcp8+JJSyn1p5nntrKkVAah+8vgberbVtIaUdml1VG3GPzqqjjhJ3ChBUKavdb+HnHTFt/9K2XYPoW0y2tq7BiQ0ZIsSLqVVHHcdJG1HJqsbh5avdNs7aQK1fCvmJTnW/L0K90xENR99gcJyokB+dP3m7h54JZlmf7wWK47E3TLW3OgGZpDiyadMPmOFEiTyFfLadW1WKe2tJaM3Af32wryDNBFrZ7iMjnRWSNiHSKyLQjXHeRiLwjIutF5MZU3tsNm+NEjaRu6fQtFn5eOReemgi3L4YJezPzGdnx2N4E/hx4/nAXiEg+8BNgFnACcJmInHC0N/Ycm+NEkaRuaTzP8mytBebBZapZNwuFAVVdCyByxM+aDqxX1Q3BtQ8Dc4G3jvRN2TVse9/axX+d2ATsyurnpkc10TlvlM4K0TpvuM/6JDAMgPHpv9mrT4FUp3hxsYis6PJ4oaouTP8MHzIa2NzlcR0w42jflFXDpqrDRGSFqh42ng4bUTpvlM4K0TpvlM6aLqp6UabeS0SWACO7eekfVPWxVN6im+eOGgR7KOo4Tp+hquen+RZ1wNguj8cAW4/2TV48cBwnzCwHJonIBBEpAi4FHj/aN/WHYctk/J0NonTeKJ0VonXeKJ01EojI/xGROuAs4H9F5Kng+VEi8gSAqsaB+cBTwFrgEVVdc9T3Vg1xl53jOE4v8FDUcZycww2b4zg5R9YMW2/GIrKJiIwVkWdFZG0w5rEgeL5KRBaLyLrga2V/nzWJiOSLyGsi8j/B4zCfdYiI/FpE3g7+jc8K63lF5BvBz8CbIvJLESkO61md7smKYevtWESWiQPXq+pk4EzgmuCMNwJPq+ok4OngcVhYgCVUk4T5rD8Gfq+qxwOnYOcO3XlFZDTwdWCaqp4I5GOVuNCd1Tk82fLYPhyLUNV2IDkWERpUtV5VVwb3G7FfvNHYOR8KLnsIuKR/TviniMgY4DPA/V2eDutZK4CzgQcAVLVdVfcR0vNi/Z0lIlIAlGJ9U2E9q9MN2TJs3Y1FjM7SZ/cYEakFTgOWASNUtR7M+AHD++9kf8KPgG8BXYcDw3rWY4CdwH8EofP9IlJGCM+rqluAO4BNQD2wX1UXEcKzOocnW4atV2MR/YGIDAJ+A1yrqg39fZ7uEJGLgR2q+mp/nyVFCoCpwL2qehrQREhDuSB3NheYAIwCykTkS/17KqenZMuw9WosItuISCFm1H6hqo8GT28XkZrg9RpgR3+drwszgTki8j4W1p8rIj8nnGcF+/9fp6rLgse/xgxdGM97PrBRVXeqagfwKPBxwnlW5zBky7D1aiwim4jtTnkAWKuqd3V56XHg8uD+5UAqg7t9iqr+vaqOUdVa7N/yGVX9EiE8K4CqbgM2i8hxwVPnYWtnwnjeTcCZIlIa/Eych+Vbw3hW5zBkbfJARGZjeaF84EFV/aesfHCKiMgngBeANziYt7oJy7M9AozDfug/r6p7+uWQ3SAi5wDfVNWLRWQoIT2riJyKFTqKgA3AV7A/rKE7r4h8F/giVil/DbgKGEQIz+p0j49UOY6Tc/jkgeM4OYcbNsdxcg43bI7j5Bxu2BzHyTncsDmOk3O4YXMcJ+dww+Y4Ts7x/wHDEp6jmIeGYQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(A,cmap='winter')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "80000"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.nbytes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "80000"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "8*100*100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2384"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "8*(100 + 2*99)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}