Predavanje 06 - Sistemi linearnih enačb 2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<font size=\"2\"> [github.com/jankoslavic/pypinm](http://github.com/jankoslavic/pypinm)</font>\n",
    "\n",
    "<center>\n",
    "<img width=300 src=\"./fig/ladisk.png\">\n",
    "<br>\n",
    "<font size=\"7\" color=\"f00e0e\" face=\"garamond\"><b>Sistemi linearnih enačb (2)</b></font>\n",
    "<br><br>\n",
    "<font size=\"4\">\n",
    "[izr. prof. dr. Janko Slavič](http://www.ladisk.si/~slavic)\n",
    "</font>\n",
    "<br><br>\n",
    "<font size=\"2\">November 2017</font>\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "toc": "true"
   },
   "source": [
    "# Table of Contents\n",
    " <p><div class=\"lev1 toc-item\"><a href=\"#Razcep-LU\" data-toc-modified-id=\"Razcep-LU-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Razcep LU</a></div><div class=\"lev2 toc-item\"><a href=\"#Razcep-LU-matrike-koeficientov-mathbf{A}\" data-toc-modified-id=\"Razcep-LU-matrike-koeficientov-mathbf{A}-11\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>Razcep LU matrike koeficientov $\\mathbf{A}$</a></div><div class=\"lev2 toc-item\"><a href=\"#Numerična-implementacija-razcepa-LU\" data-toc-modified-id=\"Numerična-implementacija-razcepa-LU-12\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>Numerična implementacija razcepa LU</a></div><div class=\"lev1 toc-item\"><a href=\"#Pivotiranje\" data-toc-modified-id=\"Pivotiranje-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Pivotiranje</a></div><div class=\"lev2 toc-item\"><a href=\"#Gaussova-eliminacija-z-delnim-pivotiranjem\" data-toc-modified-id=\"Gaussova-eliminacija-z-delnim-pivotiranjem-21\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Gaussova eliminacija z delnim pivotiranjem</a></div><div class=\"lev2 toc-item\"><a href=\"#Razcep-LU-z-delnim-pivotiranjem\" data-toc-modified-id=\"Razcep-LU-z-delnim-pivotiranjem-22\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span>Razcep LU z delnim pivotiranjem</a></div><div class=\"lev1 toc-item\"><a href=\"#Modul-SciPy\" data-toc-modified-id=\"Modul-SciPy-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Modul <code>SciPy</code></a></div><div class=\"lev1 toc-item\"><a href=\"#Računanje-inverzne-matrike\" data-toc-modified-id=\"Računanje-inverzne-matrike-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Računanje inverzne matrike</a></div><div class=\"lev1 toc-item\"><a href=\"#Reševanje-predoločenih-sistemov\" data-toc-modified-id=\"Reševanje-predoločenih-sistemov-5\"><span class=\"toc-item-num\">5&nbsp;&nbsp;</span>Reševanje predoločenih sistemov</a></div><div class=\"lev1 toc-item\"><a href=\"#Iterativne-metode\" data-toc-modified-id=\"Iterativne-metode-6\"><span class=\"toc-item-num\">6&nbsp;&nbsp;</span>Iterativne metode</a></div><div class=\"lev2 toc-item\"><a href=\"#Gauss-Seidlova-metoda\" data-toc-modified-id=\"Gauss-Seidlova-metoda-61\"><span class=\"toc-item-num\">6.1&nbsp;&nbsp;</span>Gauss-Seidlova metoda</a></div><div class=\"lev2 toc-item\"><a href=\"#Zgled\" data-toc-modified-id=\"Zgled-62\"><span class=\"toc-item-num\">6.2&nbsp;&nbsp;</span>Zgled</a></div><div class=\"lev1 toc-item\"><a href=\"#Nekaj-vprašanj-za-razmislek!\" data-toc-modified-id=\"Nekaj-vprašanj-za-razmislek!-7\"><span class=\"toc-item-num\">7&nbsp;&nbsp;</span>Nekaj vprašanj za razmislek!</a></div><div class=\"lev2 toc-item\"><a href=\"#Dodatno\" data-toc-modified-id=\"Dodatno-71\"><span class=\"toc-item-num\">7.1&nbsp;&nbsp;</span>Dodatno</a></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Razcep LU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Za rešitev sistema linearnih enačb zahteva Gaussov eliminacijski postopek najmanjše število računskih operacij. \n",
    "\n",
    "V primeru, ko se matrika koeficientov $\\mathbf{A}$ ne spreminja in se spreminja zgolj vektor konstant $\\mathbf{b}$, se je mogoče izogniti ponovni Gaussovi eliminaciji matrike koeficientov. Z razcepom matrike $\\mathbf{A}$ lahko pridemo do rešitve z manj računskimi operacijami. V ta namen si bomo pogledali razcep LU!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Poljubno matriko lahko zapišemo kot produkt dveh matrik:\n",
    "$$\\mathbf{A}=\\mathbf{B}\\,\\mathbf{C}.$$\n",
    "Pri tem je možnosti za zapis matrik $\\mathbf{B}$ in $\\mathbf{C}$ neskončno veliko.\n",
    "\n",
    "Pri razcepu LU zahtevamo, da je matrika $\\mathbf{B}$ spodnje trikotna in matrika $\\mathbf{C}$ zgornje trikotna:\n",
    "\n",
    "$$\\mathbf{A}=\\mathbf{L}\\,\\mathbf{U}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Vsaka od matrik $\\mathbf{L}$ in $\\mathbf{U}$ ima $(n+1)\\,n/2$ neničelnih elementov; skupaj torej $n^2+n$ neznank. Znana matrika $\\mathbf{A}$ definira $n^2$ vrednosti. Za enolično določitev matrik $\\mathbf{L}$ in $\\mathbf{U}$ torej manjka \n",
    "$n$ enačb. Tukaj bomo uporabili razcep LU, ki dodatne enačbe pridobi s pogojem $L_{ii}=1$, $i=0, 1,\\dots,n-1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Sistem linearnih enačb:\n",
    "$$\\mathbf{A}\\mathbf{x}=\\mathbf{b}$$\n",
    "torej zapišemo z razcepom matrike $\\mathbf{A}$:\n",
    "$$\\mathbf{L}\\,\\underbrace{\\mathbf{U}\\,\\mathbf{x}}_{\\mathbf{y}}=\\mathbf{b}.$$\n",
    "\n",
    "Do rešitve sistema $\\mathbf{A}\\mathbf{x}=\\mathbf{b}$ sedaj pridemo tako, da rešimo dva trikotna sistema enačb."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Najprej izračunamo vektor $\\mathbf{y}$:\n",
    "$$\\mathbf{L}\\,\\mathbf{y}=\\mathbf{b}.\\qquad \\textrm{(direktno vstavljanje)}$$\n",
    "Ko je $\\mathbf{y}$ izračunan, lahko iz:\n",
    "$$\\mathbf{U}\\,\\mathbf{x}=\\mathbf{y}\\qquad \\textrm{(obratno vstavljanje)}$$\n",
    "določimo $\\mathbf{x}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Razcep LU matrike koeficientov $\\mathbf{A}$\n",
    "\n",
    "V nadaljevanju bomo pokazali, da Gaussova eliminacija dejanjsko predstavlja razcep LU matrike koeficientov $\\mathbf{A}$. Pri tem te si bomo pomagali s simbolnim izračunom, zato uvozimo paket `sympy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "import sympy as sym # uvozimo sympy\n",
    "sym.init_printing() # za lep prikaz izrazov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Prikaz začnimo na primeru simbolno zapisanih matrik $\\mathbf{L}$ in $\\mathbf{U}$ dimenzije $3\\times 3$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "L21, L31, L32 = sym.symbols('L21, L31, L32')\n",
    "U11, U12, U13, U22, U23, U33 = sym.symbols('U11, U12, U13, U22, U23, U33')\n",
    "L = sym.Matrix([[  1,   0,  0],\n",
    "                [L21,   1,  0],\n",
    "                [L31, L32,  1]])\n",
    "U = sym.Matrix([[U11, U12, U13],\n",
    "                [  0, U22, U23],\n",
    "                [  0,   0, U33]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Matrika koeficientov $\\mathbf{A}$ zapisana z elementi matrik $\\mathbf{L}$ in  $\\mathbf{U}$ torej je:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAABLCAMAAACbU3QDAAAAPFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAo1xBWAAAAE3RSTlMA\nMquZdlQQQOkwRCLd77tmzYlspWksOAAAAAlwSFlzAAAOxAAADsQBlSsOGwAACVZJREFUeAHtXWuD\noygQJOaxe5PMJDn//389HlZDt4DIwHrO6odFgS66uuIjmUpWnUa7DerYflgFHk5ZpU7j+aK36w/j\nd9BRT6PrMBqJT0c5fm4FnhmJ78Nj/Px6avJf4/i6/PkibJ5AK8qbEslJrNR7/HAsLxsIbFbePIFW\nIm9IJC/xS5/jdjvfWnFdh7N5AuvSTc/ekEhe4vFrSnqrp+3NE0iLtm5kQyJZiT/G6fp8e68j1Gr2\n5gn8BCJZid941H5Ot+RWjEtxNk+gNNGleVsSyUr8oFuxpXB6OSZol3h9fzyewPv8tdWzQS2lBJHL\neej/kJOT+MZuIKfz8DAM0dayXREXT+B9V2qYXm0rwLacGidy1u9Hz3ja6ZdeTuLraN4T6+3ubslX\nK7FSaN1gx3/jCTy0xKdR/7OfLU7krJ9wTP07bzmJ6QZizhu9QVq0nVMz74qnj91YAg/9YHDHG/bu\nOTRZIE7EQG98Fr/xCpveMkFatE3Y50BSCehX29j/HpbLbOVYksjts/9nx7mz+DpdDYfpoghp0a7k\nuX56KgGlHht92raeg41IEbn8iY+FcxKr11kneBrwQoO0aCv5rghLJKDeJrE9bSkiW1+o9c3wazj7\n8wXSov0DNY4n8NybwqlKmoeK7n/DzZ7FQkNIi1YM9z+cFv7QT6IfuLT0X7XDCpaIvQ/f8a6lwyoT\n5CqJP13QdWr7ZZVAdgufXtfr9WtXj1uSjyNi3v1dPrsTKZf4fn6Mg75qo5Vpdz/Gwp/WrdJ9uX4L\ngMjpfXnjUbbfavat9+H66Fjg7aHLz+Ltcz0yqKrAIXFV2fYUdEi8J7Wqcj0krirbnoIOifekVlWu\nh8RVZdtT0CHxntSqylVI/By0KX5IfGzayPDdZ41Gya0q4j6YCImVGidrR5RrI8N3nzUaJRclnurc\nAxMp8ceY89O2MXw3WuPq/wZmFWiTXErMaH8bJpKIastESkynQpQS8+tFZ5R0NlpDVqZNciUEaE4b\nJpKIastESkwvIKIR7DRyrjdaQ1SmUXIB3eXdNkwEEdWYiZR4zLlXyWX2Pet8ozVEZRolt6xrMKMN\nE0HEuxK/V2bkKSQOby7kSCdrfNzwrWgcoPm2Zg3KJYAWleHJIQBtENdut4rJ3B4viKgEk3lgGRMh\ncXBzgSPdW+Pjhm8/Xrai/0KptrvoP4ob17vHiK+BeWwFXhkeiAC0LLDZQUW1YvZ4TkTFmcQCy4gI\nienmclPekQ4jT9zw7f3VZSv658XyNXwudo33y2yPT9vAAQyHjPX1IwBtYWorp1VUi9vjY0S0gTj8\nhgIY8MA1iRqJf/3+ByF4lnvelXekQ2K625lzQ2/oRwuQhbZiDZ9LgM1f/Dw5BKANwhruVjAxq0t7\nPCfib8W2zCEDGVhG5d/f5tsj8LrRs9xkcXSOdEiYMnxjvGzFujUi7nheGZmcD+hmq69kMrPHcyJa\n4qmQ9KXuicEssKzg/EKNU+EyXSqcIx0SpgzfGC9bsW6NiDueV0Ym5wO62errmMzt8ZyIvlC7i6S3\ndTkG88CygnOJp2e59/R1ksmRThImDN80XrRk3RoRd7yojEiOAvrZ6iuZLF2o5TcUiMkssKjeoT3v\nNrzGcRjsvzYYjnQvYdy57scX16xeA7kEKwiJha8fAWiDwCa71Uzm9nhJJMFkHlhGhJ/FLIYc6SkJ\n0Y+WRZcdlK5B8wLYWWWCMYUAtOFYl31aKFUN2x+1x2eJTEyigWVE0hJ7R3rKGo9+tGUrhrNK1/Dz\nguibu2UFPX4XAWj9SKc9v1CqGq7fvAeS9vgcEf25kvtiQCywjEtaYjjSYeyWeOhHK8dLjkvXwLwS\nTDsHAWiLA2snYqFUNdC/2h4P4NWBxCQtMU05dvZdgUPifetXkP0hcUGR9j3lkHjf+hVkf0hcUKR9\nTzkk3rd+BdkfEhcUad9TDon3rV9B9kLiLubvHOh3HO45XFULnAOtxTQ65HCrky0KFBJ3ssrnHOXk\njqn45focrjcQrQTOgf7vki1gKSUOHWfziwA5Wdb9kmwWtBhz/ml9FtcbiH5UsrIKy+WTEtPLdK6v\n7oGTRZEhITpNdmZBizElOfMbmbnfyS4G5ulmQYsxuyUrgZczkhLTi4LzdkfkZFn5K/M50HJMSU75\n8/QvSlZUoaB8UuI25m9R8hwo3DFq0RcuyOk1crj6FJ8caYvAPNscaDlmt2QFcEFGQuLw7vaGN5us\n8NzDTRZ5Gue1oqMsaBwz5nAX5PQv6AVfsUOyPpADox8tJSd3YqDElGOiG2uHUG2ShQIBsgDmGSET\ntDZOSBzciODNXrKx+/EgEbabA+W+cGBFHe6CHLsVI1kfyIHRj5blxw7WJ4u1GUyLZAuAOUsEoHUZ\nCYnprnlT3psNqwr3cJf7qHOgcUz4w6eqRR3lwa3YJ+sDOTD60TI52MH6ZH2hLFC7ZAuAOUsEoHXE\nhMR4PtNWebM5bzYkpuu+ORf0hn60tjPyTw40jhn6wwlwdmJEcH0gB0Y/WgKd7URAiSnHpG6NMTOx\nt0jW5LYALDPyAT6QS0zPZ84qP3mzIaH0cKMf7axcriMLmsLU9ZM/HC+rFsdFoARGv2+j6cZBJ4YS\nk4jPTeyNkl0ClhkpBKDVJLnEeFE4qzy82WAiPdzUn/vtCHPLdI+2UdAUpne6kxKyanFcBEpg9PuW\ngMOdOOjEVGKiAChUANQm2UVgmREC0JqMuMTT8xms8uJCLT3cYIg2YBju5kGFw52w5g53WbU4LhnL\nBTD1z4G/n+zS9VS/rty3WERlKalUsovAMpACohfqmfkb3mwqu/Bwox9tWCrsL4PG7fcRhzuTOInr\nAzkw+tEiwbBNghJDjknd8x+Hb5MsFAiSZMC6n2fkA/zP1fOzOICyF3P3m+eeSTCud9GPlo/GjpZB\nJywyngcoklww5HFjgWYi+tEGwaldD0pMxVSbbDAtGG6QbDWwCExKbL9f7EzdefO3lrr8V+bNm5Ys\nqMOCPzyomVI5Rzlwo4EaBf1oGXDqAKBJhmn3e5Nk/fpBhiVV4IFpieHNhsk7WMbuoh+tHI8eL4AC\nC/7wKEasE7ipQPSjjWHM+gCKrMQEdGOaGE4fIiCVDPoxL40kRhCA1g2nJRbhx+FeK3BIvFflivM+\nJC4u1V4nOont/5my7s/8eyX8V+X9sMLqt+S3i91yDoq/qjA/h+zTKav+AyqklcsSqZ8bAAAAAElF\nTkSuQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}U_{11} & U_{12} & U_{13}\\\\L_{21} U_{11} & L_{21} U_{12} + U_{22} & L_{21} U_{13} + U_{23}\\\\L_{31} U_{11} & L_{31} U_{12} + L_{32} U_{22} & L_{31} U_{13} + L_{32} U_{23} + U_{33}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡  U₁₁           U₁₂                   U₁₃          ⎤\n",
       "⎢                                                   ⎥\n",
       "⎢L₂₁⋅U₁₁    L₂₁⋅U₁₂ + U₂₂         L₂₁⋅U₁₃ + U₂₃     ⎥\n",
       "⎢                                                   ⎥\n",
       "⎣L₃₁⋅U₁₁  L₃₁⋅U₁₂ + L₃₂⋅U₂₂  L₃₁⋅U₁₃ + L₃₂⋅U₂₃ + U₃₃⎦"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = L*U\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Izvedimo sedaj Gaussovo eliminacijo nad matriko koeficientov $\\mathbf{A}$. \n",
    "\n",
    "S pomočjo prve vrstice izvedemo Gaussovo eliminacijo v prvem stolpcu:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAREAAABLCAMAAABgKEOZAAAAPFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAo1xBWAAAAE3RSTlMA\nMquZdlQQQOkwRCLd77tmzYlspWksOAAAAAlwSFlzAAAOxAAADsQBlSsOGwAABcZJREFUeAHtXOu6\nmyAQxGvbxMRY3v9dy21QCCiLm/brOfojqwFmNiOgcVTRSLO04rsvkxVCiEZ2vVqG7y6ImLUMrdSK\nNN9ejFWAeavIvZ3k4zmr0qeUr36tRVjjwCDQCcFNGCgixCJvNp2+Tg/dmAODpAkvYaTIS3UZs3Qj\nKaltZQ6MLd7hOi9hpIh8Ov4TBx4OjEMVthV4CUNFbtINlnHZUpLWOTD+JWGoyIKjzuymE1JqtjIH\nBomWmTBUZPLTiMmpednUEIsSTWMs3fPE1LRLnCHsu7ZqLgwUGYMR2XTtpFNB3E3LF6YxlrsQrdPX\nV+VZSRN26hyiw6RIIgoUGaQ+F1HL3U4ng1FECERbuP+ZxpiUIo1UH/xLmrBTE6H+bfQlUMSPSL1L\n1QIlEEvg0xiTmpfuONcpgSmvkybU7Rn6yAJR3bEXSiCWpJnDUPrKqnF9RJolHB9Vf06CPjK4ft26\n7g0lEI+S0+U5DCGm+tPgPeIcYV/5NyRQRLw6xd200BZKIO4l5ssyGGLR2J9YcoQco0b9J3m23bor\noQRi2e9JY8yfEiSXtJ64ai5whH0k+slQAjEqLtp0bW9q7r+h7xU1rK1kCM0ccsehkwS1r8jDYg0u\nkpBdZdu2eQ3D8PzIzBonZQn14b5/1BDuKHLvJtmqIYQYU5dso+3DXK8raXGyDgibpV9wgKBB7ihC\nA/oytS9F4l15KXIpEisQb1995FIkViDevvrIpUisQLyd7SNjt8xLxWk3t6EUJ/y2zU2YVeSl/iWN\n+lyYvPAaSgX0vIQ5RWbzV2apuTTKaygVKMJLmFPkaS6jVV32Ci5fF/yg01V4CXOKSKPIreIKw9dy\nsPzuGqW5wON/nS84XvFXgk+4YMcsmxrMhJk+cpfG5mxs2NAfr4aGEpwrxOP25BoZQhYHa03mbvtI\nhSKhoQTnCnFlYFtLEzI5WGuW9aMmNJTgXCGuDGxracJzDtaPn7/e8rMz60CfWf2oNi4YnCvEN57z\nX6QJNW6dg/X7Z/o+tKfxTGe66RQbSqtzVXUoP1YsS8jhYG3o7RlaRz9Diw2l1bn6Sw4WCHkcrI0k\nkz6Lr/EJI0PJO1d/y8HyhJWjJnP0VWos6tbOin96saEE5wpxIzrXamiZrUQfcLAYUoZzhcgAuQ9h\niT7lYO1zl5TCuUIsaXOqDog+4mCdysw1hnOFyIG5iwGiy8HalYlUmJ1ZSShfqfKlSLw3L0UuRWIF\n4u2rj1yKxArE21cfKVaE4mDNrXpkq03fBldrMH0CU90vdfyUWbaP0Bws6e4XjwXX27UG0ycwC7LJ\nKUJzsG57V6grDSYezGG9F9XurcNscorQHCzfDVJ9pNJg4sF8U+Qwm5wiNAfLC59QxHs+xOe6eDBj\nRY6zyShCvBYvd64++ivDREeLBzNW5DibjCI0B2s75BfnHHnLKmMwHT2TVYOZMK1iRaJskOymZVYR\n43IWOlibIQ/nyFtWaYPJFyeGmf2KjgnqADJSJMwGLRBNy4witFHjh/wo4Bx5yyptMPniIP3tBh0T\n1A5leellepiAh3XDbNACcU8RQXKwMH/P9gYc7Rx5y8oPXKKjRcfUv+fNtIr6SJjNtoVvmekjguJg\n+fnbPjDiPQ1jWWUNpv1nsuowPbXZ2/rjTRFXgk6DFoj26b7UuyUoDhaE781zj6tzZCyrSkerCnOl\nzikSZ4MWiKpdro8IgoPl5u8FlqjrgM6yqnO06jAPR038lNk6zg5HTbGDNbYvKdvWfLo9Y50j7ySl\nDSZf7HfnulKLmXjsKho1sb+2tvB2V7aPrPlR1sxwNM8+ZSwrfI1YAF6IuVbbYr4psilEC0RbxKyI\ngHMEJ2mTgF7F14hRcXqzEBPVApBx7/5TtEA0LbkVgXMEJylITwh8jRgVpzcLMVEtDZL6Fi0QTR1u\nRVLE/9d3lyLx/roUSStyvSHO6OLfEDfqF6P1ff1raWKZ/9dt84Y49ZavP2qxVg8F9nyTAAAAAElF\nTkSuQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}U_{11} & U_{12} & U_{13}\\\\0 & U_{22} & U_{23}\\\\0 & L_{32} U_{22} & L_{32} U_{23} + U_{33}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡U₁₁    U₁₂         U₁₃     ⎤\n",
       "⎢                           ⎥\n",
       "⎢ 0     U₂₂         U₂₃     ⎥\n",
       "⎢                           ⎥\n",
       "⎣ 0   L₃₂⋅U₂₂  L₃₂⋅U₂₃ + U₃₃⎦"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[1,:] -= L21 * A[0,:]\n",
    "A[2,:] -= L31 * A[0,:]\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Nadaljujemo v drugem stolpcu:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJkAAABLCAMAAABUdNWsAAAAPFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAo1xBWAAAAE3RSTlMA\nMquZdlQQQOkwRCLd77tmzYlspWksOAAAAAlwSFlzAAAOxAAADsQBlSsOGwAABDtJREFUaAXtmumW\noyAQhXHtnmg0Du//rsN22aS0yOlJ26fjjxQKF66A6GcUjTRbK66yzdaQEI3serUNVzEmVm2nldpZ\ncxlTwcgaO5vaWd6XVeUuUj76UOoo9YyIo0mcCbHJm3XRM33p0s+IzjWZs4fqQrN1o0swwjOic03m\nTC7OSc2F+ozoXJM6u0k3iOPG6CtX5BkRQ5M623CVrm66cfw9I2JoUmezn2bGUvOwzhDLPsuirVsO\n5iqh6bvWz+/E2ZiMftO1s/aCWPYlyqJtEqJ1J7YXljWdWq86THSROBukXsvUNtnpNhhnQiDazOy3\nLJqVs0aqn+JW1nRqcms/dkuc+dHXZ6w2OEK0kuy3LJrVRJ2wOGYKvQK6u07WENlnGwy7NQOOEHct\n6AOUSJ2Y9JMmU5Ka8e5vlEmfDa77WzcKcISY1W93KZEQM3kfoTR9dEtMnIlHp1prWviGI8SiM0ok\nNl0ZsRENkaOphmZpu3CmcIRINFMWrQfGqIb03MTDWNpnWdNwhJhll3dd4Zu60G7o/HLJcNRozByb\nsDykq0YoalPDPY15fnHfiprHMAwLdQXkQqvRK01/h+agz6Zulq0aWsS8uuI+Ct/NA3OxyO4gNM3W\nb7j47Mr2A55pd2fznQcORvM7bam2387qB+DdZ7+iz8ZuWzfO3YVDtXmPcTTkPHuoO+uo7xfn2znV\n7us411DOVnPL3Mgn+bitc6qNS9v0uYZytpjHWvqpNG4r4Zo44yB9rqGcSePs5p+WDhphUO1OzdAQ\nzkZpnvt8BbuqowMeNyrwmaEhnE3SvD5obIhsFJIp1YKAEQsCdYjQkCQcaplsn3GcpVQLAkYMVcap\nsuaAhIO4YjRTqgUBI4Yq41RZsyfhj88/scyk7RUwMK4AP2UM1YKAEXcVmwNljc4Kbw/+fpbf0y7m\n/cJKsmxoMafaQMD0mkNqKBIOzQm70naMlTan2kDAfBKGhibhyNqs706B5aOcPJlRrSfgChL2mmg0\niVVDudrUK3nOHT2nWhAwYn4idj+l51CWScLlSo+PgoARj0vbXFu2hoQ5tWZlQMCIWXZxF2XZJFys\n5ewgCBjxrLzOR9k3CXN666gMeW0eiV6S93ZW383vPvsVfcYm4bg3uKJXkXDwxsfn15Cwd1aBz68h\nYe+sAp9fQsLemODjswdZ+s9nYj2rYKdgrELkEYWmZ8JZDQl7axWijIQdAW9fS8LemPrPyLxyYOBz\nSsIgYERbI9FnFQMTnPFFKQmDgBEPnbnJzCHh4Iwv8tPM/SccCDiQMNFngk/CkTO2aEfC4EdEVSnl\njE/CkTO2KKdnEDCirpNyJvgkHFlji3J69gR8Ppo1JBys8fE5JeHwX/B/JOFgk5sCASM6HTma3Hq/\noBwIGNFWeQVnIGDE6zgrd/sV+uxnOrvql76j/rC27ys+0yv3/5cdNV/6qq9T/wEc4T3rbuDb7AAA\nAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}U_{11} & U_{12} & U_{13}\\\\0 & U_{22} & U_{23}\\\\0 & 0 & U_{33}\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡U₁₁  U₁₂  U₁₃⎤\n",
       "⎢             ⎥\n",
       "⎢ 0   U₂₂  U₂₃⎥\n",
       "⎢             ⎥\n",
       "⎣ 0    0   U₃₃⎦"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[2,:] -= L32 * A[1,:]\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Iz zgornje eliminacije ugotovimo:\n",
    "1. matrika $\\mathbf{U}$ je enaka matriki, ki jo dobimo, če izvedemo Gaussovo eliminacijo nad matriko koeficientov $\\mathbf{A}$.\n",
    "* izven diagonalni členi $\\mathbf{L}$ so faktorji, ki smo jih uporabili pri Gaussovi eliminaciji."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Numerična implementacija razcepa LU\n",
    "\n",
    "Numerično implementacijo si bomo pogledali na sistemu, ki je definiran kot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "A = np.array([[8, -6, 3],\n",
    "              [-6, 6, -6],\n",
    "              [3, -6, 6]], dtype=float) \n",
    "b = np.array([-14, 36, 6], dtype=float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Izvedimo Gaussovo eliminacijo in koeficiente, s katerim množimo pivotno vrsto `m`, shranjuje v matriko `L` na mesto z indeksi, kot jih ima v matriki $\\mathbf{A}$ eliminirani element."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Korak: 0\n",
      "[[ 8.    -6.     3.   ]\n",
      " [ 0.     1.5   -3.75 ]\n",
      " [ 0.    -3.75   4.875]]\n",
      "Korak: 1\n",
      "[[ 8.   -6.    3.  ]\n",
      " [ 0.    1.5  -3.75]\n",
      " [ 0.    0.   -4.5 ]]\n"
     ]
    }
   ],
   "source": [
    "(v, s) = A.shape # v=število vrstic, s=število stolpcev\n",
    "U = A.copy()     # pripravimo matriko U kot kopijo A\n",
    "L = np.zeros_like(A) # pripravimo matriko L dimenzije enake A (vrednosti 0)\n",
    "# eliminacija\n",
    "for p, pivot_vrsta in enumerate(U[:-1]):\n",
    "    for i, vrsta in enumerate(U[p+1:]):\n",
    "        if pivot_vrsta[p]:\n",
    "            m = vrsta[p]/pivot_vrsta[p]\n",
    "            vrsta[p:] = vrsta[p:]-pivot_vrsta[p:]*m\n",
    "            L[p+1+i, p] = m\n",
    "    print('Korak: {:g}'.format(p))\n",
    "    print(U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.   ,  0.   ,  0.   ],\n",
       "       [-0.75 ,  0.   ,  0.   ],\n",
       "       [ 0.375, -2.5  ,  0.   ]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dopolnimo diagonalo $\\mathbf{L}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(v):\n",
    "    L[i, i] = 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.   ,  0.   ,  0.   ],\n",
       "       [-0.75 ,  1.   ,  0.   ],\n",
       "       [ 0.375, -2.5  ,  1.   ]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Sedaj rešimo spodnje trikotni sistem enačb $\\mathbf{L}\\,\\mathbf{y}=\\mathbf{b}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# direktno vstavljanje\n",
    "y = np.zeros_like(b)\n",
    "for i, b_ in enumerate(b):\n",
    "    y[i] = (b_ - np.dot(L[i, :i], y[:i]))    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-14. ,  25.5,  75. ])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Nadaljujemo z reševanjem zgornje trikotnega sistema $\\mathbf{U}\\,\\mathbf{x}=\\mathbf{y}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 8.  , -6.  ,  3.  ],\n",
       "       [ 0.  ,  1.5 , -3.75],\n",
       "       [ 0.  ,  0.  , -4.5 ]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-14. ,  25.5,  75. ])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-14.        , -24.66666667, -16.66666667])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# obratno vstavljanje\n",
    "x = np.zeros_like(b)\n",
    "for i in range(v-1, -1,-1):\n",
    "    x[i] = (y[i] - np.dot(U[i, i+1:], x[i+1:])) / U[i, i]\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Kakor smo navedli zgoraj, ob spremembi vektorja konstant $\\mathbf{b}$ ponovna Gaussova eliminacija ni potrebna. Izvesti je treba samo direktno in nato obratno vstavljanje. Poglejmo primer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-4.33333333, -7.88888889, -4.55555556])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = np.array([-1., 6., 7.])\n",
    "y = np.zeros_like(b)\n",
    "for i, b_ in enumerate(b):#direktno vstavljanje\n",
    "    y[i] = (b_ - np.dot(L[i, :i], y[:i]))    \n",
    "x = np.zeros_like(b)\n",
    "for i in range(v-1, -1,-1):# obratno vstavljanje\n",
    "    x[i] = (y[i] - np.dot(U[i, i+1:], x[i+1:])) / U[i, i]\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.,  6.,  7.])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.dot(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Pivotiranje"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Poglejmo si spodnji sistem enačb:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[0, -6, 6],\n",
    "              [-6, 6, -6],\n",
    "              [8, -6, 3]], dtype=float) # poskusite tukaj izpustiti dtype=float in preveriet rezultat\n",
    "b = np.array([6, 36, -14], dtype=float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Če bi izvedli Gaussovo eliminacijo v prvem stolpcu matrike $\\mathbf{A}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Janko\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:1: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n",
      "C:\\Users\\Janko\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in multiply\n",
      "  \"\"\"Entry point for launching an IPython kernel.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ nan, -inf,  inf])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[1,:] - A[1,0]/A[0,0] * A[0,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Opazimo, da imamo težavo z deljenjem z 0 v prvi vrstici. Elementarne operacije, ki jih nad sistemom lahko izvajamo, dovoljujejo zamenjavo poljubnih vrstic. Sistem lahko preuredimo tako, da pivotni element ni enak 0. Vseeno se lahko zgodi, da ima pivotni element, s katerim delimo, zelo majhno vrednost $\\varepsilon$. Ker bi to povečevalo zaokrožitveno napako, izmed vseh vrstic za pivotno vrstico izberemo tisto, katere pivot ima največjo absolutno vrednost.\n",
    "\n",
    "Če med Gaussovo eliminacijo zamenjamo vrstice tako, da je pivotni element največji, to imenujemo **pivotiranje vrstic** ali tudi **delno pivotiranje**. Tako dosežemo, da je Gaussova eliminacija numerično stabilna."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Pokazati je mogoče, da pri reševanju sistema enačb $\\mathbf{A}\\,\\mathbf{x}=\\mathbf{b}$, pri katerem je matrika $\\mathbf{A}$ diagonalno dominantna, pivotiranje po vrsticah ni potrebno. Reševanje je brez pivotiranja numerično stabilno.\n",
    "\n",
    "Pravokotna matrika $\\mathbf{A}$ dimenzije $n$ je **diagonalno dominantna**, če je absolutna vrednost diagonalnega elementa vsake vrstice večja od vsote absolutnih vrednosti ostalih elementov v vrstici:\n",
    "$$\\left|A_{ii}\\right|> \\sum_{j=1, j\\ne i}^n \\left|A_{ij}\\right|$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Gaussova eliminacija z delnim pivotiranjem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pogledali si bomo Gaussovo eliminacijo z **delnim pivotiranjem**. Brez delnega pivotiranja v $i$-tem koraku eliminacije izberemo vrstico $i$ za pivotiranje. Pri delnem pivotiranju pa najprej preverimo, ali je $i$-ti diagonalni element po absolutni vrednosti največji element v stolpcu $i$ na ali pod diagonalo; če ni, zamenjamo vrstico $i$ s \n",
    "tisto vrstico pod njo, v kateri je v stolpcu $i$ po absolutni vrednosti največji element. Z delnim pivotiranjem zmanjšamo vpliv zaokrožitvene napake na rezultat. \n",
    "\n",
    "V koliko bi izvedli **polno pivotiranje**, bi poleg zamenjave vrstic uporabili tudi zamenjavo vrstnega reda spremenljivk (zamenjava stolpcev). Polno pivotiranje izboljša stabilnost, se pa redko uporablja in ga tukaj ne bomo obravnavali."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Algoritem za Gaussovo eliminacijo z delnim pivotiranjem torej je:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def gaussova_eliminacija_pivotiranje(A, b, prikazi_korake=False):\n",
    "    \"\"\" Vrne Gaussovo eliminacijo razširjene matrike koeficientov, \n",
    "    uporabi delno pivotiranje.\n",
    "\n",
    "    :param A: matrika koeficientov\n",
    "    :param b: vektor konstant\n",
    "    :param prikazi_korake: ali izpišem posamezne korake\n",
    "    :return Ab: trapezna razširjena matrika koeficientov\n",
    "    \"\"\"\n",
    "    Ab = np.column_stack((A, b))\n",
    "    for p in range(len(Ab)-1):\n",
    "        p_max = np.argmax(np.abs(Ab[p:,p]))+p\n",
    "        if p != p_max:\n",
    "            Ab[[p], :], Ab[[p_max], :] = Ab[[p_max], :], Ab[[p], :]\n",
    "        pivot_vrsta = Ab[p, :]                \n",
    "        for vrsta in Ab[p + 1:]:\n",
    "            if pivot_vrsta[p]:\n",
    "                vrsta[p:] -= pivot_vrsta[p:] * vrsta[p] / pivot_vrsta[p]\n",
    "        if prikazi_korake:\n",
    "            print('Korak: {:g}'.format(p))\n",
    "            print('Pivot vrsta:', pivot_vrsta)\n",
    "            print(Ab)\n",
    "    return Ab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0., -6.,  6.],\n",
       "       [-6.,  6., -6.],\n",
       "       [ 8., -6.,  3.]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Korak: 0\n",
      "Pivot vrsta: [  8.  -6.   3. -14.]\n",
      "[[  8.    -6.     3.   -14.  ]\n",
      " [  0.     1.5   -3.75  25.5 ]\n",
      " [  0.    -6.     6.     6.  ]]\n",
      "Korak: 1\n",
      "Pivot vrsta: [ 0. -6.  6.  6.]\n",
      "[[  8.    -6.     3.   -14.  ]\n",
      " [  0.    -6.     6.     6.  ]\n",
      " [  0.     0.    -2.25  27.  ]]\n"
     ]
    }
   ],
   "source": [
    "Ab = gaussova_eliminacija_pivotiranje(A, b, prikazi_korake=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Razcep LU z delnim pivotiranjem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Podobno kakor pri Gaussovi eliminaciji lahko tudi razcep LU razširimo z delnim pivotiranjem. Reševanje tako postane numerično **stabilno**. Pri tem moramo shraniti informacijo o zamenjavi vrstic, ki jo potem posredujemo v funkcijo za rešitev ustreznih trikotnih sistemov."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def LU_razcep_pivotiranje(A, prikazi_korake=False):\n",
    "    \"\"\" Vrne razcep LU matriko in vektor zamenjanih vrstic pivotiranje, \n",
    "    uporabi delno pivotiranje.\n",
    "\n",
    "    :param A: matrika koeficientov\n",
    "    :param prikazi_korake: izpišem posamezne korake\n",
    "    :return LU: LU matrika\n",
    "    :return pivotiranje: vektor zamenjave vrstic (pomembno pri iskanju rešitve)    \n",
    "    \"\"\"\n",
    "    LU = A.copy()\n",
    "    pivotiranje = np.arange(len(A))\n",
    "    for p in range(len(LU)-1):\n",
    "        p_max = np.argmax(np.abs(LU[p:,p]))+p\n",
    "        if p != p_max:\n",
    "            LU[[p], :], LU[[p_max], :] = LU[[p_max], :], LU[[p], :]\n",
    "            pivotiranje[p], pivotiranje[p_max] = pivotiranje[p_max], pivotiranje[p]\n",
    "        pivot_vrsta = LU[p, :]                \n",
    "        for vrsta in LU[p + 1:]:\n",
    "            if pivot_vrsta[p]:\n",
    "                m = vrsta[p] / pivot_vrsta[p]\n",
    "                vrsta[p:] -= pivot_vrsta[p:] * m\n",
    "                vrsta[p] = m\n",
    "            else:\n",
    "                raise Exception('Deljenje z 0.')\n",
    "        if prikazi_korake:\n",
    "            print('Korak: {:g}'.format(p))\n",
    "            print('Pivot vrsta:', pivot_vrsta)\n",
    "            print(LU)\n",
    "    return LU, pivotiranje"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Poglejmo si primer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "A = np.array([[0, -6, 6],\n",
    "              [-6, 6, -6],\n",
    "              [8, -6, 3]], dtype=float)\n",
    "b = np.array([-14, 36, 6], dtype=float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Korak: 0\n",
      "Pivot vrsta: [ 8. -6.  3.]\n",
      "[[ 8.   -6.    3.  ]\n",
      " [-0.75  1.5  -3.75]\n",
      " [ 0.   -6.    6.  ]]\n",
      "Korak: 1\n",
      "Pivot vrsta: [ 0. -6.  6.]\n",
      "[[ 8.   -6.    3.  ]\n",
      " [ 0.   -6.    6.  ]\n",
      " [-0.75 -0.25 -2.25]]\n"
     ]
    }
   ],
   "source": [
    "lu, piv = LU_razcep_pivotiranje(A, prikazi_korake=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "V zgornjem primeru smo uporabili kompaketen način zapisa trikotnih matrik $\\mathbf{L}$ in $\\mathbf{U}$; vsaka je namreč definirana s 6 elementi, pri matriki $\\mathbf{L}$ pa vemo, da so diagonalni elementi enaki 1. Matrika `lu` tako vsebuje $3\\times3=9$ elementov:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 8.  , -6.  ,  3.  ],\n",
       "       [ 0.  , -6.  ,  6.  ],\n",
       "       [-0.75, -0.25, -2.25]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Na diagonali in nad diagonalo so vrednosti zgornje trikotne matrike $\\mathbf{U}$, pod diagonalo pa so poddiagonalni elementi matrike $\\mathbf{L}$. Pri izračunu rešitve bomo upoštevali, da so diagonalne vrednosti $\\mathbf{L}$ enake 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numerični seznam `piv` nam pove, kako so bile zamenjane vrstice, kar je treba upoštevati pri izračunu rešitve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 0, 1])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "piv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Določimo sedaj še rešitev (koda za računanje rešitve je v modulu `orodja.py`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "from moduli import orodja"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -3.66666667, -14.11111111, -16.44444444])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = orodja.LU_resitev_pivotiranje(lu, b, piv)\n",
    "r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Preverjanje rešitve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-14.,  36.,   6.])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A@r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Modul ``SciPy``"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Modul ``SciPy`` temelji na ``numpy`` modulu in vsebuje veliko različnih visokonivojskih programov/modulov/funkcij. Teoretično ozadje modulov so seveda različni numerični algoritmi; nekatere spoznamo tudi v okviru tega učbenika. Dober vir teh numeričnih algoritmov v povezavi s `SciPy` predstavlja [dokumentacija](https://docs.scipy.org/doc/scipy/reference/tutorial/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Kratek pregled hierarhije modula:\n",
    "\n",
    "* Linearna algebra ([scipy.linalg](http://docs.scipy.org/doc/scipy/reference/linalg.html))\n",
    "* Integracija ([scipy.integrate](http://docs.scipy.org/doc/scipy/reference/integrate.html))\n",
    "* Optimizacija ([scipy.optimize](http://docs.scipy.org/doc/scipy/reference/optimize.html))\n",
    "* Interpolacija ([scipy.interpolate](http://docs.scipy.org/doc/scipy/reference/interpolate.html))\n",
    "* Fourierjeva transformacija ([scipy.fftpack](http://docs.scipy.org/doc/scipy/reference/fftpack.html))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "* Problem lastnih vrednosti (redke matrike) ([scipy.sparse](http://docs.scipy.org/doc/scipy/reference/sparse.html))\n",
    "* Statistika ([scipy.stats](http://docs.scipy.org/doc/scipy/reference/stats.html))\n",
    "* Procesiranje signalov ([scipy.signal](http://docs.scipy.org/doc/scipy/reference/signal.html))\n",
    "* Posebne funkcije ([scipy.special](http://docs.scipy.org/doc/scipy/reference/special.html))\n",
    "* Večdimenzijsko procesiranje slik ([scipy.ndimage](http://docs.scipy.org/doc/scipy/reference/ndimage.html))\n",
    "* Delo z datotekami ([scipy.io](http://docs.scipy.org/doc/scipy/reference/io.html))\n",
    "\n",
    "V sledečih predavanjih si bomo nekatere podmodule podrobneje pogledali."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Poglejmo si, kako je znotraj ``SciPy`` implementiran *razcep LU*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import lu_factor, lu_solve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Funkcija ``scipy.linalg.lu_factor`` ([dokumentacija](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu_factor.html)):\n",
    "\n",
    "```python\n",
    "lu_factor(a, overwrite_a=False, check_finite=True)\n",
    "```\n",
    "\n",
    "zahteva vnos matrike koeficientov (ali seznama matrik koeficientov) `a`, `overwrite_a` v primeru `True` z rezultatom razcepa prepiše vrednost `a` (to je lahko pomembno, da se prihrani spomin in poveča hitrost). Funkcija `lu_factor` vrne terko `(lu, piv)`:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "* `lu` - $\\mathbf{L}$ \\ $\\mathbf{U}$ matrika (matrika, ki je enake dimenzije kot `a`, vendar pod diagonalo vsebuje elemente $\\mathbf{L}$, preostali elementi pa definirajo $\\mathbf{U}$; diagonalni elementi $\\mathbf{L}$ imajo vrednosti 1.\n",
    "* `piv` - pivotni indeksi, predstavljajo permutacijsko matriko `P`: vrsta  `i` matrike `a` je bila zamenjana z vrsto `piv[i]` (v vsakem koraku se upošteva predhodno stanje, malo drugačna logika kot v naši funkciji `LU_razcep_pivotiranje`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "``lu_factor`` uporabljamo v paru s funkcijo ``scipy.linalg.lu_solve`` ([dokumentacija](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lu_solve.html)), ki nam poda rešitev sistema:\n",
    "\n",
    "```python\n",
    "lu_solve(lu_and_piv, b, trans=0, overwrite_b=False, check_finite=True)\n",
    "```\n",
    "\n",
    "`lu_and_piv` je terka rezultata `(lu, piv)` iz ``lu_factor``, `b` je vektor (ali seznam vektorjev) konstant. Ostali parametri so opcijski.\n",
    "\n",
    "Poglejmo si uporabo:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "A = np.array([[0, -6, 6],            \n",
    "              [-6, 6, -6],\n",
    "              [8, -6, 3]], dtype=float)\n",
    "b = np.array([-14, 36, 6], dtype=float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "lu, piv = lu_factor(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 8.  , -6.  ,  3.  ],\n",
       "       [ 0.  , -6.  ,  6.  ],\n",
       "       [-0.75, -0.25, -2.25]])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 2, 2], dtype=int32)"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "piv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pridobimo rešitev:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -3.66666667, -14.11111111, -16.44444444])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = lu_solve((lu, piv), b)\n",
    "r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Preverimo ustreznost rešitve:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-14.,  36.,   6.])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A@r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Računanje inverzne matrike"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inverzno matriko h kvadratni matriki $\\textbf{A}$ reda $n\\times n$ označimo z $\\textbf{A}^{-1}$. Je matrika reda $n\\times n$,  takšna, da velja \n",
    "$$\\mathbf{A}\\,\\mathbf{A}^{-1} = \\mathbf{A}^{-1}\\,\\mathbf{A} = \\mathbf{I},$$\n",
    "kjer je $\\mathbf{I}$ enotska matrika.\n",
    "\n",
    "Najbolj učinkovit način za izračun inverzne matrike od matrike $\\mathbf{A}$ je rešitev matrične enačbe:\n",
    "$$\\mathbf{A}\\,\\mathbf{X}=\\mathbf{I}.$$ \n",
    "Matrika $\\mathbf{X}$ je inverzna matriki $\\mathbf{A}$: $\\mathbf{A}^{-1}=\\mathbf{X}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Izračun inverzne matrike je torej enak reševanju $n$ sistemov $n$ linarnih enačb:\n",
    "$$\\mathbf{A}\\,\\mathbf{x}_i=\\mathbf{b}_i,\\qquad i=0,1,2,\\dots,n-1,$$\n",
    "kjer je $\\mathbf{b}_i$ $i$-ti stolpec matrike $\\mathbf{B}=\\mathbf{I}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "*Numerična zahtevnost:* izvedemo razcep LU nad matriko $\\mathbf{A}$ (računski obseg reda $n^3$) in nato \n",
    "poiščemo rešitev za vsak $\\mathbf{x}_i$ ($2\\,n^2$ računskih operacij za vsak $i$). Skupni računski obseg je torej $3\\,n^3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[ 8.  , -6.  ,  3.  ],\n",
       "        [ 0.  , -6.  ,  6.  ],\n",
       "        [-0.75, -0.25, -2.25]]), array([2, 2, 2], dtype=int32))"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lu_piv = lu_factor(A)\n",
    "lu_piv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tukaj se sedaj pokaže smisel razcepa LU; razcep namreč izračunamo samo enkrat in potem za vsak vektor konstant poiščemo rešitev tako, da rešimo dva trikotna sistema. Če bi sisteme reševali po Gaussovi metodi, bi rabili $(n^3+n^2)\\,n$ računskih operacij. \n",
    "\n",
    "Izračunajmo inverzno matriko od matrike $\\textbf{A}$. Najprej z uporabo `np.identity()` pripravimo enotsko matriko:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.,  0.],\n",
       "       [ 0.,  1.,  0.],\n",
       "       [ 0.,  0.,  1.]])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.identity(len(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nato rešimo sistem enačb za vsak stolpec enotske matrike:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_inv = np.zeros_like(A) #zeros_like vrne matriko oblike in tipa matrike A z vrednostmi 0\n",
    "for b, a_inv in zip(np.identity(len(A)), A_inv):\n",
    "    a_inv[:] = lu_solve(lu_piv, b)\n",
    "A_inv = A_inv.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Rešitev torej je:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ -1.66666667e-01,  -1.66666667e-01,   1.38777878e-17],\n",
       "       [ -2.77777778e-01,  -4.44444444e-01,  -3.33333333e-01],\n",
       "       [ -1.11111111e-01,  -4.44444444e-01,  -3.33333333e-01]])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_inv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Preverimo rešitev:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.00000000e+00,   0.00000000e+00,   1.11022302e-16],\n",
       "       [ -1.11022302e-16,   1.00000000e+00,  -1.11022302e-16],\n",
       "       [  1.66533454e-16,   0.00000000e+00,   1.00000000e+00]])"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A@A_inv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Rešitev z uporabo Numpy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ -1.66666667e-01,  -1.66666667e-01,   1.38777878e-17],\n",
       "       [ -2.77777778e-01,  -4.44444444e-01,  -3.33333333e-01],\n",
       "       [ -1.11111111e-01,  -4.44444444e-01,  -3.33333333e-01]])"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#%%timeit\n",
    "np.linalg.inv(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Z ustrezno uporabo funkcij iz modula `SciPy` lahko do rešitve pridemo še hitreje. V funkcijo ``lu_solve`` lahko vstavimo vektor $\\mathbf{b}$ ali matriko $\\mathbf{B}$, katere posamezni stolpec $i$ predstavlja nov vektor konstant $\\mathbf{b}_i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ -1.66666667e-01,  -1.66666667e-01,   1.38777878e-17],\n",
       "       [ -2.77777778e-01,  -4.44444444e-01,  -3.33333333e-01],\n",
       "       [ -1.11111111e-01,  -4.44444444e-01,  -3.33333333e-01]])"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#%%timeit\n",
    "lu_solve(lu_piv, np.identity(len(A)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Reševanje predoločenih sistemov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kadar rešujemo sistem $m$ linearnih enačbami z $n$ neznankami ter velja $m>n$ in je rang $n+1$, imamo predoločeni sistem.\n",
    "\n",
    "Predoločeni (tudi **nekonsistenten sistem**):\n",
    "$$\\mathbf{A}\\,\\mathbf{x}=\\mathbf{b},$$\n",
    "nima rešitve. Lahko pa poiščemo najboljši približek rešitve z metodo **najmanjših kvadratov** . "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Vsota kvadratov preostankov je definirana s skalarnim produktom:\n",
    "$$\\left\\lVert r\\right\\rVert ^2=(\\mathbf{A}\\,\\mathbf{x}-\\mathbf{b})^T\\,(\\mathbf{A}\\,\\mathbf{x}-\\mathbf{b})$$\n",
    "kar preoblikujemo v:\n",
    "$$\\left\\lVert r\\right\\rVert^2=\\mathbf{x}^T\\mathbf{A}^T\\mathbf{A}\\mathbf{x}-2\\mathbf{b}^T\\mathbf{A}\\mathbf{x}+\\mathbf{b}^T\\mathbf{b},$$\n",
    "kjer smo upoštevali, da zaradi skalarne vrednosti velja: $\\mathbf{b}^T\\mathbf{A}\\mathbf{x}=(\\mathbf{b}^T(\\mathbf{A}\\mathbf{x}))^T=(\\mathbf{A}\\mathbf{x})^T\\mathbf{b}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Rešitev enačbe, gradient vsote kvadratov, določa njen minimum:\n",
    "$$\\nabla_x\\,\\left\\lVert r\\right\\rVert^2=2\\mathbf{A}^T\\mathbf{A}\\,\\mathbf{x}-2\\mathbf{A}^T\\,\\mathbf{b}=0$$\n",
    "Tako iz normalne enačbe:\n",
    "$$\\mathbf{A}^T\\mathbf{A}\\,\\mathbf{x}=\\mathbf{A}^T\\,\\mathbf{b}$$\n",
    "določimo najboljši približek rešitve:\n",
    "$$\\mathbf{x}=\\left(\\mathbf{A}^T\\mathbf{A}\\right)^{-1}\\mathbf{A}^T\\,\\mathbf{b}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Z vpeljavo **psevdo inverzne matrike**: \n",
    "$$\\mathbf{A}^+=\\left(\\mathbf{A}^T\\mathbf{A}\\right)^{-1}\\mathbf{A}^T$$\n",
    "k matriki $\\mathbf{A}$ je rešitev predoločenega sistema zapisana:\n",
    "$$\\mathbf{x}=\\mathbf{A}^+\\,\\mathbf{b}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Psevdo inverzno matriko lahko izračunamo z uporabo funkcij:\n",
    "\n",
    "* `numpy.linalg.pinv` iz modula `numpy` ([dokumentacija](https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.pinv.html)),\n",
    "* Funkcij `pinv`, `pinv2` ali `pinvh` iz modula `scipy.linalg` (izbira je odvisna od obravnavanega problema; glejte [dokumentacijo](https://docs.scipy.org/doc/scipy/reference/linalg.html))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Primer sistema z enolično rešitvijo:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.,  2.])"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# število enačb enako številu neznank\n",
    "A = np.array([[1., 2],\n",
    "              [2, 3]])\n",
    "b = np.array([5., 8])\n",
    "np.linalg.solve(A, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Naredimo sedaj predoločeni sistem (funkcija `numpy.vstack` ([dokumentacija](https://docs.scipy.org/doc/numpy/reference/generated/numpy.vstack.html)) sestavi sezname po stolpcih, `numpy.random.seed` ([dokumentacija](https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.seed.html)) ponastavi generator naključnih števil na vrednost semena `seed`, `numpy.random.normal` ([dokumentacija](https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.normal.html)) pa generira normalno porazdeljeni seznam dolžine `size` in standardne deviacije `scale`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.01764052,  2.00400157],\n",
       "       [ 2.00978738,  3.02240893],\n",
       "       [ 1.01867558,  1.99022722],\n",
       "       [ 2.00950088,  2.99848643],\n",
       "       [ 0.99896781,  2.00410599],\n",
       "       [ 2.00144044,  3.01454274]])"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vA = np.vstack([A,A,A]) \n",
    "np.random.seed(seed=0)\n",
    "vA += np.random.normal(scale=0.01, size=vA.shape) # pokvarimo rešitev -> sistem je predoločen\n",
    "vA # matrika koeficientov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 5.00761038,  8.00121675,  5.00443863,  8.00333674,  5.01494079,\n",
       "        7.99794842])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vb = np.hstack([b,b,b])\n",
    "vb += np.random.normal(scale=0.01, size=vb.shape)\n",
    "vb # vektor konstant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Rešimo sedaj predoločen sistem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.94512271,  2.02675939])"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ap = np.linalg.pinv(vA)\n",
    "Ap.dot(vb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vidimo, da predoločeni sistem z naključnimi vrednostmi (simulacija šuma pri meritvi) poda podoben rezultat kakor rešitev brez šuma. V kolikor bi nivo šuma povečevali, bi se odstopanje od enolične rešitve povečevalo."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Psevdo inverzno matriko lahko določimo tudi sami in preverimo razliko z vgrajeno funkcijo:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  8.88178420e-15,  -8.21565038e-15,   9.65894031e-15,\n",
       "         -8.65973959e-15,   1.04360964e-14,  -7.10542736e-15],\n",
       "       [ -5.66213743e-15,   4.71844785e-15,  -5.66213743e-15,\n",
       "          4.49640325e-15,  -6.21724894e-15,   3.55271368e-15]])"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#%%timeit # preverite hitrost!\n",
    "Ap2 = np.linalg.inv(vA.T@vA) @ vA.T\n",
    "\n",
    "Ap2 - Ap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Iterativne metode"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pogosto se srečamo z velikimi sistemi linearnih enačb, katerih matrika koeficientov ima malo od nič različnih elementov (take matrike imenujemo **redke** ali tudi **razpršene**, angl. *sparse*). \n",
    "\n",
    "Pri reševanju takih sistemov linearnih enačb se zelo dobro izkažejo iterativne metode; prednosti v primerjavi z direktnimi metodami so:\n",
    "\n",
    "* računske operacije se izvajajo samo nad neničelnimi elementi (kljub iterativnemu reševanju jih je lahko manj)\n",
    "* zahtevani spominski prostor je lahko neprimerno manjši."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Gauss-Seidlova metoda"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "V nadaljevanju si bomo pogledali idejo *Gauss-Seidelove* iterativne metode. Najprej sistem enačb $\\mathbf{A}\\,\\mathbf{x}=\\mathbf{b}$ zapišemo kot:\n",
    "$$\\sum_{j=0}^{n-1} A_{ij}x_j=b_i\\qquad{}i=0,1,\\dots,n-1.$$\n",
    "\n",
    "Predpostavim, da smo v $k-1$ koraku iterativne metode in so znani približki $x_{j}^{(k-1)}$ ($j=0,1,\\dots, n-1$). Iz zgornje vsote izpostavimo člen $i$:\n",
    "$$A_{ii}\\,x_{i}^{(k-1)} +\\sum_{j=0, j\\ne i}^{n-1} A_{ij}x_{j}^{(k-1)}=b_i\\qquad{}i=0,1,\\dots,n-1,$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Ker približki $x_{j}^{(k-1)}$ ne izpolnjujejo natančno linearnega problema, lahko iz zgornje enačbe določimo nov približek $x_{i}^{(k)}$:\n",
    "$$x_{i}^{(k)} =\\frac{1}{A_{ii}}\\left(b_i-\\sum_{j=0}^{i-1} A_{ij}x_{j}^{(k)}-\\sum_{j=i+1}^{n-1} A_{ij}x_{j}^{(k-1)}\\right)$$\n",
    "\n",
    "Vsoto smo razdelili na dva dela in za izračun $i$-tega člena upoštevali v $k$-ti iteraciji že določene člene z indeksom manjšim od $i$.\n",
    "\n",
    "Iterativni pristop prekinemo, ko dosežemo želeno natančnost rešitve $\\epsilon$:\n",
    "\n",
    "$$||x_{i}^{(k)}-x_{i}^{(k-1)}||<\\epsilon$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Zgled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[8, -1, 1],\n",
    "              [-1, 6, -1],\n",
    "              [0, -1, 6]], dtype=float) \n",
    "b = np.array([-14, 36, 6], dtype=float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Začetni približek:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.,  0.,  0.])"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.zeros(len(A))\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Pripravimo matriko $\\mathbf{A}$ brez diagonalnih elementov (Zakaj? Poskusite odgovoriti spodaj, ko bomo izvedli iteracije.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = A.copy() #naredimo kopijo, da ne povozimo podatkov\n",
    "np.fill_diagonal(K, np.zeros(3)) #spremenimo samo diagonalne elemente"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0., -1.,  1.],\n",
       "       [-1.,  0., -1.],\n",
       "       [ 0., -1.,  0.]])"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "K"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Izvedemo iteracije:\n",
    "$$x_{i}^{(k)} =\\frac{1}{A_{ii}}\\left(b_i-\\sum_{j=0}^{i-1} A_{ij}x_{j}^{(k)}-\\sum_{j=i+1}^{n-1} A_{ij}x_{j}^{(k-1)}\\right)$$\n",
    "(Ker bomo vrednosti takoj zapisali v `x`, ni treba razbiti vsote na dva dela)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-----iteracja 0-----\n",
      "Približek za element 0 [-1.75  0.    0.  ]\n",
      "Približek za element 1 [-1.75        5.70833333  0.        ]\n",
      "Približek za element 2 [-1.75        5.70833333  1.95138889]\n",
      "Norma 6.281360365408393\n",
      "-----iteracja 1-----\n",
      "Približek za element 0 [-1.28038194  5.70833333  1.95138889]\n",
      "Približek za element 1 [-1.28038194  6.11183449  1.95138889]\n",
      "Približek za element 2 [-1.28038194  6.11183449  2.01863908]\n",
      "Norma 0.6227976321231091\n",
      "-----iteracja 2-----\n",
      "Približek za element 0 [-1.23835057  6.11183449  2.01863908]\n",
      "Približek za element 1 [-1.23835057  6.13004808  2.01863908]\n",
      "Približek za element 2 [-1.23835057  6.13004808  2.02167468]\n",
      "Norma 0.04590845211691733\n"
     ]
    }
   ],
   "source": [
    "for k in range(3):\n",
    "    xk_1 = x.copy()\n",
    "    print(5*'-' + f'iteracja {k}' + 5*'-')\n",
    "    for i in range(len(A)): #opazujte kaj se dogaja, ko to celico poženete večkrat!\n",
    "        x[i] = (b[i]-K[i,:].dot(x))/A[i,i]\n",
    "        print(f'Približek za element {i}', x)\n",
    "    e = np.linalg.norm(x-xk_1)\n",
    "    print(f'Norma {e}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Preverimo rešitev"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-14.01517799,  35.9969644 ,   6.        ])"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A@x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Metoda deluje dobro, če je matrika diagonalno dominantna (obstajajo pa metode, ki delujejo tudi, ko matrika ni diagonalno dominantna, glejte npr.: J. Petrišič, Reševanje enačb, 1996, str 149: Metoda konjugiranih gradientov)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Nekaj vprašanj za razmislek!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Na sliki je prikazano paličje. Ob delovanju sil $F_1$ in $F_2$ se v palicah razvijejo notranje sile $P_i$. Dimenzije paličja zagotavljata kota $\\alpha$ in $\\beta$.\n",
    "<img src=\"./fig/palicje_2.2_15.png\" width=300>\n",
    "Sile v palicah izračunamo s pomočjo sistema linearnih enačb.\n",
    "\n",
    "1. V simbolni obliki zapišite ravnotežje sil za točko 1 v $x$ in $y$ smeri (namig: naloga je posplošitev naloge 15 na strani 81 v knjigi Numerical methods in Eng with Py 3 z nastavkom za rešitev):\n",
    "<img src=\"./fig/palicje_2.2_15_resitev.png\" width=300>\n",
    "    Zgornji nastavek ima napako v predzadnji vrstici. Rešitev sistema v knjigi je: $P_{1} = -42000, P_{2}=-12000, P_{3}= 42426, P_{4} = -12000, P_{5} = 16971, P_{6} = 12000$\n",
    ").\n",
    "2. V simbolni obliki zapišite ravnotežje sil za točko 2 v $x$ in $y$ smeri.\n",
    "3. Najdite simbolno rešitev za sile $P_i$.\n",
    "4. Uporabite podatke: $\\alpha=\\beta=\\pi/4$, $F_1=18$ kN in $F_2=12$ kN ter najdite številčno rešitev.\n",
    "5. Pripravite si funkcijo, ki bo za poljubne podatke (npr: ``podatki = {a: pi/4, b: pi/4, F1: 18000, F2: 12000}``) vrnila numerično matriko koeficientov $\\mathbf{A}$ in vektor konstant $\\mathbf{b}$. Če ne uspete tega narediti avtomatizirano, delajte \"na roke\" (splača se vam potruditi, saj bomo to večkrat rabili).\n",
    "6. Razširite zgornjo funkcijo, da vam vrne rešitev linearnega sistema (uporabite kar ``numpy`` knjižnico)\n",
    "7. Predpostavite $F_1=F_2=10$ kN. V vsaj petih vrednostih kota $\\alpha=\\beta$ od 10$^{\\circ}$ do 80$^{\\circ}$ izračunajte sile v palicah.\n",
    "8. Za primer iz predhodne naloge narišite sile v palicah.\n",
    "9. S pomočjo funkcije ``np.linalg.solve`` izračunajte inverz poljubne matrike ``A`` (nato izračunajte še inverz s pomočjo funkcije ``np.linalg.inv``).\n",
    "10. Na primeru poljubnih podatkov (npr: ``podatki = {a: pi/4, b: pi/4, F1: 18000, F2: 12000}``) pokažite Gaussovo eliminacijo z delnim pivotiranjem.\n",
    "11. Na primeru poljubnih podatkov (npr: ``podatki = {a: pi/4, b: pi/4, F1: 18000, F2: 12000}``) pokažite Gauss-Seidlov iterativni pristop k iskanju rešitve."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Dodatno"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analizirajte, koliko e-mailov dobite na dan:\n",
    "<img src=\"./fig/povprecje_emailov.jpg\" width=300>\n",
    "https://plot.ly/ipython-notebooks/graph-gmail-inbox-data/\n",
    "(Nasvet: sledite kodi in uporabite ``matplotlib`` za prikaz)."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "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.6.3"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  },
  "livereveal": {
   "theme": "sky",
   "transition": "default",
   "width": 1200
  },
  "nav_menu": {},
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "281px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": "2",
   "toc_cell": true,
   "toc_section_display": "none",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}