math_linear_algebra.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use the scipy.linalg module, for solving linear systems of equations\n",
    "# use the scipy.optimize module, for solving nonlinear equations\n",
    "\n",
    "from scipy import linalg as la\n",
    "from scipy import optimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "sympy.init_printing()\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# When solving a system of linear equations,it is important to look at the \n",
    "# condition number to estimate the accuracy of the solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAA0AAAASCAYAAACAa1QyAAAABHNCSVQICAgIfAhkiAAAAMJJREFUKJHF0j9LQlEYB+DH6BNYhJ/HxcmpD9CgQxCImxAIbi5u7eLSJji5NzQGDRItDtHQ3J0aSnS4LxJytBMN/eBw4XCe+77nD3/MCdqYYYkPFLhHC0cpdIk13nCLIcZ4j/kpKruojmbijzW8Bjz/TevXgW7s6zORz/h+5VY5xiIqNXLRKMA8F3QCPKOaA64CPClP8Md0AyxwlgN6AR5xmgP6AR4c2MP3Z3GBCVbKSywS619izTaDqHJo3OW0+8/ZALBUMPfsVQVVAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle 2$"
      ],
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGMAAAA2CAYAAADTeCfRAAAABHNCSVQICAgIfAhkiAAABMxJREFUeJztm12IVVUUx39+jYijhFPmRzQgMaA00AeWKOrNStTyA9EmIoRMHxrzzRd98QqigiiTZA8JcR+DCFEf/MAXC6NBBKOH0oKsB8sPpCRE09CHde547nafs9e5d597z1z2Dw7M2WedvdfZ6+y91z73PxAoDCMy2k8E/snDkUB2NgG9rXaiXRmZ0f554Mc8HAlkC0YncCsvRwLZgrEEOJWXI4FsHADGtNqJdkY7MkYhmdc9T+2+Cjxo8VE4Rivt5gLfeWx3NtnT6rZHOzKWAsc9tTkKfyOsrdAG4wngb09tzsPvKGsbzGBsBl4yynqAix7b7CXsVazEg/E2sAF4x7BZDhy13LsVOIfsPa4Dx5BNoQvb4rkguv9KdH2NxeYy9oX4oGHXD/wG3AHOA/MVPtWLqw+0Plt5HfjJKEu68STwQdR4L3AY+AuYlFJ/T9SGyVJgJ7Ca5GA8BUyJHW9EtqWYTR+yHm0EZiLp+L/Asyk+mVSAstLW1QcanxMZDdxEOg2gC9ihdKwT+B8ZSVXWAW/Fzjfi3qskBcNkAPiV2qxsEDhk2P0C7FbUV6WCPhgmtj6IY/N5CHPNuI9kTSuj82Xos6gJUX03o/OxwDXgvZjNGPxkUh3A+8AXPJr2OoCXefwrwSkkNW8GZh/Esflcgy2bOgKsiP5+BXnbNAwAF4Dvo/O7wAmgGxlxXcANZV0uViEZXiVW9iSSNl81bK8iU0QzMPsgjs1nJxOA28B0YL/ynr3An8Bzlmu7gUXAu5EzLjTT1ElksYwzLbrXXLC3Az+n1LUNWVeqxz3gP6NMkwSk9UGSzypOAF/yaISksQ95+2YlXJ8DfAJ8pGzbFYxuZF5eaZR3INPsWqP8IHAmpb5JSAdWj6+RhT9eNs7hs6sPknxW0Y+MDpcTBxxOgCxWg8B6ZduuYJSRN9D2KWcQ+Nwou0S+C7imD8ok++zkGWTtSOMzJL9eRG361mmxPQTMSKmrE3ghOh4AW6K/zZR0JPA7sCehnj5kitmApLYDyDTT7XiWOBX0wdD0gctnFeMd15O+hpYttl2OukoJdVUMu8VReQ/J9CObrbvIpm+Bo22TCvpgaPpA43MgEAgEApmofiMp5M+QgUAgELS2w5qgtc2RoLUtEEFrWyCC1naYErS2OdMqrW2RKIzut1Va2yJRGN1vK7S2RaJQut9WaG2LRKF0v41obTXyzss0IG/MQL1SU5vut2VS00a0tiXkN+C5yG/A94HT1Mo7ZwNTY8ebUflXaU7VgcYXG7av1eOBH4CPU+7TPFcfoorZBbwInEWmerXUNIvW1sQlbQSHvNEjGl+SdL9xmio1NaepM8DT1GpttSrANGkjKOSNHrH5Yup+FwLfeGjLm9TUp9Y2TdoIdcob68T0JU/drzepqS+t7V5k4VuLTA82PkQCeyWlnjLuDVKpDl/y1P2mPZc5A4ywlA1h2/QdR6I8HXlzXFPKPmSYvobMmza6kf9NWO2o61NEVprGHw348i0SqMlIcBol6bluIC+COQom8/hocaLV2mqkjdCgvFGJxhffut8y+UpNAZ3WVivv9CJvdKD1RaP7LZLUFNBpbbXyzmbIG7NITV2631JCXRXDrhlS0yFcWtvhikv3GwgEAoFA4XkIMLfCFiPzP/MAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\frac{\\sqrt{2 \\sqrt{170} + 27}}{\\sqrt{27 - 2 \\sqrt{170}}}$"
      ],
      "text/plain": [
       "  _____________\n",
       "╲╱ 2⋅√170 + 27 \n",
       "───────────────\n",
       "  _____________\n",
       "╲╱ 27 - 2⋅√170 "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMYAAAASCAYAAAAANKFkAAAABHNCSVQICAgIfAhkiAAABgNJREFUaIHt2muMnkUVB/AfWEJBoFVBNygRWpGLJEQwKCqlLdKUcBfhEyDGCkTBghAgKLqYkBREAxIvEG4BDB8gAl6432+SaBRQAbn1hRYo0DYWkApilw9n3uyzz77Pdd/dfnn/yZvZnZlzZubMmZlzeRhggAFKcQxGKn7/b8CvU8JnRQHN/rgdy7EWz+M67Nmj70ewCDfg2dR/DR7EN7FhzXkelZnXooI+n8DleBnvpLVdgA8V9P8aLsIDeCPxvqZiHufiLiwTa1mNv+FHYq15HGPi+1Vn7W1pJltmHc31q/YY0zJ/P4qzCxjuhfm4pWSivbBGCCOPt3rUnYvTsAo3YiU+hYNxGI42dhGH41d4BffgRXwMX8Wl2C/1GSmZ3zZCUG9hs4I+s/EwPoqb8BT2wGIsxJfSnLP4AXZNfJdjx5I5dHEy/oo78Bo+iC9gGMemv5dl+k90v+qsvS3NVMmsiX61HaMUfxIKdlADmk761cGQuN1WCGFmMS+N/Xyufj4ONP5lGBKHZEQcqCJsgDvxHH6i+Aa8LbWdmKv/War/dQ+aedg+jTFXvRdjekH9OYn+lxX0WVTtV921t6WZCpl11NevtmOUYpfEYDk+0ICuo/7EP5/GuKmg/Q282WDsMxO/i0r6LMY6zBG3cq+NnpXqlxp/ADcXN89/xO1ehLkmtgG7Jvo7avavs1911t6WZqpk1tH8YNQeo44dflwqL9PMx4CNcaRQ1MXixPbarGfwrnhut8y1zRECvbPBuP9L5XsF7TthCS7E/SV85qfydqEUWbyJh7CpMHMmCwem8vGa/av2q+7a29JMpczq6ldjTKto3yQNvE7Y7U0xhKtzdUvxDdyXqVuN08VT+4TwMVYJW/UgcVsepx6mCX8Ebi1ov1qYW2dW8NohlU8XtD+DBfi0cJz7gVOF/T4Dn8OXxaFYUoO2ar+arL0tzVTKrK5+NUbVwTgCM/FHYx2/OrhCeP//FDfFLJwgHMlbRKTpsUz/C8TTeDm+lal/FlcKh7QOlghz4mZh6+bxQ3xWKNzaCl4zUrmmoL1bP7Pm3OrgVBFE6OJWEYF6vQZt1X41WXtbmqmSWVP9aoQqU+rYVF7cgvfZuBuv4m38A8eLV2ETYadmcRquF4dgtrBBdxdO929wXo0xv4tTRBTkqB7te4hb76fCQZ0oNkhlWeSrKYYS3yERYZslwra71aAt2682a++3vOifzJrqVyOUHYyd8UXhxN08kUFy6EYk5mTq5opw7e/wPXEY3hbhy0PxklD4WSV8vyNs4CeErbk61941CZ7GWTXn2r3dZhS0b5Hr10+8KnI0C0Qe46qK/mX71WbtbWhYvzKjt371FReKUz3cZ75bJL7/zdSdr3d4r4vfKg+/npTa/258uLeLmaoTYt1fNza+KP1f9GJ2w5L7FLTTh7CgeDFGjA9MZFG2X23W3oaG9S+zXvrVeIwiH2O6MEXWiehGP9HNYmfzEhuncqsCmm79uz3aThd+xaPYVyQGe+EdxWvZTdjRD+JfRs2Ge1K5QLyu2SjL5iJRtRaPFPDtF7ZOZVFUsGq/2qy9DQ3rX2a99Ktv6Kb9f1/Rb7bIHm6Uq/8MPtyj/ydFVGLE2AjHEUZT+R/P0ewnhLvW+E8jzkp0fykYry6G9TfBl8Vc1bffjsKnyGNDowm+h0ro6+5XLwxr9klIHZrJlllT/Wo8RtGL0XXiLqlgfleazHbGJlsOxxni9lgqogazxbdQ04UNfH6m//UiT/EVPCls6xUifn6AcNjOMPYzgq/jx+IWfUA43nl0hDM/EXxbfN7wc/H8PykSkvOE/f39HjSHpB+jCr9nZi4rRfSpi4Uim3y/yCyvEpGpvYVftcLYSF0edfdrqjDZMmuqX23GGIedxElapjpZ0kl9t83V741rRXTo3yLh9rrIRxxtNDKRxUbCV3hEZLrfEyHaP4hnOY9h1bbvvRXzz/MqugG3EeHBV4Q594Kw6Yteqaq5dXL9d8EvhDm4Uqx9Df6ceJW9hk32q2yu/XwxmFyZtdGvpmMMMMAAAwwwwAADDNAO7wMuFWmq1gECmQAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle 7.58240137440151$"
      ],
      "text/plain": [
       "7.58240137440151"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACwAAAAbCAYAAAAH+20UAAAABHNCSVQICAgIfAhkiAAAAmBJREFUWIXt1tuLzVEUB/CP4xKhMG5NSh6Uy5RLHph4EjLyIJd/A5FSRBkmlygNRfLkDQ8iI5JbKMo9hdxSSBhEbjPGw96T4+ec8/sdnZPSfOvX3q219netVvv3XZsuVBfdqszfHx+qnKOiWIQ5/7qIctCMXpUkzFWSLIHu8ftWSdJqFlyPS1XkrziaUFNp0mp2uAZvKk2aLHgzTuMZPuMtrmOd8ro1GvcyxM3AYbzA17ieREOxA0kd/oZruItX6IupmILncf8sQyErcBT3S8SswQa8xrFY7GBMwhmsypBH7yL2jejA7iwk2JPiXxz5TgnDJYmeGfMUxYS8BJ3IoRFzE7GDsKkEVw6P8AlDyi2kR8a4+XG9lWebjWmoRUuevQHHS3DVYxQOoRXzUIcvuILLGWv6DSuxHjtwQejuTX92ZKJw1/N/3l1Kq8/yyNcsNKAj8Z0rkCcVLxMkLRhWJPYJpsd9T+xM4W6KnG14gJnoh/E4EX1nyy24E8OwQJCo55hcIGYntsb9LCxM4dwSi2oX/o189BFUqEO4bn+NkYJO3ingm+mXfG0VulUKq2NBD4r490X/0kLOrJPuqaDN4wWtzEfnnRsr6PbHFK7OgfKuiL81rn0KOcsZzbVxbU/Y2wRVWIvbGXjOxzOjFX561sX1SRrRGAwvYM/5NTguFjm7JPpHpCWJOBDjGxP2WfghdH9AGskyfBfeEnuFv3k/HkbyFxhX5Gx/5T0lhwp3uEPo+DYcFDr/XZiEqagTNPSGMN/b8B5XBU0elHJ+YBkFi3zb8Vh4w7zBEeG90oUu/Lf4CX0Kich8YfSLAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle 3 \\sqrt{6}$"
      ],
      "text/plain": [
       "3⋅√6"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.matrix([[2, 3], [5, 4]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.array([4, 3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.matrix_rank(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANMAAAASCAYAAADBs+vIAAAABHNCSVQICAgIfAhkiAAABjRJREFUaIHt2mmMnlUVB/AfUEJRoHXDCUrQ1gWUhIAGBLG0FRsIgit8EsRYwQBaQQKEzcHEpCgYlgBqZIlo/CBRcQFkX5VEo6jILn21BQq0DWUra8cP5z6ZZ555n3Xemcbk/Sdv7sy999x7zzl3OcvDEEMMMXAcibGa3+stxutVjLO6hOYgXI9V2IBH8Qvs3afvW7AUv8Ijqf963IkvY/OG6zw8t66lJX3eicvwOF5OvJ2HN5X0/zwuxB14No3905p1nI2bsFLwsg5/w7cEr0Ucaer6asJ7V5rplllP+/3VRS9FlPI/K/f3PTirZICPYTGubTnxeiHAIp7vU3c2TsJa/Bpr8B58Cp/DESYyfiguwRO4Bf/F2/FZ/BgHpj5jFevbUQj3eWxT0mc+/ojtcTUewJ5YhgPw0bTmPE7HbmncVdi5Yg0ZjsdfcQOewhvxEYziqPT3ylz/qeqrCe9daWZKZm32V9c58ugis0n4k9iUh7Sg6aVfE4yIW3S1UEAei9LcjxbqF+Ngk1+gEXGwxsQhLMNmuBH/xveU37R/SG1fK9R/P9X/oA/NIrw3zbFQsxtwdkn9dxL9xTX0edTpqynvXWlmQmY9zfdX1zny6CKzSdg1Ea7CFi3oepozu1ea4+qS9mfxXIu5T03jXVjRZxk2YoG4/fsJZ16qX2Hyod1W3FAviFekDAt1Mycy7Jbob2jYv4m+mvDelWamZNbT/jC1nSOPWv6b+BVHp/JS7Xwm2ApfEJt7mbgZ+in4YbwiTIG3FtoWCCXc2GLeV1P5Wkn7LliO83F7xTiLU3m9EGQez+EuvEGYYNOFg1P5j4b96/TVlPeuNDMps6b7a6poxP+ssoaErcViNwo/pC1GcGWhbgW+hNtydetwsjAD7hM+01phex8ibuWjNcMs4V/BdSXtVwpT8NSasd6fyodK2h/GErxPBA8GgROFPT4HH8a+4iAtb0Bbp682vHelmUmZNd1fU0Fj/usO02GYi9+b6Pw2weUiavIvcSPNw3HCmb5WROj+nut/nni2L8NXcvWP4ArhlDfBcmHqXCNs9yLOxO5ik26oGWtOKteXtGf1cxuurQlOFIGUDNeJyN3TDWjr9NWG9640MyWztvurKxrzX2fmHZXKH3ZYxFm4GU/iRdyLr4rXZ2thd+ZxEq4SB2e+sKk/JAIPP8N3G8z5dXxTRI8O79O+p7hdzhVO+lSxWSqrIoZtMZLGHRGRyXkiRL5HA9oqfXXhfdDyYnAya7u/uqAV/1WH6QPYRziy1wxgYRmySM6CXN1CERr/DU4QB+hFESr+DB4Th2RexbjHCpv2PmE7ryu0Z8/1Qzij4VqzW3ROSft2hX6DxJMih7ZE5Jl+UtO/Sl9deO9Cw6aVGf33Vxd05b8vzhe3x+hUBypguzTuS7m6c/QPpWb4pepQ9zdS+z9NDq1nmKs+yZn9stzF0vR/2cuchYA/XtLO1KN5xMs0ZnJwJo8qfXXhvQsNm15m/fZXlzla81/mM80WZtJGERUaJLKvGfJ5o61S+bYSmqz+lT5tJws/6R58QiR7++Fl5bzsIeziO/Gg8Sf9llQuEa94Pjq1rUg+bsDdJeMOCjuksiyaWqevLrx3oWHTy6zf/uqCrvxPQvbJxG9rJpwvsshbFuo/iDf36b+TiOaMmRgZOcz4ZyDvKNAcKBSyweTPas5IdH8pma8pRg02aZvHQvW37M7CRypic+NJ27sq6Jvqqx9GtU9A1tFMt8za7q8uc1RhVM3nRHlkjuyPaga9STDwbhMTaIfiFHFLrRDRlvni27vZwqY/J9f/KpFH2h/3C19htYjvf1I4raeY+AnKF/FtcVvfIYIPRfREQGMqOEZ8GnOBME3uF0nmRcKePq0PzafTj/FDsnduLWtE1C7DASKrfrvIsK8VEb39hJ+42sQIZxFN9TVTmG6Ztd1fXeYYCHYRp26l+gRYL/V9V6F+P/xcRNWeEUnUp0W+6AjjEZ08thS+z93ii4fXRDj8d8JkKGJUvS17a836i2OV3bQ7ilDsE8LU/I/wUcpew7q19Qr9d8VFwlRdI3hfjz+nsape3Tb6qlrrIF8mpldmXfZX2zmqkI3T+nOiIYYYYoghhhhiiCH+7/A/50uR0wPVRi4AAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle 7.582401374401514$"
      ],
      "text/plain": [
       "7.582401374401514"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.cond(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAN8AAAASCAYAAADbjwtGAAAABHNCSVQICAgIfAhkiAAABrtJREFUaIHt2musHVUVB/Bf60ULKhqottEoIJhQ2kSsio8EiggqaknB4Acij6io0Qj4wne48REEiUpUNA0mSogmihRD8VEkxQKRIApCxQdIbyzBprZaQCiPluuHtSdn7tyZc/bMOb39cv7Jydw7e6+9XnvWXnvtzRhjjLHXcRamB/x2txjvItyAzdiJ/+AOXIADM+hPL/F9X59+b8c6PJD43I+f4vWZcubyORo/w7/wRHquw9tq+s7De3ArHsFjQvdz8Iya/gcm3mtwX9LjIdyM92L+iGhob69h/UiejdvabErzPN2yF3XJlmte6e8jsaqB4dE4DtfhHZlCPok/4h5sxbPxOrwaD6a/NzfQvgR3C6M/B2fj8pp+F+F8bMc12IbDcBImcAau7CNjLp/P40tp/LXiw1uIV2J9kqGMK4STtuJaPIrjcYT4gE8VzijwQXw3jbse/8QinILnjZCmi72G8SP5Nm5rsyk8H9+sGet/uKTm/Vzo0kWuvvidUPykFjQLGt5/JY11WUP7PPwG/8DXNEeYxWIl3oIXVtremOju7yNfLp/C6dfjuTXt+1T+X1XivbDSb01qO6tCcxxWmr1aLRYf1TTeOSRNV3t19SP5Nu5is6n0a4O50KWLXI1Ylhg9oH75b4tX6E3mOpyLp3EMJjUr+drU9vOGcR4W6UsTcvjMFxPiUbygz1hlXJHG+nBNW2HLP2SOBZ9NNN8akmZYe1UxyI/k+7KLzaaMbpKPUpdsuSYy+nwgPb+v3Z6vCSvT866atiX4Ki7FBhHdm3CvSCOOEtFyW6ntGLFKXdNAm8vnDTgEV+G/Yr+0DI/jNpERVLE4PetWkeLdcpGa7GjgW8ZT6bkro28/mmHsVYd+fqSdL7va7Fl4N14qAuRdiVfbeTpKXUYm175i0u0W+W4XfEJEim/gJhEx/mT2SjKB2/G3xJf+EQbOE9FoK1bjQvxEfBzrzE6v2vL5aHr/bWHA6gb6tzV6/Ci1faiGdxHFp8X+YhAmxB5jGm/J6D+Ipou9CuT6sZChjS+72GzKbH8UqeuKPnrsaV2GkWsGzkyEa9sQVbClIsQvRWGgii+Kj7xcdZs0uAq5SlStyjzuxWkN/dvwuTC935XGfJPYaC/Fr1LbjRWa09L7+3BA6f2EKBwUMp7YR6cCl6S+12X0zaVpa68CuX6kvS+72OwCsQItwn7iI/2eCC6PiVRyb+gyjFwzcEtisnJQxwwswskigjwo0ogCR4kJfnGFZlL/j+/8RPd1vEwouxy/TnTV8dryuVjviKVqtH1FZWzaTMfMxy/0SsurReVroyhv/z21vblBpwLnpH5/MXNCDkPT1l516OdHuvlyVDajF3zWZPTdE7qMQi5HpM6bjabQUuAgcU62Mf0/IZS/R+TKZUxqVvLY1HZ1Tdt+okC0W0yyrnw+o7cy1OHy1H5u5f0EPo47xeR5WKyUr9KrHB/ZMCZReJjGn/X2Q4MwiOZY7ew1CFU/0t2XBe0wNitwWOq7PaNvgVHrMrRcl6bOky0Y5OKONPZCsYmuy5HrfuWzkyKSfKSBx9WpvSi3d+FzSvr/9w08ipLzp3OUFqvlTpF+VI8oCpyXxrxb/z1YW5q29spB2Y9092U/5NisjP3T+I9njl9gT+syS66maucCceD5tKhyjhovSs/dIuI08VguDrJvFlGoXF0solHTEUDx/sn07MJng0g7Xo5nlsYqsCw9pxrGreJ0Ydsf6lUky/iUqKrdiRPMrEg2IZemrb1yUPYj3X3ZD4NsVkWxBeh3xluHPa1LtlzFFZprB/Q7FIebHZEOV5/6zNc70LxlkBD6L+/v0tsjvLjSdqIIHDvlXRvqx+fK1PblyvsTEo8dIkqWsX/NOK8RhY5H1Kd2X0h8bpe/x2tD08Veo/Ijg1O1NjZbql7fg8QWYVqcdZYxF7q0kqtp5Xt/eq4eIMQNaeBDzIz+bxUp2QZxI2C72NyuEEbcIq7mDIOrxI2D40WBYU0ad4m4AjdPpINtcv86fEwcUH9OnIfdJnQ+WUTJs80+r7teTOSNYuIsFXdAnxCpbDX6nalXVbtJFE6qmMIPhqDpYq+58GOBNjY7Ncm6HptS/0PFOewCUbypXuOaC126yDUDS+QXWqZS34Mr75fhOyIV2iZSt4fE3mlSfmSf1D9a7iP2O7eKDfoucYa1Vl5lLJfPAaJCuEmkZdvFTZGms7pPihsZO8Tk2STKzQcP4N/vd+MIaNraa1R+LMvbZOM2NluBH+Ovqf9T+Lf4gM8w887yXOrSRa4xxhhjjDHGGGOMMcbYk/g/IslsdYqHmvgAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle 7.3484692283495345$"
      ],
      "text/plain": [
       "7.3484692283495345"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# direct approach to solve the linear problem is to compute the inverse of the matrix A and multiply it with the vector b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[-1.,  2.]])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.linalg.inv(A).dot(b)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# However, this is not the most efficient computational method to find the solution vector x.\n",
    "# A better method is LU factorization of the matrix A \n",
    "# A = LU where L is a lower triangular matrix and U is an upper triangular matrix.\n",
    "# Given L and U, the solution vector x can be efficiently constructed by \n",
    "# first solving Ly = b with forward substitution and \n",
    "# then solving Ux = y with backward substitution.\n",
    "# Owning to the fact that L and U are triangular matrices, these two procedures are computationally efficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# In SymPy we can perform a symbolic LU factorization by using the LUdecomposition method of the sympy.Matrix class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = sympy.Matrix([[2, 3], [5, 4]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = sympy.Matrix([4, 3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAA0AAAASCAYAAACAa1QyAAAABHNCSVQICAgIfAhkiAAAAMJJREFUKJHF0j9LQlEYB+DH6BNYhJ/HxcmpD9CgQxCImxAIbi5u7eLSJji5NzQGDRItDtHQ3J0aSnS4LxJytBMN/eBw4XCe+77nD3/MCdqYYYkPFLhHC0cpdIk13nCLIcZ4j/kpKruojmbijzW8Bjz/TevXgW7s6zORz/h+5VY5xiIqNXLRKMA8F3QCPKOaA64CPClP8Md0AyxwlgN6AR5xmgP6AR4c2MP3Z3GBCVbKSywS619izTaDqHJo3OW0+8/ZALBUMPfsVQVVAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle 2$"
      ],
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.rank()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGMAAAA2CAYAAADTeCfRAAAABHNCSVQICAgIfAhkiAAABMxJREFUeJztm12IVVUUx39+jYijhFPmRzQgMaA00AeWKOrNStTyA9EmIoRMHxrzzRd98QqigiiTZA8JcR+DCFEf/MAXC6NBBKOH0oKsB8sPpCRE09CHde547nafs9e5d597z1z2Dw7M2WedvdfZ6+y91z73PxAoDCMy2k8E/snDkUB2NgG9rXaiXRmZ0f554Mc8HAlkC0YncCsvRwLZgrEEOJWXI4FsHADGtNqJdkY7MkYhmdc9T+2+Cjxo8VE4Rivt5gLfeWx3NtnT6rZHOzKWAsc9tTkKfyOsrdAG4wngb09tzsPvKGsbzGBsBl4yynqAix7b7CXsVazEg/E2sAF4x7BZDhy13LsVOIfsPa4Dx5BNoQvb4rkguv9KdH2NxeYy9oX4oGHXD/wG3AHOA/MVPtWLqw+0Plt5HfjJKEu68STwQdR4L3AY+AuYlFJ/T9SGyVJgJ7Ca5GA8BUyJHW9EtqWYTR+yHm0EZiLp+L/Asyk+mVSAstLW1QcanxMZDdxEOg2gC9ihdKwT+B8ZSVXWAW/Fzjfi3qskBcNkAPiV2qxsEDhk2P0C7FbUV6WCPhgmtj6IY/N5CHPNuI9kTSuj82Xos6gJUX03o/OxwDXgvZjNGPxkUh3A+8AXPJr2OoCXefwrwSkkNW8GZh/Esflcgy2bOgKsiP5+BXnbNAwAF4Dvo/O7wAmgGxlxXcANZV0uViEZXiVW9iSSNl81bK8iU0QzMPsgjs1nJxOA28B0YL/ynr3An8Bzlmu7gUXAu5EzLjTT1ElksYwzLbrXXLC3Az+n1LUNWVeqxz3gP6NMkwSk9UGSzypOAF/yaISksQ95+2YlXJ8DfAJ8pGzbFYxuZF5eaZR3INPsWqP8IHAmpb5JSAdWj6+RhT9eNs7hs6sPknxW0Y+MDpcTBxxOgCxWg8B6ZduuYJSRN9D2KWcQ+Nwou0S+C7imD8ok++zkGWTtSOMzJL9eRG361mmxPQTMSKmrE3ghOh4AW6K/zZR0JPA7sCehnj5kitmApLYDyDTT7XiWOBX0wdD0gctnFeMd15O+hpYttl2OukoJdVUMu8VReQ/J9CObrbvIpm+Bo22TCvpgaPpA43MgEAgEApmofiMp5M+QgUAgELS2w5qgtc2RoLUtEEFrWyCC1naYErS2OdMqrW2RKIzut1Va2yJRGN1vK7S2RaJQut9WaG2LRKF0v41obTXyzss0IG/MQL1SU5vut2VS00a0tiXkN+C5yG/A94HT1Mo7ZwNTY8ebUflXaU7VgcYXG7av1eOBH4CPU+7TPFcfoorZBbwInEWmerXUNIvW1sQlbQSHvNEjGl+SdL9xmio1NaepM8DT1GpttSrANGkjKOSNHrH5Yup+FwLfeGjLm9TUp9Y2TdoIdcob68T0JU/drzepqS+t7V5k4VuLTA82PkQCeyWlnjLuDVKpDl/y1P2mPZc5A4ywlA1h2/QdR6I8HXlzXFPKPmSYvobMmza6kf9NWO2o61NEVprGHw348i0SqMlIcBol6bluIC+COQom8/hocaLV2mqkjdCgvFGJxhffut8y+UpNAZ3WVivv9CJvdKD1RaP7LZLUFNBpbbXyzmbIG7NITV2631JCXRXDrhlS0yFcWtvhikv3GwgEAoFA4XkIMLfCFiPzP/MAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\frac{\\sqrt{2 \\sqrt{170} + 27}}{\\sqrt{27 - 2 \\sqrt{170}}}$"
      ],
      "text/plain": [
       "  _____________\n",
       "╲╱ 2⋅√170 + 27 \n",
       "───────────────\n",
       "  _____________\n",
       "╲╱ 27 - 2⋅√170 "
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.condition_number()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMYAAAASCAYAAAAANKFkAAAABHNCSVQICAgIfAhkiAAABgNJREFUaIHt2muMnkUVB/AfWEJBoFVBNygRWpGLJEQwKCqlLdKUcBfhEyDGCkTBghAgKLqYkBREAxIvEG4BDB8gAl6432+SaBRQAbn1hRYo0DYWkApilw9n3uyzz77Pdd/dfnn/yZvZnZlzZubMmZlzeRhggAFKcQxGKn7/b8CvU8JnRQHN/rgdy7EWz+M67Nmj70ewCDfg2dR/DR7EN7FhzXkelZnXooI+n8DleBnvpLVdgA8V9P8aLsIDeCPxvqZiHufiLiwTa1mNv+FHYq15HGPi+1Vn7W1pJltmHc31q/YY0zJ/P4qzCxjuhfm4pWSivbBGCCOPt3rUnYvTsAo3YiU+hYNxGI42dhGH41d4BffgRXwMX8Wl2C/1GSmZ3zZCUG9hs4I+s/EwPoqb8BT2wGIsxJfSnLP4AXZNfJdjx5I5dHEy/oo78Bo+iC9gGMemv5dl+k90v+qsvS3NVMmsiX61HaMUfxIKdlADmk761cGQuN1WCGFmMS+N/Xyufj4ONP5lGBKHZEQcqCJsgDvxHH6i+Aa8LbWdmKv/War/dQ+aedg+jTFXvRdjekH9OYn+lxX0WVTtV921t6WZCpl11NevtmOUYpfEYDk+0ICuo/7EP5/GuKmg/Q282WDsMxO/i0r6LMY6zBG3cq+NnpXqlxp/ADcXN89/xO1ehLkmtgG7Jvo7avavs1911t6WZqpk1tH8YNQeo44dflwqL9PMx4CNcaRQ1MXixPbarGfwrnhut8y1zRECvbPBuP9L5XsF7TthCS7E/SV85qfydqEUWbyJh7CpMHMmCwem8vGa/av2q+7a29JMpczq6ldjTKto3yQNvE7Y7U0xhKtzdUvxDdyXqVuN08VT+4TwMVYJW/UgcVsepx6mCX8Ebi1ov1qYW2dW8NohlU8XtD+DBfi0cJz7gVOF/T4Dn8OXxaFYUoO2ar+arL0tzVTKrK5+NUbVwTgCM/FHYx2/OrhCeP//FDfFLJwgHMlbRKTpsUz/C8TTeDm+lal/FlcKh7QOlghz4mZh6+bxQ3xWKNzaCl4zUrmmoL1bP7Pm3OrgVBFE6OJWEYF6vQZt1X41WXtbmqmSWVP9aoQqU+rYVF7cgvfZuBuv4m38A8eLV2ETYadmcRquF4dgtrBBdxdO929wXo0xv4tTRBTkqB7te4hb76fCQZ0oNkhlWeSrKYYS3yERYZslwra71aAt2682a++3vOifzJrqVyOUHYyd8UXhxN08kUFy6EYk5mTq5opw7e/wPXEY3hbhy0PxklD4WSV8vyNs4CeErbk61941CZ7GWTXn2r3dZhS0b5Hr10+8KnI0C0Qe46qK/mX71WbtbWhYvzKjt371FReKUz3cZ75bJL7/zdSdr3d4r4vfKg+/npTa/258uLeLmaoTYt1fNza+KP1f9GJ2w5L7FLTTh7CgeDFGjA9MZFG2X23W3oaG9S+zXvrVeIwiH2O6MEXWiehGP9HNYmfzEhuncqsCmm79uz3aThd+xaPYVyQGe+EdxWvZTdjRD+JfRs2Ge1K5QLyu2SjL5iJRtRaPFPDtF7ZOZVFUsGq/2qy9DQ3rX2a99Ktv6Kb9f1/Rb7bIHm6Uq/8MPtyj/ydFVGLE2AjHEUZT+R/P0ewnhLvW+E8jzkp0fykYry6G9TfBl8Vc1bffjsKnyGNDowm+h0ro6+5XLwxr9klIHZrJlllT/Wo8RtGL0XXiLqlgfleazHbGJlsOxxni9lgqogazxbdQ04UNfH6m//UiT/EVPCls6xUifn6AcNjOMPYzgq/jx+IWfUA43nl0hDM/EXxbfN7wc/H8PykSkvOE/f39HjSHpB+jCr9nZi4rRfSpi4Uim3y/yCyvEpGpvYVftcLYSF0edfdrqjDZMmuqX23GGIedxElapjpZ0kl9t83V741rRXTo3yLh9rrIRxxtNDKRxUbCV3hEZLrfEyHaP4hnOY9h1bbvvRXzz/MqugG3EeHBV4Q594Kw6Yteqaq5dXL9d8EvhDm4Uqx9Df6ceJW9hk32q2yu/XwxmFyZtdGvpmMMMMAAAwwwwAADDNAO7wMuFWmq1gECmQAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$\\displaystyle 7.58240137440151$"
      ],
      "text/plain": [
       "7.58240137440151"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.N(_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACwAAAAbCAYAAAAH+20UAAAABHNCSVQICAgIfAhkiAAAAmBJREFUWIXt1tuLzVEUB/CP4xKhMG5NSh6Uy5RLHph4EjLyIJd/A5FSRBkmlygNRfLkDQ8iI5JbKMo9hdxSSBhEbjPGw96T4+ec8/sdnZPSfOvX3q219netVvv3XZsuVBfdqszfHx+qnKOiWIQ5/7qIctCMXpUkzFWSLIHu8ftWSdJqFlyPS1XkrziaUFNp0mp2uAZvKk2aLHgzTuMZPuMtrmOd8ro1GvcyxM3AYbzA17ieREOxA0kd/oZruItX6IupmILncf8sQyErcBT3S8SswQa8xrFY7GBMwhmsypBH7yL2jejA7iwk2JPiXxz5TgnDJYmeGfMUxYS8BJ3IoRFzE7GDsKkEVw6P8AlDyi2kR8a4+XG9lWebjWmoRUuevQHHS3DVYxQOoRXzUIcvuILLGWv6DSuxHjtwQejuTX92ZKJw1/N/3l1Kq8/yyNcsNKAj8Z0rkCcVLxMkLRhWJPYJpsd9T+xM4W6KnG14gJnoh/E4EX1nyy24E8OwQJCo55hcIGYntsb9LCxM4dwSi2oX/o189BFUqEO4bn+NkYJO3ingm+mXfG0VulUKq2NBD4r490X/0kLOrJPuqaDN4wWtzEfnnRsr6PbHFK7OgfKuiL81rn0KOcsZzbVxbU/Y2wRVWIvbGXjOxzOjFX561sX1SRrRGAwvYM/5NTguFjm7JPpHpCWJOBDjGxP2WfghdH9AGskyfBfeEnuFv3k/HkbyFxhX5Gx/5T0lhwp3uEPo+DYcFDr/XZiEqagTNPSGMN/b8B5XBU0elHJ+YBkFi3zb8Vh4w7zBEeG90oUu/Lf4CX0Kich8YfSLAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$\\displaystyle 3 \\sqrt{6}$"
      ],
      "text/plain": [
       "3⋅√6"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.norm()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "L, U, _ = A.LUdecomposition()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0\\\\\\frac{5}{2} & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡ 1   0⎤\n",
       "⎢      ⎥\n",
       "⎣5/2  1⎦"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}2 & 3\\\\0 & - \\frac{7}{2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡2   3  ⎤\n",
       "⎢       ⎥\n",
       "⎣0  -7/2⎦"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}2 & 3\\\\5 & 4\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡2  3⎤\n",
       "⎢    ⎥\n",
       "⎣5  4⎦"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L*U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A == L*U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-1\\\\2\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡-1⎤\n",
       "⎢  ⎥\n",
       "⎣2 ⎦"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = A.solve(b)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-1\\\\2\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡-1⎤\n",
       "⎢  ⎥\n",
       "⎣2 ⎦"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = A.LUsolve(b)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For numerical problems we can use the la.lu function form SciPy’s linear algebra module.\n",
    "# It returns a permutation matrix P and the L and U matrices, such that A = PLU.\n",
    "# Like with SymPy, we can solve the linear system Ax = b \n",
    "# without explicitly calculating the L and U matrices\n",
    "# by using the la.solve function, which takes the A matrix and the b vector as arguments.\n",
    "# This is in general the preferred method for solving numerical linear equation systems using SciPy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.matrix([[2, 3], [5, 4]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = np.array([4, 3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "P, L, U = la.lu(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1. , 0. ],\n",
       "       [0.4, 1. ]])"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[5. , 4. ],\n",
       "       [0. , 1.4]])"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 1.],\n",
       "       [1., 0.]])"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2., 3.],\n",
       "       [5., 4.]])"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P.dot(L.dot(U))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ True,  True],\n",
       "        [ True,  True]])"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A == P.dot(L.dot(U))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.,  2.])"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = la.solve(A, b)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = sympy.symbols(\"p\", positive=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = sympy.Matrix([[1, sympy.sqrt(p)], [1, 1/sympy.sqrt(p)]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = sympy.Matrix([1, 2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{- 2 \\sqrt{p} + \\frac{1}{\\sqrt{p}}}{- \\sqrt{p} + \\frac{1}{\\sqrt{p}}}\\\\\\frac{1}{- \\sqrt{p} + \\frac{1}{\\sqrt{p}}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡        1 ⎤\n",
       "⎢-2⋅√p + ──⎥\n",
       "⎢        √p⎥\n",
       "⎢──────────⎥\n",
       "⎢       1  ⎥\n",
       "⎢ -√p + ── ⎥\n",
       "⎢       √p ⎥\n",
       "⎢          ⎥\n",
       "⎢    1     ⎥\n",
       "⎢ ──────── ⎥\n",
       "⎢       1  ⎥\n",
       "⎢ -√p + ── ⎥\n",
       "⎣       √p ⎦"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = A.solve(b)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "# condition number in the numerical case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Symbolic problem specification\n",
    "p = sympy.symbols(\"p\", positive=True)\n",
    "A = sympy.Matrix([[1, sympy.sqrt(p)], [1, 1/sympy.sqrt(p)]])\n",
    "b = sympy.Matrix([1, 2])\n",
    "\n",
    "# Solve symbolically\n",
    "x_sym_sol = A.solve(b)\n",
    "Acond = A.condition_number().simplify()\n",
    "\n",
    "p_vec = np.linspace(0.9, 1.1, 200)\n",
    "for n in range(2):\n",
    "    x_sym = np.array([x_sym_sol[n].subs(p, pp).evalf() for pp in p_vec])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, '$p$')"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEYCAYAAACwQCa4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de3hc1Xnv8e87ulq2jC2MhWMbcMAhkFBuCjiXJgJSbm1q2kJLmwafhNbpCWnTS9qGtqe0JPQkOecpSc5J0voEJ4YkEEJDoC0NdcHTliZgIBDuxMZgW9hgY/kmy7rNvOePvUbakkfSbI1Gmm39Ps8zz+xZe+0972yN3lmz9pq9zN0REZGZITPdAYiIyNRR0hcRmUGU9EVEZhAlfRGRGURJX0RkBlHSFxGZQZT05ahjZv/NzB6KPe4yszePUf9ZM2ufkuAqxMyyZvZb0x2HVD8lfZkyZvYbZvZYSMI7zexfzOw9lX5ed5/j7ltCDN8ws8+MWP82d89WOg6RaqCkL1PCzP4Q+ALwN0ArcALwFWDldMYlw1lEeeEopj+uVJyZHQPcCFzn7t9z90Pu3u/u/+jufxzqNJjZF8xsR7h9wcwawrp2M+swsz8ys13hW8KHY/s/1szuNbMDZrYROHnE87uZnWJmq4EPAn8Svm38Y1j/ipm9v9w4irzurJl92sz+y8wOmtm/mtmC+L5G1I/H8Vdm9l0z+2bY9mkze4uZXR+ee7uZXTziKU82s41mtt/M7jGzlti+V5jZD81sn5n9JN6dFeK8ycz+C+gGRu0Kk/RT0pep8E6gEbh7jDp/DqwAzgLOBM4D/iK2/njgGGAxcC3wZTObH9Z9GegBFgEfCbcjuPsa4FvA50OXzwcmOY5ifgP4MLAQqAc+OUbdkT4A3AbMB54A7if6n11M9CH69yPqX0P02t8EDABfAjCzxcA/A58BWkIM/2Bmx8W2/RCwGmgGtiaIUVJGSV+mwrHAG+4+MEadDwI3uvsud98N/DVRIiroD+v73f0+oAs41cxqgF8B/jJ8g3gGWFdGrBOKY4z9fd3df+ruh4E7iT5MSvWf7n5/OG7fBY4DPuvu/cAdwElmNi9W/zZ3f8bdDwH/A/jVcHx+E7jP3e9z97y7rwceAy6PbfsNd3/W3QfC/uUopaQvU2EPsMDMaseo8yaGtzC3hrLBfYz40OgG5hAlwlpg+4htJ2qicYzmtQR1R3o9tnyY6IMzF3vMiP2NPAZ1wALgROCq0LWzz8z2Ae8h+mZUbFs5iinpy1T4EVH3yxVj1NlBlJwKTghl49lN1JWxdMS2oxnvsrITjSOpQ0BT4UFokR83evWSjDwG/cAbRAn9NnefF7vNdvfPxurrcrszhJK+VJy77wf+kqj/+wozazKzOjO7zMw+H6rdDvyFmR0XTnb+JfDNEvadA74H/FXY7+nAqjE2eZ2xT1ROKI4J+CnQaGY/b2Z1ROcNGsrc52+a2elm1kTU539XOD7fBD5gZpeYWY2ZNYYTyUvKfD5JISV9mRLu/rfAHxIlt91Erc+PA98PVT5D1M/8FPA08ONQVoqPE3VzvAZ8A/j6GHVvAU4P3RzfL7K+nDhKFj4IPwZ8DXiVqOXfMeZG47uN6PW/RnTi/PfCc20nGhr7Zwwd+z9G//8zkmkSFRGRmUOf9CIiM4iSvojIDKKkLyIygyjpi4jMIGP9WGbaLViwwE866aQJb3/o0CFmz549eQFNEsWVjOJKRnElczTG9fjjj7/h7sV/9+HuVXs799xzvRwbNmwoa/tKUVzJKK5kFFcyR2NcwGM+Sl5V946IyAyipC8iMoMo6YuIzCBK+iIiM0hJSd/M/iBMHv2Mmd0eLti0zMweMbNNZvYdM6sPdRvC481h/Umx/Vwfyl80s0sq85JERGQ04yb9MOvO7wFt7v52oAa4GvgccLO7Lwf2Es0iRLjf6+6nADeHeoSrH14NvA24FPhKuJysiIhMkVK7d2qBWWESjCZgJ3AhcFdYv46ha6WvZGjmoruAi8zMQvkd7t7r7i8Dm4mmohMRkSky7o+z3P1VM/vfwDai2Xr+FXgc2OdDMwh1EM3bSbjfHrYdMLP9RNPlLQYeju06vs2gMHn1aoDW1lay2WzyVxV0dXWVtX2lKK5kqi2urQdyPL8nT29fLwP5DdRmbLpDGqbajleB4kqmUnGNm/TDpM8rgWXAPqK5Oi8rUrVwjeZi/wE+Rvnwgmjy6jUAbW1t3t7ePl6Io8pms5SzfaUormSqLa5f+eoPeXzrXsD44CVn0nZSy3SHNEy1Ha8CxZVMpeIqpXvn/cDL7r7bowmTvwe8C5gXm/N0CUNTynUQpm0L648BOuPlRbYRSY2e/hyz66PTUX25/DRHI5JMKUl/G7AiTEVnwEXAc8AG4MpQZxVwT1i+l6Hp6q4EHgw/C74XuDqM7lkGLAc2Ts7LEJk6ubxTXxv96+SV8yVlSunTf8TM7iKaNm4AeIKo++WfgTvM7DOh7JawyS3AbWa2maiFf3XYz7NmdifRB8YAcJ1H83eKpErenbqaKOnnNPOcpExJV9l09xuAG0YUb6HI6Bt37wGuGmU/NwE3JYxRpKrk8kNJP59X0pd00S9yRRLKO4PdOzklfUkZJX2RhKKWfjQYTd07kjZK+iIJxbt3XElfUkZJXyShYSdyNXpHUkZJXyShXN6p1+gdSSklfZGE8u7U1UZ9+hq9I2mjpC+SULxPX6N3JG2U9EUSGpb01b0jKaOkL5JQ3hns01f3jqSNkr5IQhqnL2mmpC+SUM51GQZJLyV9kYTyeadOl2GQlFLSF0ko5/Fx+tMcjEhCSvoiCbg77gz26at7R9JGSV8kgUJ3joZsSlop6YskUEjytfpxlqTUuEnfzE41sydjtwNm9vtm1mJm681sU7ifH+qbmX3JzDab2VNmdk5sX6tC/U1mtmr0ZxWpToXpEevVvSMpNW7Sd/cX3f0sdz8LOBfoBu4GPgU84O7LgQfCY4DLiOa/XQ6sBr4KYGYtRLNvnU8049YNhQ8KkbQotPTVvSNplbR75yLgJXffCqwE1oXydcAVYXklcKtHHgbmmdki4BJgvbt3uvteYD1wadmvQGQKFbpzajJq6Us6lTRHbszVwO1hudXddwK4+04zWxjKFwPbY9t0hLLRyocxs9VE3xBobW0lm80mDHFIV1dXWdtXiuJKppri6uqLkvyWl14ig7Plla1kszunOarhqul4xSmuZCoVV8lJ38zqgV8Erh+vapEyH6N8eIH7GmANQFtbm7e3t5ca4hGy2SzlbF8piiuZaorrja5eePDfOPUty8m8+CxLlp5Ae/tbpzusYarpeMUprmQqFVeS7p3LgB+7++vh8euh24ZwvyuUdwBLY9stAXaMUS6SGvnQh5/JGBkbeiySFkmS/q8z1LUDcC9QGIGzCrgnVn5NGMWzAtgfuoHuBy42s/nhBO7FoUwkNQqjd2osSvoasilpU1L3jpk1AT8HfDRW/FngTjO7FtgGXBXK7wMuBzYTjfT5MIC7d5rZp4FHQ70b3b2z7FcgMoUKo3VqMmBK+pJCJSV9d+8Gjh1RtodoNM/Iug5cN8p+1gJrk4cpUh0Ko3Uypu4dSSf9IlckgfiQTXXvSBop6YskMNS9Y2TM1NKX1FHSF0lgWPcOaulL+ijpiyQwvKUfzZcrkiZK+iIJ5EaeyFXWl5RR0hdJYHCcfuFErvr0JWWU9EUS0Dh9STslfZEEjujeUUtfUkZJXySBfPxELmrpS/oo6YskMPjjLIvG6efy0xyQSEJK+iIJDI7T11U2JaWU9EUSGDlOX907kjZK+iIJ6ESupJ2SvkgCebX0JeWU9EUSyMUmUTGU9CV9lPRFEhjs3smg7h1JpZKSvpnNM7O7zOwFM3vezN5pZi1mtt7MNoX7+aGumdmXzGyzmT1lZufE9rMq1N9kZqtGf0aR6qTuHUm7Ulv6XwR+4O5vBc4Engc+BTzg7suBB8JjiCZQXx5uq4GvAphZC3ADcD5wHnBD4YNCJC3i4/TNjJxyvqTMuEnfzOYC7wVuAXD3PnffB6wE1oVq64ArwvJK4FaPPAzMM7NFwCXAenfvdPe9wHrg0kl9NSIVVmjpD47TV0tfUqaUOXLfDOwGvm5mZwKPA58AWt19J4C77zSzhaH+YmB7bPuOUDZa+TBmtproGwKtra1ks9kkr2eYrq6usravFMWVTDXF9cyr/QA8unEj+dwABw4crJrYCqrpeMUprmQqFVcpSb8WOAf4XXd/xMy+yFBXTjFWpMzHKB9e4L4GWAPQ1tbm7e3tJYRYXDabpZztK0VxJVNNce1+bDs8/RTvWrGCO1/8D5rqm2hvf+90hzVMNR2vOMWVTKXiKqVPvwPocPdHwuO7iD4EXg/dNoT7XbH6S2PbLwF2jFEukhpD3TvoRK6k0rhJ391fA7ab2amh6CLgOeBeoDACZxVwT1i+F7gmjOJZAewP3UD3Axeb2fxwAvfiUCaSGoUcr0lUJK1K6d4B+F3gW2ZWD2wBPkz0gXGnmV0LbAOuCnXvAy4HNgPdoS7u3mlmnwYeDfVudPfOSXkVIlNk+FU2dSJX0qekpO/uTwJtRVZdVKSuA9eNsp+1wNokAYpUk/joHVNLX1JIv8gVSWBYSx8bnDNXJC2U9EUSyI24nr5O5EraKOmLJHDEZRjUvSMpo6QvkkD8Kps6kStppKQvksDIcfq6yqakjZK+SAIjh2yqT1/SRklfJIHBpJ8xMgz9WEskLZT0RRLIu2MGVri0srK+pIySvkgCubxTY9G1AzV6R9JISV8kgZw7mcxQ0tfoHUkbJX2RBPJq6UvKKemLJJDLRydxIZogwh1ciV9SRElfJIG8OyHnD97rZK6kiZK+SAK5vA+29AeTvlr6kiJK+iIJ5PzIpK8rbUqaKOmLJJDPO5nBE7nRvVr6kiYlJX0ze8XMnjazJ83ssVDWYmbrzWxTuJ8fys3MvmRmm83sKTM7J7afVaH+JjNbNdrziVSrot076tOXFEnS0r/A3c9y98IMWp8CHnD35cAD4THAZcDycFsNfBWiDwngBuB84DzghsIHhUha5DzW0g9lGqsvaVJO985KYF1YXgdcESu/1SMPA/PMbBFwCbDe3TvdfS+wHri0jOcXmXL5WEvfdCJXUqjUidEd+Fczc+Dv3X0N0OruOwHcfaeZLQx1FwPbY9t2hLLRyocxs9VE3xBobW0lm82W/mpG6OrqKmv7SlFcyVRTXDtf66G3J082m6W/rxcwHnrov5jXWD2nx6rpeMUprmQqFVepSf/d7r4jJPb1ZvbCGHWtSJmPUT68IPpAWQPQ1tbm7e3tJYZ4pGw2SznbV4riSqaa4vruqz9md+4A7e3tZLevB/o4/53vZNExs6Y7tEHVdLziFFcylYqrpOaJu+8I97uAu4n65F8P3TaE+12hegewNLb5EmDHGOUiqTHygmuFMpG0GDfpm9lsM2suLAMXA88A9wKFETirgHvC8r3ANWEUzwpgf+gGuh+42MzmhxO4F4cykdTQOH1Ju1K6d1qBuy1q3dQC33b3H5jZo8CdZnYtsA24KtS/D7gc2Ax0Ax8GcPdOM/s08Giod6O7d07aKxGZAu4apy/pNm7Sd/ctwJlFyvcAFxUpd+C6Ufa1FlibPEyR6jBsnH6sTCQtqmfIgUgK5Jxh19MHTY4u6aKkL5JAdD39aNl0IldSSElfJAFdhkHSTklfJIFhl2FQ946kkJK+SAJ5tfQl5ZT0RRIYNk4/lCnnS5oo6YskUOx6+urekTRR0hdJoNgvctW9I2mipC+SQC7PkSdylfQlRZT0RRKITuRGy7qevqSRkr5IAsVO5Kp7R9JESV8kgeEnckOZWvqSIkr6IgkUP5E7jQGJJKSkL5JAfBIVXXtH0khJXySBfN5jV9nUOH1JHyV9kQRyrukSJd1KTvpmVmNmT5jZP4XHy8zsETPbZGbfMbP6UN4QHm8O60+K7eP6UP6imV0y2S9GpNJy+dj19EOZWvqSJkla+p8Ano89/hxws7svB/YC14bya4G97n4KcHOoh5mdDlwNvA24FPiKmdWUF77I1Mr70Dh9tfQljUpK+ma2BPh54GvhsQEXAneFKuuAK8LyyvCYsP6iUH8lcIe797r7y0Rz6J43GS9CZKrET+Qq6UsalTIxOsAXgD8BmsPjY4F97j4QHncAi8PyYmA7gLsPmNn+UH8x8HBsn/FtBpnZamA1QGtrK9lsttTXcoSurq6ytq8UxZVMNcXV19/Pjh2vks2+weHubsB47oUXyHa9NN2hDaqm4xWnuJKpVFzjJn0z+wVgl7s/bmbtheIiVX2cdWNtM1TgvgZYA9DW1ubt7e0jq5Qsm81SzvaVoriSqaa47MEfcOLSpbS3n87dP3gQOMzy5afSfv4J0x3aoGo6XnGKK5lKxVVKS//dwC+a2eVAIzCXqOU/z8xqQ2t/CbAj1O8AlgIdZlYLHAN0xsoL4tuIpEIuNmRT196RNBq3T9/dr3f3Je5+EtGJ2Afd/YPABuDKUG0VcE9Yvjc8Jqx/0N09lF8dRvcsA5YDGyftlYhMgbwXuZ6++vQlRUrt0y/mT4E7zOwzwBPALaH8FuA2M9tM1MK/GsDdnzWzO4HngAHgOnfPlfH8IlMuF7vKpi64JmmUKOm7exbIhuUtFBl94+49wFWjbH8TcFPSIEWqRd45YvSOxulLmugXuSIlKnTjZDRzlqSYkr5IiQonbI8Yp6+WvqSIkr5IiXKjtPR1IlfSRElfpESFvntdT1/STElfpESFlv7g9fQL5erekRRR0hcpUT606Id+nGVkTN07ki5K+iIlGjqRO1RWkzEN2ZRUUdIXKdFg905mKOubmbp3JFWU9EVKVGjRZ2JJv8ZM3TuSKkr6IiUaeSIXola/Ru9Imijpi5Ro5Dh9iIZtqk9f0kRJX6REeR+tpa+kL+mhpC9SomIncmsyOpEr6aKkL1KiYidyMzqRKymjpC9SosIJW3XvSJop6YuUaKh7Z6gso3H6kjLjJn0zazSzjWb2EzN71sz+OpQvM7NHzGyTmX3HzOpDeUN4vDmsPym2r+tD+YtmdkmlXpRIJQx274xo6at7R9KklJZ+L3Chu58JnAVcamYrgM8BN7v7cmAvcG2ofy2w191PAW4O9TCz04mmTnwbcCnwFTOrmcwXI1JJo5/Ina6IRJIrZWJ0d/eu8LAu3By4ELgrlK8DrgjLK8NjwvqLzMxC+R3u3uvuLwObKTLdoki1yhU9kasLrkm6lDRHbmiRPw6cAnwZeAnY5+4DoUoHsDgsLwa2A7j7gJntB44N5Q/HdhvfJv5cq4HVAK2trWSz2WSvKKarq6us7StFcSVTLXFt2psD4Nmnn8J21tLV1UXP4Qyv7eqpivgKquV4jaS4kqlUXCUlfXfPAWeZ2TzgbuC0YtXCvY2ybrTykc+1BlgD0NbW5u3t7aWEWFQ2m6Wc7StFcSVTLXHN2rIHHnmYs886i3efsoBsNkvznAwtLU20t7dNd3iDquV4jaS4kqlUXIlG77j7PiALrADmmVnhQ2MJsCMsdwBLAcL6Y4DOeHmRbUSqXk4ncuUoUMroneNCCx8zmwW8H3ge2ABcGaqtAu4Jy/eGx4T1D7q7h/Krw+ieZcByYONkvRCRShucRGXE9fQ1ZFPSpJTunUXAutCvnwHudPd/MrPngDvM7DPAE8Atof4twG1mtpmohX81gLs/a2Z3As8BA8B1odtIJBVyfuTonYzpx1mSLuMmfXd/Cji7SPkWioy+cfce4KpR9nUTcFPyMEWmX77IVTY1c5akjX6RK1KiolfZVEtfUkZJX6RExX6clckM9fWLpIGSvkiJBidRibX0azMZBpT1JUWU9EVKdKgvGncwu2Ho6iFN9TV092k8gqSHkr5IiQ729APQ3Fg3WDansZaDPQOjbSJSdZT0RUrUFZJ7c+PQoLe5jXUcCB8GImmgpC9SooO9AzTWZaiLXVC/ubGWrt4B/SpXUkNJX6REB3v6mdNQN6ysubEWd+juV7++pIOSvkiJDvQMMLdx+O8ZC/37B9XFIymhpC9SooM9A8P68wHmNNQOrhNJAyV9kRJ19fQPG7kDQyd11dKXtFDSFynRwZ6BwZZ9QeFD4IBa+pISSvoiJSrWvVPo4+9S0peUUNIXKdHBIt07cxrVpy/poqQvUoJc3jnUlzuipa/RO5I2SvoiJejqPfLXuACz62vImFr6kh6lTJe41Mw2mNnzZvasmX0ilLeY2Xoz2xTu54dyM7MvmdlmM3vKzM6J7WtVqL/JzFaN9pwi1WboujvDk76ZMaehVi19SY1SWvoDwB+5+2lEE6JfZ2anA58CHnD35cAD4THAZUTz3y4HVgNfhehDArgBOJ9oxq0bCh8UItXu4OB1d+qOWNfcWMfBXrX0JR3GTfruvtPdfxyWDxJNir4YWAmsC9XWAVeE5ZXArR55GJhnZouAS4D17t7p7nuB9cClk/pqRCrkYJGLrRU060qbkiKlTIw+yMxOIpov9xGg1d13QvTBYGYLQ7XFwPbYZh2hbLTykc+xmugbAq2trWSz2SQhDtPV1VXW9pWiuJKphrie3BUl9U3PPkXu1eh6+oW48r2H2f7aoWmPsaAajlcxiiuZSsVVctI3sznAPwC/7+4HLDZ70MiqRcp8jPLhBe5rgDUAbW1t3t7eXmqIR8hms5SzfaUormSqIa79T74KP36S977rfE5ZOGdYXLe+8iivH+ihvf1npzXGgmo4XsUormQqFVdJo3fMrI4o4X/L3b8Xil8P3TaE+12hvANYGtt8CbBjjHKRqlf4xe3IC66BunckXUoZvWPALcDz7v63sVX3AoUROKuAe2Ll14RRPCuA/aEb6H7gYjObH07gXhzKRKpesVmzCuY01A4O6RSpdqV077wb+BDwtJk9Gcr+DPgscKeZXQtsA64K6+4DLgc2A93AhwHcvdPMPg08Gurd6O6dk/IqRCrsYM8AtRmjse7IdlJzYx0He/pxd8bo9hSpCuMmfXd/iOL98QAXFanvwHWj7GstsDZJgCLVoCtcd6dYUm9urKU/5/QO5GmsqymytUj10C9yRUpwsKd/8Do7IxX6+TVXrqSBkr5ICQ72DNDccGR/PsSvv6N+fal+SvoiJSh2WeWCwjX2dXllSQMlfZESHChyWeWCZl1eWVJESV+kBHsO9TG/qXjSP3ZOAwC7DvZMZUgiE6KkLzKOw305dh/s5cRjm4quXzJ/FmawrbN7iiMTSU5JX2QchWS+tKV40m+sq+H4uY1K+pIKSvoi4ygk8xOPnT1qnaUtTWxX0pcUUNIXGUch6Z8wSksf4MSWJrbuUdKX6qekLzKObXsO0dxQO+qJXIg+EHYd7OVwX24KIxNJTklfZBzbOrtZ2tI05nV1TggnebfvVWtfqpuSvsg4tnZ2jzpyp6DQ9bNNXTxS5ZT0RcaQzzsdnYfH7M+HoZO8W3UyV6qckr7IGF470ENfLj/YfTOa+U11NDfUagSPVD0lfZExDA7XbBl9uCaAmbG0pYmtew5NRVgiE6akLzKGl3Z3AYzbpw+wbMFsNof6ItWqlOkS15rZLjN7JlbWYmbrzWxTuJ8fys3MvmRmm83sKTM7J7bNqlB/k5mtKvZcItXm8Vf2smBOPUvmzxq37tknzGN752Fe269r8Ej1KqWl/w3g0hFlnwIecPflwAPhMcBlwPJwWw18FaIPCeAG4HzgPOCGwgeFSDV75OVOzlvWUtI0iOcvOzZss6fSYYlM2LhJ393/Axg5l+1KYF1YXgdcESu/1SMPA/PMbBFwCbDe3TvdfS+wniM/SESqSsfebl7dd5jzTmopqf5pi5qZ01DLxpc19bNUr1ImRi+m1d13Arj7TjNbGMoXA9tj9TpC2WjlRzCz1UTfEmhtbSWbzU4wROjq6ipr+0pRXMlMV1w/3BFdHz+zZwvZ7NYj1heLa1mzs+HZ7WTnT19rX3/HZGZaXBNN+qMp9h3Yxyg/stB9DbAGoK2tzdvb2yccTDabpZztK0VxJTNdcd3/vaeY27iTD/7ChdRkjnwLF4vrOTbz+R+8yBlt7xy8zv5U098xmZkW10RH77weum0I97tCeQewNFZvCbBjjHKRquTu/OilPZy3rKVowh/N+cuirqAfbVG/vlSniSb9e4HCCJxVwD2x8mvCKJ4VwP7QDXQ/cLGZzQ8ncC8OZSJV6Scd+3llTzcXvrU10XZnLpnHwuYGvv/EqxWKTKQ8pQzZvB34EXCqmXWY2bXAZ4GfM7NNwM+FxwD3AVuAzcD/Az4G4O6dwKeBR8PtxlAmUpXufGw7jXUZPnDmokTb1dZk+OVzlrDhxd2aPlGq0rh9+u7+66OsuqhIXQeuG2U/a4G1iaITmQaH+3L845M7uOzti0adDH0sV7Ut4e/+/SXu/vGrfPR9J1cgQpGJ0y9yRUb4p6d2cLB3gKvalkxo+5OPm8O5J87n9o3bGMjlJzk6kfIo6YvE9A7k+OIDmzht0VxWhB9bTcRv/+ybeWVPN3c+1jGJ0YmUT0lfJObWH26lY+9h/vzy08gkGLUz0iVva6XtxPn87fqfcqh3YBIjFCmPkr5IsL2zmy89uIn3veU43rN8QVn7MjP+7OdP442uXj7/gxcmKUKR8inpiwB9A3l+9/YnwOHGlW+blH2ec8J8rn3PMtb9aCv3Pb1zUvYpUi4lfZnx8nnn+u89zZPb9/G5K39mcBasyfCnl76Vs5bO44+/+xMe36pRyjL9lPRlRhvI5bn+e0/zDz/u4Pffv5zLz0g2Ln889bUZ1nzoXBbObWTV2kd1MTaZdkr6MmN1Hupj1dc38p3HtvN7F57CJy5aXpHnWTi3kW//9vkc19zAB7/2MLc9vJXoJy0iU09JX2Ycd+efntrBxTf/O4++spf/deXP8IcXn1rSNfMnatExs/j+x97Nu09ZwP/4/jNcs3Yj2/ZoPl2ZepN9lU2RqpXPOw9tfoMvPrCJx7fu5YzFx3DrR36G0980d0qe/5imOtauegfffGQrn/uXF7job7Nc/Y4T+Mh7lrFsweSdRxAZi5K+HPX2dffxvR+/yjcf3sqWNw6xsLmBv/mlM/jVtiXU1kztl91MxrjmnXVSVOkAAAySSURBVCdxyduO5wv/tonbN27jtoe38rPLF/CbK07kfW85jsa6mimNSWYWJX056rg7W/d0s+HFXax/7nUeebmTXN45+4R53PxrZ3L5GYtoqJ3exNo6t5H/+ctn8AfvX87tG7fz7Y1b+ehtjzO7vob3nXocP3d6KyvefCyLjhl/bl6RJJT0JfU6D/Xxws4DPLfzAE9s38fGlzvZfbAXgFMWzuGj730zl5+xiLcvPmaaIz3SwrmNfOL9y/nYBSfz0OY3WP/c6/zbc69z39OvAbBk/izecVILb3vTXE5bFN1aZtdPc9SSZkr6kgqH+3J07O1m655utnVGty1vHOKFnQfYFRI8wJuOaeRdJx/LO05q4V0nH8ubj5szjVGXrq4mwwWnLuSCUxfymZVv57mdB9j4cicbX+7koc1vcHfs+vwLmxtY3jqHpfObWNrSxJL5s1ja0sTS+U0smKMPBBmbkr5Mi/5cngOH+9kfu+3t7mPXgV52HQy3Az3s7upl595DHP7BD4ZtP6ehlhOPbeI9yxdw2vFRC/jU45s5rnl6piicTJmM8fbFx/D2xcfwkfcsA+CNrl5e2HmQ53ce4PmdB3jpjUOsf+519hzqG7ZtfW2G5lpnyTMPsWBOA8c1N7BgTgML5tSzoLmB+U31zG2sY+6sWpob62hurKVuis9ryPRS0pdRuTu9A3n6cnn6BmK3XJ7uvhzdfQN09+bo7s9xuG8glIXyvtywdQd6BoYl+e6+3KjPO6uuhoVzG1jY3MBpx8/l5KZezj7tZBbPm8WJx87mhJYm5jfVVXSIZbVZMKeB9yxvOOKaQId6B+jYe5jtnd1s39vNzv09PLN5GzWz6tixv4enXt1P56E+cvnRfxfQVF8z+EEwt7GO2Q21NNXXMKuuhln1NbHlUB4ex5frazM01Gaor4mWB281GepqbEb9rardlCd9M7sU+CJQA3zN3T87zibTyt1xh5w7+cJyPlrOezQMcHDZxyjPDy1vPZDj6Y795N3JuePu9OecXN4ZyDu5fH7sx7k8A4PrnP5cPlZ3xOOc058Pj3NDSby3Pzc8mefyHOg6jP3n+mFlE1GTMZrqamhqqKGpvpZZdTXMaaxlaUsTb59VxzEjbnNn1XLMrDrmN9WzcG4js+trhiWJaILoUybrT3pUmd1Qy6nHN3Pq8c2DZdns67S3nz/4OJ939h3u542uXvZ193PgcD8HeqL7gz0DYTm633+4n32H+9m5/zDdfTkOhw/yw/2jf0iXor42Qw15mv5z/bAPhMJyXU2G2oxRW7jPGLU1Rk0mQ13GqImvq4nW12SiD5SajFFXk4nqhFtNzdB2hVvGoltNJrogXo1F5U/vHqBm0+7B9RmL3sOZsE2NGRbKov0Q21ehHqFerE4m2jZjRibDsH1N54fglCZ9M6sBvkw0xWIH8KiZ3evuz03m87zw2gE+/u0n6OrqpvHRDeRDovaQhAuJtpCEo3VFlkOSr4gfPjSpuzODukzsjR/+YWozhX+K6L6+NrTKajLMaailvmnoH69zdx8nLj1+WGutodCCG/FPGrX+opbf7IbQCgyJvr4mo5ZdFclkjJbZ9WWdAHZ3evrzg9/ievpzg9/sevpzRb4RDm9U9ObybHllGwuPP36wQRFvXBSWu/tyDOTzDOSGGjHxxwOhwVNoxAzk84zxJaZ0j2+chJ0kE//wwIY/NuCsBdDePvnPO9Ut/fOAze6+BcDM7gBWApOa9GfV1XBqazO7M4dZdPy86CCGT+L4p27hU33Yp7MNfXIP1Rn65I63BKzIcvzTvNASiC9nDJ599lnOPOMMMqHFkTE7ojVTM9iiiSfuoaReU2PDknw5134viFrUZ0zCX0CONmYWdeXU1zDRqWWibyCT//7Kx7/l5vNHfLvN5Ye+UefyQ427wjfxxx57nLPOPjtWL/5t3snnwzf9/PBGYy4/1EjMDdb1wUZmYfvCc+Zj+3d3HIb1BuDDH9cc2DHpxwrApvIaIGZ2JXCpu/9WePwh4Hx3/3iszmpgNUBra+u5d9xxx4Sfr6urizlzqm/0huJKRnElo7iSORrjuuCCCx5397aiKz18ak3FDbiKqB+/8PhDwP8Zrf65557r5diwYUNZ21eK4kpGcSWjuJI5GuMCHvNR8upUj9XqAJbGHi8BKvMdRkREjjDVSf9RYLmZLTOzeuBq4N4pjkFEZMaa0hO57j5gZh8H7icasrnW3Z+dyhhERGayKR+n7+73AfdN9fOKiIgmURERmVGU9EVEZhAlfRGRGWRKf5yVlJntBraWsYsFwBuTFM5kUlzJKK5kFFcyR2NcJ7r7ccVWVHXSL5eZPeaj/SptGimuZBRXMoormZkWl7p3RERmECV9EZEZ5GhP+mumO4BRKK5kFFcyiiuZGRXXUd2nLyIiwx3tLX0REYlR0hcRmUFSk/TN7FIze9HMNpvZp4qsP9HMHjCzp8wsa2ZLYutWmdmmcFsVKz/XzJ4O+/ySTWCOv4nGZWZnmdmPzOzZsO7XYtt8w8xeNrMnw+2sqYorrMvFnvveWPkyM3skHMfvhCulTklcZnZBLKYnzazHzK4I6ybjeK01s11m9swo6y28RzaH2M6Jravk+2tCcU3B+6uc41XJ99dEj9d0v7/eGv5evWb2yRHriv7PTPh4jXah/Wq6EV2R8yXgzUA98BPg9BF1vgusCssXAreF5RZgS7ifH5bnh3UbgXcCBvwLcNkUxvUWYHlYfhOwE5gXHn8DuHI6jld43DXKfu8Erg7Lfwf896mMK1anBegEmibjeIV9vBc4B3hmlPWXh/eIASuARyr9/iozroq9v8qJq5Lvr3Ljmub310LgHcBNwCdL+Z+Z6PFKS0t/cG5dd+8DCnPrxp0OPBCWN8TWXwKsd/dOd98LrAcuNbNFwFx3/5FHR+1W4Iqpisvdf+rum8LyDmAXUPQXdBNQzvEqKrRSLwTuCkXrmMLjNcKVwL+4e3fC5x+Vu/8H0T/6aFYCt3rkYWBeeA9V8v014bgq/P4q53gVNUnvr8mKa8rfX+6+y90fBfpHrCr6P1PO8UpL0l8MbI897ghlcT8BfiUs/xLQbGbHjrHt4rA81j4rGdcgMzuP6FP8pVjxTeHr581m1jDFcTWa2WNm9nDhKy5wLLDP3QfG2Gel4yq4Grh9RFk5x6sUY72PKvX+KieuQRV4f5UbV6XeX+XGVTAd76/RjBbvhI9XWpJ+sb7QkWNNPwm8z8yeAN4HvAoMjLFtKfusZFzRDqJWxm3Ah909H4qvB95K9HWvBfjTKY7rBI9+/v0bwBfM7OQS91npuArH6wyiiXgKyj1epUj6PpqM41WKMZ+nQu+vcuOq1Pur3Lim8/01mkl/f6Ul6Y87t66773D3X3b3s4E/D2X7x9i2IyyPus8Kx4WZzQX+GfiL8FWzsM3O8PWzF/g60Ve8KYsrdAfg7luALHA20YWf5plZ7Wj7rHRcwa8Cd7t7f2ybco9XObFX8v1VTlyVfH+VFVcF319lxRVM1/trNKPFO+HjlZakP+7cuma2wMwKr+d6YG1Yvh+42Mzmm9l84GLgfnffCRw0sxWhf+wa4J6piivUv5uof/G7I7ZZFO6NqJ+u6Bn/CsU1v/D11cwWAO8Gngv90huI+jsBVjGFxyvm1xnx1XsSjlcp7gWuCaM/VgD7w3uoku+vCcdV4fdXOXFV8v014bhi66fr/TWaov8zZR2vkWd2q/VGdNb9p0T9kn8eym4EfjEsXwlsCnW+BjTEtv0IsDncPhwrbyP6A74E/F/CL5SnIi7gN4lO2jwZu50V1j0IPB1i+yYwZwrjeld47p+E+2tj+3wz0YiUzUSjbBqmKq6w7iSi7p7MiH1OxvG6nWiESz9R6+pa4HeA3wnrDfhyiPtpoG2K3l8TimsK3l8TjavS769y/o7T+f46PpQfAPaF5bmj/c+Uc7x0GQYRkRkkLd07IiIyCZT0RURmECV9EZEZRElfRGQGUdIXEZlBlPRFRGYQJX0RkRmkdvwqIhJnZncQ/cjnJKIf1XzM3f95WoMSKZFa+iLJnQlscffzgQ8CN0xzPCIl0y9yRRIws1nANmCpu/eYWQvRRBzLpzk0kZKopS+SzNuBTe7eEx6fQ3QdGZFUUJ++SDJnAieYWSPRVHZ/DfzJ9IYkUjolfZFkzgS+RXQd+LnA37j7f01rRCIJqE9fJAEz+w/gt939xemORWQilPRFEjCzV4lO4ubHrSxShZT0RURmEI3eERGZQZT0RURmECV9EZEZRElfRGQGUdIXEZlBlPRFRGYQJX0RkRnk/wNp6Tyj8RFIbAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "condA_vec= [Acond.subs(p, pp).evalf() for pp in p_vec]\n",
    "\n",
    "# Plot symbolic solution\n",
    "plt.plot(p_vec, condA_vec)\n",
    "plt.title(label=\"Condition number\")\n",
    "plt.grid()\n",
    "plt.xlabel(r'$p$')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Numerical problem specification\n",
    "AA = lambda p: np.array([[1, np.sqrt(p)], [1, 1/np.sqrt(p)]])\n",
    "bb = np.array([1, 2])\n",
    "\n",
    "# Solve numerically\n",
    "x_num_sol = lambda p: np.linalg.solve(AA(p), bb)\n",
    "p_vec = np.linspace(0.9, 1.1, 200)\n",
    "\n",
    "for n in range(2):\n",
    "    x_num = np.array([x_num_sol(pp)[n] for pp in p_vec])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Error in solution\\n(numerical - symbolic)/symbolic')"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEXCAYAAACu1P9TAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZwU5b3v8c9vdoZ9UyGKgLsREMGIS1yiyXFHXKJeXI8REq/mut4TTGJ8aRKTeDQel1wV4zlGE0VNPOoJilsUzYkCOaDiDipLUGRRwHEYZvndP6qqrRl6unu6Z3qmZr7v12teU11VXc/TNTXffvqp6qfM3RERke6vpLMrICIixaHAFxHpIRT4IiI9hAJfRKSHUOCLiPQQCnwRkR5CgS89mpl93czeKXKZz5vZdwp4/udmNro96yQ9gwJfOo2ZfWhmtWGART+3FrMO7v6iu+9WzDLbIt2bg7v3cff3O6tOklxlnV0B6fGOc/dnsq1kZmXu3pBtXlu3IdKTqIUvXZKZnWNmfzWzX5vZeuDqVuaVmNmPzGyZmX1iZr8zs/7hNkaamZvZeWa2HHguTTmHmtnK2OMPzexyM3vNzDaY2Swzq2qljjub2QvhemvNbFZs2QFmNj9cNt/MDmhlG1eb2X2xx1Gdy8zsZ8DXgVvjn37C5TuH0/3D17wm3Ac/MrOS2D58ycz+1cw+NbMPzOyoNv8xpNtQ4EtXth/wPrAN8LNW5p0T/hwGjAb6AC27hQ4B9gD+Kcdyvw0cCYwCxobbT+da4ClgILA9cAuAmQ0C/gzcDAwGbgT+bGaDcywfAHf/IfAicGHYjXNhmtVuAfoTvPZDgLOAc2PL9wPeAYYAvwJ+a2bWlnpI96HAl872n2b2Wezn/NiyVe5+i7s3uHttK/OmAje6+/vu/jkwAzjNzOLdlVe7e01sG9nc7O6r3H098Diwdyvr1QM7AsPdfbO7vxTOPwZ4z93vDet5P/A2cFyO5efEzEqBU4EZ7r7J3T8EbgDOjK22zN1nunsjcA8wDNi2PeshyaHAl852grsPiP3MjC1bkWb9lvOGA8tij5cRnJuKh1q67WTycWz6C4JPDen8X8CAeWb2hpn9cyt1iur1lTbWI5shQAVbv/54OanX4u5fhJOtvR7p5hT40pWlG8q15bxVBK3syAigAVidZTsFc/eP3f18dx8OTAd+E/att6xTVK9/pNlMDVAde7xdy2IyVGEtX37KyFaOiAJfEu9+4BIzG2VmfYCfA7OKcTWOmZ1iZtuHDz8lCOdGYDawq5n9r/Dk66nAnsB/pdnMIuBgMxsRnmye0WL5aoL++a2E3TQPAj8zs75mtiNwKXBfuvVFFPjS2R5vcR3+I218/t3AvcBc4ANgM3BRe1eyFfsCr5jZ58BjwP9x9w/cfR1wLHAZsI6g6+dYd1/bcgPu/jQwC3gN+Dtbvyn8G3ByeJXNzWnqcBHBp4T3gZeAPxDsE5GtmG6AIiLSM6iFLyLSQyjwRUR6CAW+iEgPocAXEekhFPgJZ2bXmdnFnV2PXIRfTjq0wG00G3umKypk+ONwLJ8jwukrzeyuHJ5zo5l9N5/y2lOxX7e0nUbLTDAzG0owdsrOnV2XXLj7Vzu7Dkni7j/PcdXrCb7te7e7b4lmhgO2/au7px24ratqw+uWNlILP9nOAWa3YYyYTtFiXBtpZ+7+EcFYPce3WHQ0wZfARAAFftIdBbwQPYiG+jWzy8Khgj8ys3Njy5t95I6Gz409djO7wMzeM7NNZnatme1kZn8zs41m9qCZVcTWP9bMFoWDnv23mY2NLfvQzP7FzF4DasJvnMY/tpeGH92XhmX93cx2CJf9m5mtCMv8u5l9vT12lgV+He6bDRYMgbyXme1rZqvjb0xmdpKZLQqnrzazh8zsvrCur5vZrmY2I9zWCjP7VovidjKzeWE5j4YjaEbbPj7s3vos/Jvs0Up9Ww6dfFC4nz8LyzwntvrzBIO2xR0NzO7mr1vaQIGfbGMIhr6N245guNyvAOcBt5nZwDZs80hgAjCJ4BuidxKMSLkDsBdwOoCZ7UPwjc7pBEMA3wE8ZmaVsW2dThBCA9IMdXBpuPxooB/wzwQDlQHMJxihchDBN0cfslbGpG+jbwEHA7sCAwhGmlzn7vMJvhH7zdi6ZxB8gzdyXPh4ILAQmEPw//MV4BqC1x93VviahhOM7XMzgJntSjAcxMXAUIIW+OPxN9J0zGwE8ATBcMhDCfbPotgqbwHjYutHo2Iu7OavW9qgywe+md0dtiYWt9P2ngxbCv/VYv5vzezVsPXzsAXjsnR1A4BNLebVA9e4e727zwY+B9pyC79fuvtGd38DWAw8FQ49vIHgH298uN75wB3u/oq7N7r7PUAdwRtF5GZ3X9FKl9N3gB+5+zseeDUckgB3v8/d14VDC98AVLbxNbSmHugL7E7wLfO3wu4QCIYOPgNS49n/E8GbTeRFd58TvnE9RBA+v3D3euABYKSZDYitf6+7L3b3GuDHwLfty+GM/+zuT4fP/VegF5Ctn30q8Iy73x/+bde5ezz4NhEcD5GjgSc9+Cp9d37d0gZdPvCB/yBodbaX62k+XnjkEncf5+5jgeVAuptNdDWfEvwjx61r0ZrONLxvOvFRJmvTPI62tSNwmcXGsif4FDA8tn6mYYl3AJamW2BBl9RbYbfAZwSfWIZkq7iZPWFfjskzteVyd3+O4OYotwGrzexOM+sXLr4POC58o/82QdB9FHt6y/2wNhy8LHoMzfdz/LUvA8rD19Bs6GR3bwrXzTZ0cqv7K9QX+Cz2ONV/381ft7RBlw98d58LrI/Ps6Bf+UkL+ndfNLPd27C9Z9m6VYy7bwy3bQQtjyQMMvQawcf0XGUbirctVgA/azGWfXV4s49Ipn24Atip5UwL+uv/hSB8Brr7AGADwbjzGbn7UeGdofq4++9bWedmd58AfJVg310Rzv8H8DdgCkGD4N50z2+DHWLTIwha2WtpMXRyeLztQPYhjdPur5g9gFfDbZYT3P3q6WhhN37d0gZdPvBbcSdwUXgAXw78pj02amb/TnDDiN0Jb1fXxc0m+MfO1SLgRDOrtmDc9vMKKHsm8F0z2y88KdjbzI4xs5afOFpzF3Ctme0SPn+sBbcA7EvQ97sGKDOzqwj6+AsWnqTcLwzEGoKRNRtjq/yO4LzFGKCto3a2dIaZ7Wlm1QR93Q/HhjM+xswOD+txGUFX2H9n2d7vgSPM7NsWnAAfbGbxO3EdQtDlBsF9cF+LNWK68+uWNkhc4IcfPQ8gOJG3iOCk0bBw2YlmtjjNz5xctu3u5xJ89HyLoM+xq/sdcLSZ9cpx/V8DWwg+pt9D8M+UF3dfQNCPfytB19ISWr/3azo3EoTAU8BG4LcEn6zmEATXuwRdAJtp+x2rWtOP4I3q03Db6wj6kiOPELRCHwn7oAtxL0F35MdAFfB9AHd/h6DP/BaClu9xwHHx6+fTcfflBN00lxF84l1EeJI2PEG7J/Cf4eotL8fslq9b2i4RwyOb2Ujgv9x9r7Dv8R13H1bA9g4FLnf3Y1tZfghwRWvLuxIz+znwibvf1Nl16Q7MbCkw3d2f6ey65MrMbgCWuvtvwsdvAie7+5tt2EbiXre0XeK+EOPuG83sAzM7xd0fCvsCx7r7q/luM9zGTu6+JJw+juCLLF2eu1/Z2XXoLszsJILzDs91dl3awt0vi6bDyxx/18awT+Trlrbr8i18M7sfOJTgTP9q4CcEB+b/I+jKKQcecPdrctzeiwR99H0IPtqeR3By60WCj75GcPLre1EfqHR/ZvY8QbfIme6eUxdgd9BTX3dP1eUDX0RE2kfiTtqKiEh+unQf/pAhQ3zkyJGdXQ0RkcT4+9//vtbdh6Zb1qUDf+TIkSxYsKCzqyEikhhmtqy1ZerSERHpIRT4IiI9hAJfRKSHUOCLiPQQCnwRkR5CgS8i0kMo8EVEeggFvkieXnvtNWbMmNHZ1RDJmQJfJE+XX345v/jFL/jiiy+yryzSBSjwRfJUWxvc0nXz5s2dXBOR3CjwRfK0ZUtwsyYFviSFAl8kT/X19QDU1dV1ck1EcqPAF8lTFPhRS1+kq1Pgi+RJLXxJGgW+SJ4U+JI0CnyRPDU0NAAKfEkOBb5InqLAVx++JIUCXyRPjY2NgFr4khwKfJE8qYUvSaPAF8lT1MKPTt6KdHUKfJE8qUtHkkaBL5KnpqYmQF06khwKfJE8RS18Bb4khQJfJE9RC199+JIUCnyRPKlLR5JGgS+SJ3cHFPiSHO0S+GZ2pJm9Y2ZLzOwHaZafY2ZrzGxR+POd9ihXpDNFgR9djy/S1ZUVugEzKwVuA74JrATmm9lj7v5mi1VnufuFhZYn0lWohS9J0x4t/K8BS9z9fXffAjwATG6H7Yp0aVHg66StJEV7BP5XgBWxxyvDeS2dZGavmdnDZrZDO5Qr0iWohS9J0R6Bb2nmeYvHjwMj3X0s8AxwT6sbM5tmZgvMbMGaNWvaoXoi7S+6QgfUhy/J0R6BvxKIt9i3B1bFV3D3de4eff98JjChtY25+53uPtHdJw4dOrQdqifS/uI3LleXjiRFewT+fGAXMxtlZhXAacBj8RXMbFjs4fHAW+1Qrkin+fzzz1PTCnxJioKv0nH3BjO7EJgDlAJ3u/sbZnYNsMDdHwO+b2bHAw3AeuCcQssV6Uw1NTWpaXXpSFIUHPgA7j4bmN1i3lWx6RnAjPYoS6QriAe+WviSFPqmrUge1MKXJFLgi+RBLXxJIgW+SB6++OKL1HQ0TLJIV6fAF8lDPPDVwpekUOCL5KG2tjY1rT58SQoFvkge1KUjSaTAF8mDWviSRAp8kTwo8CWJFPgieYiPpaPAl6RQ4IvkIR746sOXpFDgi+Shrq4uNa3Al6RQ4IvkQV06kkQKfJE8xFv48ZuhiHRlCnyRPKhLR5JIgS+SBwW+JJECXyQPUeCbmQJfEkOBL5KHLVu2AFBaWqrAl8RQ4IvkIQr8kpISnbSVxFDgi+QhCnx16UiSKPBF8hCNgW9mauFLYijwRfJQX1+PmalLRxJFgS+Shyjw1aUjSaLAF8mDWviSRAp8kTw0NDRQUlJCSUkJ7t7Z1RHJiQJfJA/xwFcLX5JCgS+Sh8bGRkpKSnSVjiSKAl8kD1ELv7S0VF06khgKfJE8NDY2Ulpaqi4dSRQFvkge4oGfawv/+eefp7y8nOeee66DayeSngJfEuWRRx7hpptu6uxq5BX4DzzwAA0NDTzwwAMAfPHFFx1ZRZGtKPAlUaZPn85ll13W2dWgqamJsrIySktLc+7SefPNN1O/586dS+/evenbty/Tpk3js88+68jqigAKfEmY9evX09TUxAcffNCp9WhqampzC//DDz8EYPny5cyaNQuA2tpaZs6cyaBBg5g9e3ZHVVcEUOBLgixfvjw1jMETTzzRqXWJt/BzDfy1a9emfv/P//wPELyBXXvttbg7M2fO7LD6ioACXwqwbNkyli5dWrTynnzyydT0X//616KVm05TUxPl5eU5B35TUxO1tbVA0KpfunQpVVVV9OvXjyuvvBKAlStXdmidRRT40ibjxo3jgAMOAGCfffZh3LhxRSv7pZdeSk2//vrrRSs3HXenrKyMkpLc/oWi/vuKigoA1qxZw7bbbgsEN1EpLy/nk08+KbheL730EvPnzy94O9I9KfC7qFdeeSV1k418LF68mCFDhrRrv/A777zDa6+9xt/+9jfmzp3L+vXrqampKdplhlHI9+7dmxUrVhSlzNa4O+Xl5ZSVleXUwp87dy4Ae+21V2reHnvskZqurq5ulxO3Rx55JAcffHDqRPKCBQvo06dP6pxBpKGhgUMOOYQ5c+YUXKYkR7sEvpkdaWbvmNkSM/tBmuWVZjYrXP6KmY1sj3K7qyeeeIJJkyax9957572N888/n3Xr1nHBBRe0W72uvfba1PSUKVNS0zfccEO7lZHJ8uXL6d27N8OHD2fDhg1FKbM17k5FRQWlpaU5rR+1uk899dTUvEmTJqWmBwwYUPBlmhs3bqSmpobNmzdzzz33APDd736Xmpqara5s+uMf/8jcuXOZOnVqQWW2p4aGBu666y5uvPHGVPfWxo0bM+6XzZs3t7otScPdC/oBSoGlwGigAngV2LPFOhcAt4fTpwGzctn2hAkTvCcaNmyYAw74o48+2ubn/+Mf/0g9H/CXXnqpXeo1cOBALy8v95KSEge8rKzMe/Xq5X369GmX7WdTUlLiO++8sx9//PEO+JtvvlmUctMBfNKkST5hwgQP/o0ymzRpkgO+YcOG1N9l3rx5qeX77ruvA97Y2Jh3ne67777UtnfddVdfsWJFs+PgL3/5S2rdE088MTX/jjvuyLvM9nT66aen6lReXu6PP/64l5eXe2VlZdq/9U9/+lMH/NBDD202/7777vOSkhKfPHmy19bW+vDhw3369OleX1/v2223ne+9995eV1dXUF3PPvtsHzVqlL/99tsFbefDDz/0k046yefPn1/QduKABd5Kppa1w3vG14Al7v4+gJk9AEwG3oytMxm4Opx+GLjVzCysXLs75phjqK+vZ+TIkdx+++188MEH3HDDDdx0002sXbuWCy+8kM8//xwIWla33347FRUVTJ8+nTVr1nRElXK2ceNGPvroI8aPH8/ChQuZOnUq+++/f5u28c477wBwxRVXcP3113PCCScwfvz4gurV0NDAp59+yiGHHMLGjRtZuHAh++yzD/379+fpp5/msMMOo7y8vKAyspXf1NTEbrvtxqRJk3jssceYM2cOe+yxB9OmTUtd8tjeKisrue222xgxYgQQnCyOPulUVFRsdfOTW265hTFjxnDooYdy4403pk40L1q0iPLycvr160dFRQX19fVMmDAh9bwddtiB+fPn895777Hbbrs12+bixYv5/e9/z3XXXZea96tf/Yrx48fzzW9+k1tuuYXtt9+ep556CoD+/fvz7rvvcuCBBwJw5ZVX8vOf/5yTTz6Z/fbbj5tvvplXXnkldUnplVdeybRp05g/fz4//vGP8xoq4uSTT2batGlA8LeaNm0aK1eu5KCDDuKqq65K7beSkhKuueYaJk6c2Oz5TU1N/PGPf6RXr16ce+65/OY3v+G4445LLZ8wYQJvvPEGgwcPZvr06axatSrVTfb888/zrW99i6eeeooHH3yQM888E3fn0Ucf5YQTTmDVqlXceeedrF+/no8//piPP/6YESNGMHbsWM466yzOOOMMbrrpJiorK/ne974HwGeffcYll1zCj370I4YNG8b555+fyoalS5fy/vvvA0EX3f3338+JJ57IZZddxuGHH86xxx7LjBkz2GuvvZg6dSpXXXUVL7/88lb7rKGhgblz59LY2Mif/vQnDjjgAKqrqwHo06cPf/rTn9r8d8iqtXeCXH+Ak4G7Yo/PBG5tsc5iYPvY46XAkFa2Nw1YACwYMWJEvu9wqZ8DDjjA+/btm2qRbbPNNs2WA77LLrv4uHHjtprfWT9VVVW+evXqZi2etv5ss8027u7+1a9+tV3r9swzz/hDDz3kZuYPP/ywP/vss0XdN3fffbfPmzfPAT/llFP85Zdf7vAyp0+fnjq29ttvv9T8Cy64wPfff3+HL1v4JSUlvssuu7i7e3l5ebPt7L777u7uPm7cOB89enSzY/byyy93wGfNmrXV8Rwdmy+++KK7u69evdoBHzRokG/YsMHNzHv16pX6W8db+kOHDnV39zFjxqTmjR071s3MR48e7UcffbQDfvvtt/u2225b0H56/PHH3d39mGOOaTb/kksuabYvKisrffXq1c1e4w033OCAf+c733F39wsvvNDNzK+99tpUS753797N/n979erlr7/+uo8aNcoBnzBhgpuZl5aW+owZM1LrVVRUpKarqqr8qKOOala/vfbaKzU9ZcoUnzVrViozKioqfPDgwVu91vHjx/sdd9yR+rQ7aNCg1LLtttsuNT1w4MCM+6yystKvuOIK79+/f7P5JSUleWVfmH+ttvDbI/BPYevAv6XFOm+wdeAPzrbtfLt0Nm3a5Bs2bPCdd945tQP79euXmj7xxBN906ZNvmnTJj/ooINS88eNG5ea35k/9fX1qddSU1OT1zbiXQPtVa/WPgbX1dUVZb/U1ta6u3tjY2Pq73Xdddc54Nddd127l/fCCy844GeeeWbqtY4fP94Br6mpcXf3r3/9680CH/AhQ4akpseMGZP2b9LSzJkzHfCrr756q2VRYJ100knu/uWbA+BTp05NTZuZ9+/f393da2tr0x4HO+20U2r9s846y1evXu1mlirj8MMPb/N+mjdvnpuZV1VV+TnnnOOADxs2zJctW+ZVVVWp8n7961+nwrt///6+0047+ejRo32XXXbxqqoqNzP/9NNPU/WN1z0e4Keffnqz11ZbW+s77LBDKiifeuopd/fUvJkzZ/ro0aOb7d/a2lpfvHix9+nTJ/XG2PIN79BDD222r6LXG/3t3d0XLlyY2sZhhx2WeqPYb7/9UuUffPDBvm7dujb/n+arowN/f2BO7PEMYEaLdeYA+4fTZcBawLJtu9A+/OgAP/XUU33Dhg0+dOhQHzduXLOdXFdX58OGDfNBgwYVtJOluHr16uWDBw/2k08+2QFfuHBhu5fx+uuvO+Cnnnpqat7YsWObBXwUCvX19V5TU+OAV1dX+6effuqAH3HEETmVFX1SmTp1atr5USvS3X3kyJHNgsnMmrU8M3nyySdT686ePdvdv2yRm5l/9NFHOdW3pZ/85CfNWqdR3/aTTz7ppaWlPmXKlNS6U6ZMccBLS0u9rKzMS0pK3Mz82GOPzVjGL3/5S7/yyivTLqupqfHJkyenwt7d/e233/af/vSn7u6+YsUK/8lPfrLVm+7q1av9iiuu8JqaGq+rq/Orr77aL7rootSnlRdeeMHvvffejPWqqanxl19+OTX9zDPPpJatWLEi43M7QkcHfhnwPjCKL0/afrXFOv+b5idtH8xl28U8aVvIyTIpvhEjRnhZWVkqgDvi77dkyZLUx/zIXnvt5eHpJ3d3P+KIIxzwTZs2pdYvLS1NBXX800Em0ZvFwQcf3Gx+1GIeOnSoA/7uu+86kGoVA77nnnumWpPx7qfWjBo1yktLS1P7bPXq1V5dXe3nnXdeTnVtzcKFC/2hhx7aKuT0v1VcmQK/4JO27t5gZheGrfhS4G53f8PMrgkLfgz4LXCvmS0B1oeh36Xk+gUa6RpGjRrF8uXLWbZsGdXV1R3y9+vVqxdAs+9DtDxJW1ZWllrn448/Tq3z3nvvAbD99tvnVFb0GqJtRJ577jlKSkr48Y9/zPe//3323XdfAE477TTmzZvHnDlzmDZtGuXl5Vx44YWcccYZWctatGgRa9asSe2zbbbZhpqampzqmcnee++d9lJi/W91He1xlQ7uPhuY3WLeVbHpzQR9/SLtYsyYMbzwwgts2LCBHXfcsUPKiK6YqK+vT81ramrCzFKPo+vwN2/e3Cyso7Fy2lK3qqoq1q1bx+bNm7nrrrv4wx/+wPLlyxk9ejTnn38+l156KRs2bGDAgAFceumlNDU1cf3113PRRRdRUlKS83cu+vXrR79+/XKul3Qf7RL4IsU2adIkbr31VgBGjx7dIWW0Fvhx0aWodXV1zYZGiAJ/5513zrm8fv368fHHH1NdXR11hdK3b19uvPFGqqqq2LBhA2VlZanhGYBml2qKZKPPWpJIhx9+eGp6zJgxHVJGFKzxwG9sbGzWwo+6dDZv3tzsOxxLliwB2Oqa+kwmTpyImTFy5Eguvvhili1bxsaNG5k8eTIQvAHFw16krdTCl0TabrvtUveTjQZz6yiZunSiFn59fT3r1q1LzY/Cf/jw4TmX8/jjjxdaVZGM1MKXxOrfvz8Ahx12WIeWEx+Xxd3TtvDr6uqaDX62ZcuW1A1SRLoKHY2SWGPHjqW6upptttmmw8ows4wt/NYCH1D3i3Q56tKRxCrWsMzxFn7LwI9CfcuWLVuN4Bmd9BXpKtTCl8QqKSnp8C4TM8sY+NFlmXV1dWzatAn4stUfdTmJdBUKfJEMsgV+vIX/+eefU1JSQlVVFQADBw4sbmVFslDgi2RQUlLS7Nu1rV2lU1dXR21tLWVlZfTu3RuAoUOHFreyIlko8EUyyNbCj7pv6uvrqa2tpby8PNWVM2zYsOJWViQLBb5IBi1b+C0vy4x36dTV1VFRUcHgwYMBUjdNEekqFPgiGaTr0omfKI66dLZs2cKWLVuoqqpi2223BYIB3kS6EgW+SAbZWvjxb9o2NDRQXV3NQQcdBDS/SblIV6DAF8kgGr4h0rKFH+/SaWpqonfv3lx22WXU1dW1aRwdkWJQ4ItkUFpamlMffl1dHe5O3759m80X6UoU+CIZtGzhu3vaFv769esBNM68dGkKfJEMSktLM3bpRH340Vj4+natdGUKfJEMSktLUzcjgdZb+GvXrgX07Vrp2hT4IhmUlZVt1aWTrg8/6tKJrsEX6YoU+CIZZOvDj7p0opEyhwwZUtwKirSBAl8kg7KysoxdOpWVlQCpkTI7cmx+kUIp8EUyyDXwa2pqAAW+dG0KfJEMWgY+fDkGPnzZpbN582YguNeuSFelwBfJINcWfm1tLYC+XStdmgJfJINo+ONIy8CPbnYCQfi3XF+kK1Hgi2QQddnEpWvhAwwYMKAodRLJlwJfJIMo8KNLM929WR9+fMyc4cOHF7dyIm2kwBfJID7efaS1Fv5OO+1UvIqJ5EGBL5JBFPhffPEFsHULP96HP2bMmOJWTqSNFPgiGbQMfKDVLp399tuveBUTyYMCXySDKNCjyy6BtIOnARx44IHFq5hIHhT4IhmkC/x4Cz+6DLO0tJQ+ffoUt3IibaTAF8kgW+BHdOMTSQIFvkgG8cCPLs1MF/jbbrttUeslkg99LVAkg3jgNzQ0AFsH/gUXXMAxxxxT9LqJtJUCXySD6Dr7zZs3p67Fbxn4t912W9HrJZKPgrp0zGyQmT1tZu+Fv9Pe383MGs1sUfjzWCFlihRTLoEvkhSF9uH/AHjW3XcBng0fp1Pr7nuHP8cXWKZI0SjwpTspNPAnA/eE04VsUMIAAAwGSURBVPcAJxS4PZEuJV3ga0RMSapCA39bd/8IIPzd2u1+qsxsgZm9bGYZ3xTMbFq47oI1a9YUWD2RwkSBX1dXR319PaAWviRX1qaKmT0DpLuNzw/bUM4Id19lZqOB58zsdXdfmm5Fd78TuBNg4sSJnm4dkWKJxsqJt/Dj37QVSZKsge/uR7S2zMxWm9kwd//IzIYBn7SyjVXh7/fN7HlgPJA28EW6knQtfHXpSFIV2lR5DDg7nD4beLTlCmY20Mwqw+khwIHAmwWWK1IUUQu/rq5OJ20l8QoN/F8A3zSz94Bvho8xs4lmdle4zh7AAjN7FfgL8At3V+BLIsQDXy18SbqCjlx3Xwccnmb+AuA74fR/AxooXBKpV69eQHADFJ20laTT2SeRDKI+/Pr6erXwJfEU+CIZxLt0orF0FPiSVAp8kQzSdeko8CWpFPgiGSjwpTtR4ItkUF1dDagPX7oHBb5IBmrhS3eiwBfJIN7C10lbSToFvkgG6tKR7kSBL5JBNFBaPPDLy8s7s0oieVPgi+RAXTrSHSjwRXLQ0NCgFr4kngJfJAszo6GhQS18STwFvkgWZqY+fOkWFPgiWZgZjY2NauFL4inwRbJo2aWjFr4klQJfJAsFvnQXCnyRLEpKSpp16SjwJakU+CJZKPClu1Dgi2ShwJfuQoEvkoUCX7oLBb5IFiUlJTQ1NSnwJfEU+CJZlJaWqoUv3YICXySL0tJSmpqaaGxsBKCioqKTaySSHwW+SBZR4KuFL0mnwBfJorS0FHdX4EviKfBFslCXjnQXCnyRLKIWvgJfkk6BL5JFWVmZ+vClW1Dgi2RRVlamFr50Cwp8kSwqKysV+NItKPBFsogCXoEvSafAF8miqqoKINWHr8CXpFLgi2RRWVkJwObNm4FgbB2RJNKRK5JFFPi1tbWdXBORwijwRbLo1asX8GULXySpFPgiWSjwpbsoKPDN7BQze8PMmsxsYob1jjSzd8xsiZn9oJAyRYqturoagLq6uk6uiUhhCm3hLwZOBOa2toKZlQK3AUcBewKnm9meBZYrUjRRC7++vh4z6+TaiOSvrJAnu/tbQLZ/gq8BS9z9/XDdB4DJwJuFlC1SLFHgb9mypZNrIlKYYvThfwVYEXu8MpyXlplNM7MFZrZgzZo1HV45kWz69OkDfPnFK5GkytrCN7NngO3SLPqhuz+aQxnpmv/e2srufidwJ8DEiRNbXU+kWKI+/IaGBnXpSKJlDXx3P6LAMlYCO8Qebw+sKnCbIkUTtfCbmpoU+JJoxejSmQ/sYmajzKwCOA14rAjlirSLeOCLJFmhl2VOMbOVwP7An81sTjh/uJnNBnD3BuBCYA7wFvCgu79RWLVFiqdv376pabXwJckKvUrnEeCRNPNXAUfHHs8GZhdSlkhniVr4oMCXZNM3bUWyUOBLd6HAF8mif//+qWkFviSZAl8kC/XhS3ehwBfJol+/fqlpBb4kmQJfJIvojlegwJdkU+CLtIECX5JMgS+SgyjoFfiSZAp8kRxEQa/72UqS6egVyUEU9GrhS5Ip8EVyEAW+WviSZDp6RXJQWloKqIUvyabAF8mBWvjSHejoFclBWVkwzqACX5JMR69IDhT40h3o6BXJQRT46sOXJFPgi+SgvLwcUAtfkk1Hr0gOFPjSHejoFcmBAl+6Ax29IjmorKwEvrweXySJFPgiOaioqADUwpdk09ErkoOoha/AlyTT0SuSAwW+dAc6ekVyEN31Sn34kmQKfJEcRIGvFr4kmY5ekRyohS/dgQJfJAfV1dWAAl+STYEvkoNevXoB6tKRZNPRK5IDtfClO1Dgi+RAgS/dgQJfJAcKfOkOFPgiOejTpw+gwJdkU+CL5KB3796AAl+STYEvkoO+ffsCCnxJNgW+SA7UpSPdgQJfJAf9+vUDvry3rUgSFRT4ZnaKmb1hZk1mNjHDeh+a2etmtsjMFhRSpkhniAJfLXxJskKbK4uBE4E7clj3MHdfW2B5Ip1it9124xvf+AYXX3xxZ1dFJG8FBb67vwVgZu1TG5EuqqSkhGeffbazqyFSkGL14TvwlJn93cymZVrRzKaZ2QIzW7BmzZoiVU9EpPvL2sI3s2eA7dIs+qG7P5pjOQe6+yoz2wZ42szedve56VZ09zuBOwEmTpzoOW5fRESyyBr47n5EoYW4+6rw9ydm9gjwNSBt4IuISMfo8C4dM+ttZn2jaeBbBCd7RUSkiAq9LHOKma0E9gf+bGZzwvnDzWx2uNq2wEtm9iowD/izuz9ZSLkiItJ2hV6l8wjwSJr5q4Cjw+n3gXGFlCMiIoXTN21FRHoIBb6ISA9h7l33ykczWwMsy/PpQ4Cu+M1e1attVK+2Ub3apjvWa0d3H5puQZcO/EKY2QJ3b3V8n86ierWN6tU2qlfb9LR6qUtHRKSHUOCLiPQQ3Tnw7+zsCrRC9Wob1attVK+26VH16rZ9+CIi0lx3buGLiEiMAl9EpIdIROCb2ZFm9o6ZLTGzH6RZvqOZPWtmr5nZ82a2fWzZ2Wb2Xvhzdmz+hPC2i0vM7GbL4y4u+dbLzPY2s7+Ft4d8zcxOjT3nP8zsg/B2kIvMbO9i1Stc1hgr+7HY/FFm9kq4H2eZWUWx6mVmh8XqtMjMNpvZCeGy9thfd5vZJ2aWdlA/C9wc1vs1M9sntqwjj6+86lWE46uQ/dWRx1e++6uzj6/dw79XnZld3mJZ2v+ZvPeXu3fpH6AUWAqMBiqAV4E9W6zzEHB2OP0N4N5wehDwfvh7YDg9MFw2j2DQNwOeAI4qYr12BXYJp4cDHwEDwsf/AZzcGfsrfPx5K9t9EDgtnL4d+F4x6xVbZxCwHqhuj/0VbuNgYB9gcSvLjw6PEQMmAa909PFVYL067PgqpF4deXwVWq9OPr62AfYFfgZcnsv/TL77Kwkt/K8BS9z9fXffAjwATG6xzp5AdP+5v8SW/xPwtLuvd/dPgaeBI81sGNDP3f/mwR77HXBCserl7u+6+3vh9CrgEyDtN+PyUMj+SitsnX4DeDicdQ9F3F8tnAw84e5ftLH8VnlwM571GVaZDPzOAy8DA8JjqCOPr7zr1cHHVyH7K612Or7aq15FP77c/RN3nw/Ut1iU9n+mkP2VhMD/CrAi9nhlOC/uVeCkcHoK0NfMBmd47lfC6Uzb7Mh6pZjZ1wjevZfGZv8s/Mj5azOrLHK9qiy4xeTL0cdaYDDwmbs3ZNhmR9crchpwf4t5heyvXGQ6jjrq+CqkXikdcHwVWq+OOr4KrVekM46v1rRW37z3VxICP13fZ8trSS8HDjGzhcAhwD+AhgzPzWWbHVmvYANB6+Je4Fx3bwpnzwB2J/iINwj4lyLXa4QHX+n+X8BNZrZTjtvs6HpF+2sMMCf2nEL3Vy7aehy1x/7KRcZyOuj4KrReHXV8FVqvzjy+WtPux1cSAn8lsEPs8fbAqvgK7r7K3U909/HAD8N5GzI8d2U43eo2O7hemFk/4M/Aj8KPl9FzPgo/ctYB/07wsa5o9fIvb0f5PvA8MJ5gEKcBZlbW2jY7ul6hbwOPuHt97DmF7q9C6t6Rx1ch9erI46ugenXg8VVQvUKddXy1prX65r2/khD484FdwrPSFQQfuR6Lr2BmQ8wsei0zgLvD6TnAt8xsoJkNJLi94hx3/wjYZGaTwv6ws4Bcb8hecL3C9R8h6E98qMVzhoW/jaBfrq23gyykXgOjj6xmNgQ4EHgz7If+C0H/JsDZFHF/xZxOi4/b7bC/cvEYcFZ4lcckYEN4DHXk8ZV3vTr4+CqkXh15fOVdr9jyzjq+WpP2f6ag/dXyLG5X/CE4u/4uQT/kD8N51wDHh9MnA++F69wFVMae+8/AkvDn3Nj8iQR/vKXArYTfOi5GvYAzCE7QLIr97B0uew54PazbfUCfItbrgLDsV8Pf58W2OZrgypMlBFfTVBarXuGykQRdPCUtttke++t+gitZ6glaVecB3wW+Gy434Law3q8DE4t0fOVVryIcX/nWq6OPr0L+jp15fG0Xzt8IfBZO92vtf6aQ/aWhFUREeogkdOmIiEg7UOCLiPQQCnwRkR5CgS8i0kMo8EVEeggFvohID6HAFxHpIf4/T4d8DxrSSaoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for n in range(2):\n",
    "    plt.plot(p_vec, (x_num - x_sym)/x_sym, 'k')\n",
    "    plt.grid()\n",
    "plt.title(label=\"Error in solution\\n(numerical - symbolic)/symbolic\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#   x1 + 2*x2 + 3*x3 = 7\n",
    "# 4*x1 + 5*x2 + 6*x3 = 8\n",
    "# -----------------------"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_vars = sympy.symbols(\"x_1, x_2, x_3\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = sympy.Matrix([[1, 2, 3], [4, 5, 6]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = sympy.Matrix(x_vars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = sympy.Matrix([7, 8])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATQAAAAWCAYAAABOkmAWAAAABHNCSVQICAgIfAhkiAAAB/xJREFUeJztnGmMFEUUgD8QRQXxFo2oq6BBRTHeiUoG7yOoePyQaDRKvFAwHvGIholGFw8EBS9icELERNF4nxwh8T5RURAVWEUOFdEFQRF1/PFeO7291V1dPVPsjPaXTHa2qt7r1/WqXlVXVQ/k5OTk/E8YArwDtAJ/Aft1rDk5OTn/U4rAOmAZ8Aywh6uCvYC/gVXAROAWYPva2ZcT4gxgHPA6sBIoA49aZDoB5yMDzipgDTALGA5skCDXExmc7g2l3Q5MBxYBvwErVNdIYGu3W2lItgaGAk8DXyN10Aq8AVwAdE6Q7YX0jyXAWqAFGAtsablmI/ihmnqpNQVgFPAK0j/ec1UwTAWH19SsHBMfI3W9CphLuoA2Sct9DzwM3AN8rmlPIgHPxIVaZmAo7Q8kME5EGs044H0ttxjYyfWGGoyLkXtdAkwGmpG6+IXk+uyN1H8ZmTWMAmbo/1+QHIQawQ9Z68U3HyKTre4uQjchRh/jw6KcNgwEdkcaRwF7QDtVyywAtgmlb4iMpmXgvBjZl4HltJ3FbRxT9lbVdX+S8f8BjgQG0X7GsT3wLVIHpxvkXtW8yyPpd2v6gwnXbAQ/ZK0X30zWazsF+KIKFRLKvKZlToukdwJKmjfK5aJVUm/2ZKGAPaAFs7Nhhrx+mvehIW9z5LHokZS29FddU1OWd6UR/HWD2jEukr6bpi+kfYffDPgVWA10M+isNz9kIa5ewL9fAx1NLkJF7AGtP7IOMIe2I81olZ3gcsEYSiTPONa3Pb4pYA9oQYM5wZC3qeaVgS0ieUM0fVBKW27U8qNTlnel3tqPiWtUfkwkfaimPxQjF8zejjLk1ZsfshBXL+DfryUyBLSxKnRISuXn6f9B5H6c2iwaRvW7lq+1Pb4pYA9oj2mZSw15wQytDBwayZuCzBziHm2uRgayMcgGRRn4BNg2leXZKFFf7SdMF2C2yh8XybtT06+KkR2v+ZcY8urRDy4k1UtACX9+fVB19U8r0InKwpttZ7MXsvvRAlymF3oF2CiDoSZ2APoi0/Q0+LbHNwXsAS0Y4b8GtgqldwGeohLQwjO4jZFNhykJepeFZMvIOk9PJ+vdqbf2E+YutedFQ94EzRsaIxuse10fSa9XP7iQVC8BPv16neobYSt4PDItDHZWkhY1wzRTqfw3kceeKAOA55EdkzJyVMEXaewZBnyKHJNYCbwNnJRCdwttG5vtY9utjFJIIdcZeEnLLUM611jgM6QRfal5x4ZkTta0ISls6AkMBuYh/trf5QYyYPPX9UibXAn8iLSjfp5tGq72zKXtoBFgC2i3af51kXSffmjBb9sEe72Esfk1ax/cEbnXdcATyODRBDKihzmeStT7DPOCn4kfQ98vQM5ERemGTJsfQWYRPkljz3dIY/sKmZGei2y9H4BUchzzgd8dbFniUDYtfyMdYwRwjn7WAW8h9zEe2TX9ISQzGDkWkDSqBnyP7JZ+hATHSfgNIDZ/FZAdvvcRX90MTEPOSq7wYM8w5BjMHGQNzHSNVv0bN/PrESkX4NMPvttmmnoJY/Nr1j64GAmWY4AzNW0qEuTa0QM4BRnpvyP5kCbAWUgHW4pE4gcs5cHvDC2LPQErgIt8GOVAgeyjJ8AmiO/WIMc4QHy4HHl0cWWW2rONrWBGsvirO7LonHZR3YUr1I7ZwHYJ5bJsCtSzH2ykrZeArP0wTR88SXW+AexNpZ0n8qgKJY0IJyKjzWxkwXIuMlPoa9HtK6BltWcDxAF/APt4sMuFAtUFtODAZimUNlDTLsygLzg4ajv5noWs/tpBbTqsxvZcq3pnYQ8cvbXsQuKPbayh7bGNevWDDZd6gWx+demDwW7p0Sls+ZdgUXNATP7hiMMWIA0MJEiVkWljEi4BLe2ibhZ79kEa3p/Az4gjOpoC6QJaD0PaQcgItwo5JxUwDpnRmBaW+2Le+OlMpQ28acgvUd1xiGraz+NI57I9PUD69hMcJP8A+9pQgOvBWh9+8I1rvbj6NUsfnKj6+qQo+y9FFSoY8vojrz8sRUaqMMGGwhEJul0CWgl7x8lqz0ZIpRyIPJMvx/9is4lTkfssUXlXbX4o7S6DzLvATGS9rBl4DmkUq2m/lb4I2f43cQUyek5HFrqD11vmqx1LkbWqKMHh3rOTb81INe3nTpVL25hL2NvPuVrmT2Rdpmj4mOSjrz41U3n1aR7tX33y4QefuNZLFr9m6YMl1dWU6i6UkZgDWh9kZ+1nYF+D3NEq906C7loGtFrYEzANeS9yfVMkeTeqxSBzDXK05hfk1PlCZEbQFCl3kOq4Muba/YD7kPdJlyONtxVpgEXiR+VZyM6U6yNQNf4ajQQQl45dwh7Qith3BGfGyO6EbHQtRR6XvkEWzqP15ssPPimSvl5q1Q/T9MGS6trFUq4NwVmPwS5CKfG5KVANM8i+dlWvBMcHdq2hzi2QR6c7aqjTxr24B7N6wocf/ouk6YPPInXpdDbvdBV6ARmFqj3d2x35PbX9VO/V+n3nKvVmZRQy/W1CnuObkV0Z0+tEjcxcZNSvJYOQ4wHr6+ek7kdmg0fqNYOP068tdDA+/NDouPbBrsDByMy1FceYtAmV1xuCTzU/8FjAPGUtVaGzGkrI48Fa5LzWNOJf48jpWOIed4odaFNO9ZRI3weLtPX9SFMh228ZdUVGxT2R0XAC8pyck5OTsz4pIDuoPyFvFOSz3ZycnJycnJychuAf9QVDryMPQ1AAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$\\displaystyle \\left\\{ x_{1} : x_{3} - \\frac{19}{3}, \\  x_{2} : \\frac{20}{3} - 2 x_{3}\\right\\}$"
      ],
      "text/plain": [
       "{x₁: x₃ - 19/3, x₂: 20/3 - 2⋅x₃}"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sympy.solve(A*x - b, x_vars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a11=p+1\n",
    "a12=1/p\n",
    "a21=p^2-1\n",
    "a22=abs(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "sp.init_printing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = sp.symbols(\"p\", real=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}p + 1 & \\frac{1}{p}\\\\p^{2} - 1 & \\left|{p}\\right|\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡         1 ⎤\n",
       "⎢p + 1    ─ ⎥\n",
       "⎢         p ⎥\n",
       "⎢           ⎥\n",
       "⎢ 2         ⎥\n",
       "⎣p  - 1  │p│⎦"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sp.Matrix([[p+1, 1/p], [p**2-1, sp.Abs(p)]])\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "[1  1]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = sp.Array([1, 1])\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'list' object has no attribute 'rows'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-17-57fc85bd3ce3>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mA\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\sympy\\matrices\\matrices.py\u001b[0m in \u001b[0;36msolve\u001b[1;34m(self, rhs, method)\u001b[0m\n\u001b[0;32m   4500\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mmethod\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;34m'GJ'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   4501\u001b[0m             \u001b[1;32mtry\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 4502\u001b[1;33m                 \u001b[0msoln\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mparam\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgauss_jordan_solve\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mrhs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   4503\u001b[0m                 \u001b[1;32mif\u001b[0m \u001b[0mparam\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   4504\u001b[0m                     raise ValueError(\"Matrix det == 0; not invertible. \"\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\sympy\\matrices\\matrices.py\u001b[0m in \u001b[0;36mgauss_jordan_solve\u001b[1;34m(self, B, freevar)\u001b[0m\n\u001b[0;32m   2971\u001b[0m         \u001b[1;32mfrom\u001b[0m \u001b[0msympy\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmatrices\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0mMatrix\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mzeros\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2972\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 2973\u001b[1;33m         \u001b[0maug\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mhstack\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mB\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   2974\u001b[0m         \u001b[0mB_cols\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mB\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcols\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2975\u001b[0m         \u001b[0mrow\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcol\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0maug\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mB_cols\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\sympy\\matrices\\common.py\u001b[0m in \u001b[0;36mhstack\u001b[1;34m(cls, *args)\u001b[0m\n\u001b[0;32m    379\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    380\u001b[0m         \u001b[0mkls\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 381\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mreduce\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkls\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrow_join\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    382\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    383\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mreshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrows\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcols\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\ProgramData\\Anaconda3\\lib\\site-packages\\sympy\\matrices\\common.py\u001b[0m in \u001b[0;36mrow_join\u001b[1;34m(self, other)\u001b[0m\n\u001b[0;32m    481\u001b[0m             \u001b[1;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_new\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mother\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrows\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrow_join\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mother\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    482\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 483\u001b[1;33m         \u001b[1;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrows\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[0mother\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrows\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    484\u001b[0m             raise ShapeError(\n\u001b[0;32m    485\u001b[0m                 \"`self` and `rhs` must have the same number of rows.\")\n",
      "\u001b[1;31mAttributeError\u001b[0m: 'list' object has no attribute 'rows'"
     ]
    }
   ],
   "source": [
    "x = A.solve([1, 1])\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}