6-Sistemas Lineares.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "% matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sistemas lineares\n",
    "\n",
    "Nesta aula exploraremos alguns métodos de resolução de sistemas lineares. Nosso objetivo é conseguir a solução para um sistema de n equações lineares e n incógnitas similar ao sistema abaixo\n",
    "\n",
    "\\begin{eqnarray}\n",
    "x &+ 2y &+ 3z &= 4 \\\\\n",
    "x &+ 3y &+ 3z &= 3 \\\\\n",
    "x &+ 2y &+ z &= 0\n",
    "\\end{eqnarray}\n",
    "\n",
    "A abordagem intuitiva consiste em fazer simplificações a partir de combinações das diferentes linhas do sistema de equações. No caso acima, por exemplo, poderíamos avançar subtraindo a primeira equação menos a terceira ou a segunda menos a primeira. Isto gera as novas equações\n",
    "\n",
    "\\begin{eqnarray}\n",
    "2z &= 4 \\\\\n",
    "y &= -1\n",
    "\\end{eqnarray}\n",
    "\n",
    "É fácil ver que $z = 2$ e, por substituição, x = 0. Simples, né?\n",
    "\n",
    "Talvez não. O método intuitivo é eficaz em sistemas de equações pequenos. No então, é comum encontrarmos problemas que envolvam literalmente milhares de variáveis e qualquer tentativa de solução manual pode ser descartada. Precisamos de um método mecânico que possa ser transportado facilmente para um programa de computador.\n",
    "\n",
    "Esta aula irá apresentar algumas opções."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sistemas lineares e matrizes\n",
    "\n",
    "O primeiro passo para resolver sistematicamente um sistema linear, é escrevê-lo na forma matricial. Isto elimina qualquer menção às variáveis desconhecidas e permite que trabalhemos apenas com números, que é algo que o computador sabe fazer bem:\n",
    "\n",
    "\\begin{equation}\n",
    "\\left[\n",
    "    \\begin{array}{ccc}\n",
    "    1 & 2 & 3 \\\\\n",
    "    1 & 2 & 1 \\\\\n",
    "    1 & 3 & 3\n",
    "    \\end{array}\n",
    "\\right]\\times\n",
    "\\left[\\begin{array}{c}x\\\\y\\\\z\\end{array}\\right]=\n",
    "\\left[\\begin{array}{c}4\\\\3\\\\0\\end{array}\\right]\n",
    "\\end{equation}\n",
    "\n",
    "(Se você não se lembra como faz isto, use a [wikipedia](https://pt.wikipedia.org/wiki/Sistema_de_equa%C3%A7%C3%B5es_lineares))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrizes em Python\n",
    "\n",
    "O Numpy possui duas estruturas de dados semelhantes que podem ser utilizadas para representar matrizes. Uma delas é a nossa velha conhecida array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Uma matriz é criada a partir de lista de listas\n",
    "# Alinhamos o como ele estaria disposto em uma matriz\n",
    "# para facilitar a visualização\n",
    "M = np.array([[1, 2, 3], \n",
    "              [1, 3, 3],\n",
    "              [1, 2, 1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Obtemos a transposta de uma matriz utilizando o atributo `.T` da mesma. Para multiplicar duas matrizes, é necessário utilizar a função `np.dot`. A multiplicação de duas arrays é feita termo a termo e não segue as regras de multiplicação de matrizes em álgebra linear. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 1, 1],\n",
       "       [2, 3, 2],\n",
       "       [3, 3, 1]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M.T  # transposta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[14, 16,  8],\n",
       "       [16, 19, 10],\n",
       "       [ 8, 10,  6]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(M, M.T)  # multiplicação de matrizes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 2, 3],\n",
       "       [2, 9, 6],\n",
       "       [3, 6, 1]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M * M.T  # multiplicação termo a termo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Outra opção é utilizar o tipo `np.matrix`. Uma matriz do numpy é muito semelhante a um array, mas se distingue por tratar o símbolo de multiplicação `*` como sendo a multiplicação matricial. É possível criar um objeto do tipo matrix da mesma forma que criamos o array (a partir de uma lista de listas) ou chamando um array já inicializado."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "M1 = np.array([[1, 2, 3],\n",
    "               [1, 3, 3],\n",
    "               [1, 2, 1]])\n",
    "M2 = np.matrix(M.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe que a multiplicação entre M1 e M2 é feita utilizando as regras de multiplicação matricial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[14, 16,  8],\n",
       "        [16, 19, 10],\n",
       "        [ 8, 10,  6]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M1 * M2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A função np.dot continua realizando a multiplicação matricial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[14, 16,  8],\n",
       "        [16, 19, 10],\n",
       "        [ 8, 10,  6]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(M1, M2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A escolha entre matrizes e arrays é um tanto pessoal. Fora o comportamento na multiplicação, os dois objetos se comportam de maneiras quase idênticas. Eu particularmente prefiro utilizar arrays porque não acho que a função np.dot seja inconveniente o suficiente para justificar aprender um novo tipo de dados. Arrays e matrizes as vezes se distiguem de formas sutis e prefiro evitar supresas utilizando somente arrays no código inteiro. É assim que seguiremos neste tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Escalonamento de matrizes\n",
    "\n",
    "Escalonamento consiste na realização de combinações lineares entre as linhas de uma matriz para transformá-la numa forma diagonal. Começaremos com o escalonamento para os casos mais simples (já já veremos algumas situações onde o método simples pode ser problemático). \n",
    "\n",
    "Nosso objetivo é rearranjar as linhas para que possamos escrever a matriz na forma diagonal superior. Esta forma é fácil de resolver por substituição. Antes vamos nos lembrar da nossa matriz M:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  2.,  3.],\n",
       "       [ 1.,  3.,  3.],\n",
       "       [ 1.,  2.,  1.]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = M.astype(float); M  # convertemos os números para float para evitar truncamentos errados"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "É necessário também definir o vetor que ocupa o lado direito do nosso sistema linear. Pensamos no sistema linear como sendo $Mu = v$, onde $M$ é uma matriz, $u$ é o vetor solução desconhecido e $v$ é o vetor conhecido no lado direito da equação."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Mu = v\n",
    "v = np.array([4., 3., 0.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "O primeiro passo antes de combinar as linhas da matriz é fazer a matriz extendida: anexamos o vetor v como sendo a última coluna da matriz. Desta forma, quando fazemos uma operação do tipo \"somar linha 1 com a 2\" ela se aplica tanto ao lado esquerdo quanto ao lado direito das equações correspondentes.\n",
    "\n",
    "O numpy possui a função hstack que \"empilha\" as matrizes horizontalmente."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  2.,  3.,  1.,  2.,  3.],\n",
       "       [ 1.,  3.,  3.,  1.,  3.,  3.],\n",
       "       [ 1.,  2.,  1.,  1.,  2.,  1.]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.hstack([M, M])   # mas *não* é isso que queremos :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se utilizarmos o hstack com [M, v], obtemos um erro: o formato de M (3, 3)  e v (3,) não são compatíveis. Precisamos transformar v em um vetor coluna. Vetor coluna nada mais é que uma matriz n x 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 4.],\n",
       "       [ 3.],\n",
       "       [ 0.]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v[:, None]  # None cria dimensões extras quando elas são necessárias"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Então colocamos o resultado na função hstack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  2.,  3.,  4.],\n",
       "       [ 1.,  3.,  3.,  3.],\n",
       "       [ 1.,  2.,  1.,  0.]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Mext = np.hstack([M, v[:, None]]); Mext"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Agora sim!\n",
    "\n",
    "O próximo passo é zerar os termos abaixo do elemento Mext[0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  2.,  3.,  4.],\n",
       "       [ 0.,  1.,  0., -1.],\n",
       "       [ 0.,  0., -2., -4.]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pivo = Mext[0, 0]\n",
    "Mext[1] = Mext[1] - Mext[1, 0] / pivo * Mext[0]\n",
    "Mext[2] = Mext[2] - Mext[2, 0] / pivo * Mext[0]\n",
    "Mext"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A idéia básica das equações acima é subtrair de cada linha a linha do pivo multiplicada pelo termo A/pivo, onde A é o valor que queremos eliminar. Em notação matemática, isto seria \n",
    "\n",
    "$$L_{ij} \\mapsto L_{ij} - \\frac{L_{i0}}{L_{00}} L_0j,$$\n",
    "\n",
    "onde $L_{0j}$ é a linha do pivo para a primeira iteração."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  2.,  3.,  4.],\n",
       "       [ 0.,  1.,  0., -1.],\n",
       "       [ 0.,  0., -2., -4.]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Continuaríamos usando M[1, 1] como pivo, mas a matriz já está escalonada!\n",
    "pivo = Mext[1, 1]\n",
    "Mext[2] = Mext[2] - Mext[2, 1] / pivo * Mext[1]\n",
    "Mext"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Traduzimos de volta para um sistema linear e vemos que a resposta é fácil de obter:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "x &+ 2y &+ 3z &= 4 \\\\\n",
    "0 &+ y  &+ 0 &= -1 \\\\\n",
    "0 &+ 0  &- 2z &= -4\n",
    "\\end{eqnarray}\n",
    "\n",
    "Onde encontramos novamente que $z=2$, $y=-1$ e $x=0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideOutput": true
   },
   "source": [
    "### Generalizando\n",
    "\n",
    "Generalizamos o procedimento para uma matriz aleatória de dimensão m x n (onde normalmente, n = m + 1, mas na verdade basta que n > m para o algoritmo funcionar)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.48198161,  0.98422871,  0.82976294,  0.33444333],\n",
       "       [ 0.4669344 ,  0.16256605,  0.32366328,  0.41938237],\n",
       "       [ 0.47518896,  0.45080182,  0.12961868,  0.52464709]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m, n = 3, 4\n",
    "matriz = np.random.uniform(size=(m, n))\n",
    "matriz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.48198161,  0.98422871,  0.82976294,  0.33444333],\n",
       "       [ 0.        , -0.79093558, -0.48019491,  0.09538018],\n",
       "       [ 0.        ,  0.        , -0.37301607,  0.13226304]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Escalonamento\n",
    "for k in range(m - 1):\n",
    "    pivo = matriz[k, k]\n",
    "    for i in range(k + 1, m):\n",
    "        matriz[i] = matriz[i] - matriz[i, k] / pivo * matriz[k]\n",
    "matriz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A matriz triangular fica ainda mais simples de trabalhar se o pivô for sempre igual a 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  2.04204622,  1.72156558,  0.69389231],\n",
       "       [-0.        ,  1.        ,  0.60712266, -0.12059159],\n",
       "       [-0.        , -0.        ,  1.        , -0.35457733]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Normalização\n",
    "for k in range(m):\n",
    "    matriz[k] = matriz[k] / matriz[k, k]\n",
    "matriz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Daí temos imediatamente o valor de z. Podemos calcular o valor das outras coordenadas com o procedimento:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "u = np.zeros(m)\n",
    "u[2] = matriz[2, m]  # este é o valor de z\n",
    "u[1] = matriz[1, m] - sum(matriz[1, 0:m] * u)  # agora y\n",
    "u[0] = matriz[0, m] - sum(matriz[0, 0:m] * u)  # e finalmente z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A implementação genérica seria: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.11097881,  0.09468034, -0.35457733])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u = np.zeros(m)\n",
    "for k in reversed(range(m)):\n",
    "    u[k] = matriz[k, m] - sum(matriz[k, 0:m] * u)\n",
    "u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -1.11022302e-16,   0.00000000e+00,   0.00000000e+00])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Conferimos a resposta (quase zero)\n",
    "np.dot(matriz[:, 0:m], u) - matriz[:, -1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Juntando os pedaços: eliminação de Gauss\n",
    "\n",
    "Abaixo segue um código que realiza todas as etapas acima de uma vez"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def resolver_gauss(M, v):\n",
    "    \"\"\"Retorna o vetor solução x para o sistema linear Mx = v\"\"\"\n",
    "    \n",
    "    m = len(M)\n",
    "    \n",
    "    # Cria matriz extendida de floats\n",
    "    v = np.asarray(v, dtype='float')\n",
    "    M = np.hstack([M, v[:, None]])\n",
    "    \n",
    "    # Escalona matriz\n",
    "    for k in range(m - 1):\n",
    "        pivo = M[k, k]\n",
    "        for i in range(k + 1, m):\n",
    "            M[i] = M[i] - M[i, k] / pivo * M[k]\n",
    "            \n",
    "    # Normaliza pivos\n",
    "    for k in range(m):\n",
    "        M[k] = M[k] / M[k, k]\n",
    "        \n",
    "    # Resolve sistema triangular superior\n",
    "    x = np.zeros(m)\n",
    "    for k in reversed(range(m)):\n",
    "        x[k] = M[k, m] - sum(M[k, 0:m] * x)\n",
    "    return x\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Podemos comparar com a rotina utilizada pelo numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 0.2,  0. ,  0.2]),\n",
       " array([  2.00000000e-01,  -1.11022302e-16,   2.00000000e-01]))"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = np.array([[1, 2, 4], \n",
    "              [3, 5, 7], \n",
    "              [6, 8, 9]])\n",
    "v = np.array([1, 2, 3])\n",
    "\n",
    "resolver_gauss(M, v), np.linalg.solve(M, v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vemos que as duas soluções são praticamente idênticas. Vamos comparar a performance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eliminação de Gauss\n",
      "The slowest run took 4.88 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "10000 loops, best of 3: 51.9 µs per loop\n",
      "\n",
      "\n",
      "Numpy\n",
      "The slowest run took 4.15 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "100000 loops, best of 3: 16.6 µs per loop\n"
     ]
    }
   ],
   "source": [
    "print('Eliminação de Gauss')\n",
    "%timeit resolver_gauss(M, v)\n",
    "\n",
    "print('\\n\\nNumpy')\n",
    "%timeit np.linalg.solve(M, v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "O numpy é consideravelmente mais rápido pois todo processo é implementado de forma bastante otimizada em FORTRAN."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Notas finais\n",
    "\n",
    "Não cobrimos aqui o importante processo de pivoteamento mencionado em sala de aula. Seguem algumas referências:\n",
    "\n",
    "* http://www.mat.ufrgs.br/~fabio/linear.pdf\n",
    "* http://conteudo.icmc.usp.br/pessoas/andretta/ensino/aulas/sme0301-1-11/SistemasLinearesGaussPivoteamento.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Método de Jacobi\n",
    "\n",
    "Este é outro método para resolução de sistemas lineares. Trata-se de um método iterativo que produz sucessivas aproximações da solução candidata.\n",
    "\n",
    "Tudo parte da seguinte fatoração das matrizes:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "      Mu &= v, \\\\\n",
    "(D + R)u &= v,\n",
    "\\end{eqnarray}\n",
    "\n",
    "onde $D$ é uma submatriz de M com inversa conhecida e $R = M - A$ é o resto. Reorganizando os termos, temos\n",
    "\n",
    "\\begin{eqnarray}\n",
    "& Du &= v - Ru, \\\\\n",
    "& u &= D^{-1}(v - Ru),\n",
    "\\end{eqnarray}\n",
    "\n",
    "Podemos olhar a equação de baixo e tentar utilizá-la como uma regra de iteração\n",
    "\n",
    "$$u_{k+1} = D^{-1}(v - Ru_k)$$\n",
    "\n",
    "A idéia é que se $u_k$ convergir obteremos o valor da solução do sistema linear. No método de Jacobi, escolhemos\n",
    "a matriz $D$ como sendo a formada apenas pelos elementos da diagonal. Uma matriz diagonal pode ser invertida trivialmente e é possível aplicar a regra acima."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### Exemplo\n",
    "\n",
    "Definimos matrizes e o vetor v. O sistema linear é Mu = v."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Matriz\n",
    "M = np.array([[10, 5, 2],\n",
    "              [-1, 2, 1],\n",
    "              [0,  2, 5]])\n",
    "\n",
    "# Vetor\n",
    "v = np.array([1, 2, 3])\n",
    "\n",
    "# Matriz diagonal\n",
    "D = np.array([[10, 0, 0],\n",
    "              [ 0, 2, 0],\n",
    "              [ 0, 0, 5]])\n",
    "Dinv = np.array([[0.1,   0,   0],\n",
    "                 [  0, 0.5,   0],\n",
    "                 [  0,   0, 0.2]])\n",
    "\n",
    "# Resto\n",
    "R = M - D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Começamos com um chute inicial para u e repetimos a regra de iteração algumas vezes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.6  1.   0.2]\n",
      "[-0.44  0.6   0.2 ]\n",
      "[-0.24  0.68  0.36]\n",
      "[-0.312  0.7    0.328]\n",
      "[-0.3156  0.68    0.32  ]\n",
      "[-0.304   0.6822  0.328 ]\n",
      "[-0.3067   0.684    0.32712]\n",
      "[-0.307424  0.68309   0.3264  ]\n",
      "[-0.306825  0.683088  0.326764]\n",
      "[-0.3068968  0.6832055  0.3267648]\n"
     ]
    }
   ],
   "source": [
    "u = np.array([1, 1, 1])\n",
    "\n",
    "for _ in range(10):\n",
    "    u = np.dot(Dinv, v - np.dot(R, u))\n",
    "    print(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vemos que a resposta converge para um valor determinado.\n",
    "\n",
    "Agora vamos juntar tudo isso em uma função."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def jacobi(M, v, niter=10):\n",
    "    \"\"\"Realiza niter iterações do método de Jacobi.\"\"\"\n",
    "    \n",
    "    # Extraímos a diagonal da matriz como vetor.\n",
    "    diag = np.diag(M)\n",
    "    \n",
    "    # Se chamarmos np.diag num vetor, ela gera uma matriz\n",
    "    D = np.diag(diag)\n",
    "    Dinv = np.diag(1 / diag)\n",
    "    \n",
    "    # Criamos a matriz resto\n",
    "    R = M - D\n",
    "    \n",
    "    # Realizamos a n iterações do algoritmo\n",
    "    u = v\n",
    "    for _ in range(niter):\n",
    "        u = np.dot(Dinv, v - np.dot(R, u))\n",
    "    return u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.30707197,  0.68323194,  0.3266777 ])"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jacobi(M, v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe que nem sempre o método de jacobi converge (para falar a verdade, ele diverge para a maioria das matrizes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Gauss: [ 0.2  0.   0.2]\n",
      "Jacobi: [ 51189.24210402  23652.34524547  21626.36327457]\n"
     ]
    }
   ],
   "source": [
    "M = np.array([[1, 2, 4],\n",
    "              [3, 5, 7],\n",
    "              [6, 8, 9]])\n",
    "print(' Gauss:', resolver_gauss(M, v))\n",
    "print('Jacobi:', jacobi(M, v))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mencionamos o critério das linhas em sala de aula que diz quando uma matriz pode ser resolvida pelo método de Jacobi. (https://pt.wikipedia.org/wiki/M%C3%A9todo_de_Jacobi). O critério das linhas diz que o elemento da diagonal em cada linha da matriz deve ser maior em módulo que a soma dos outros elementos da linha (também em módulo). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}