03-SistemasLineares-III.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Métodos Iterativos para a solução de sitemas lineares \n",
    "\n",
    "A idéia dos métodos iterativos é produzir uma sequência de estimativas vetoriais $ \\mathbf{x}^{(k)} $, que idealmente convirja para a solução $\\mathbf{x}$ e\n",
    "cada passo é calculado a partir da estimativa anterior. Costumam ser bem mais rápidos que oe métodos diretos no caso de matrizes muito grandes.\n",
    "Para construir recursivamente a sequência que dá origem ao método iterativo fazemos uma decomposição da matriz $A$ de coeficientes do sistema como\n",
    "$$ A = N - P$$\n",
    "Onde $N$ é de uma classe de sistemas inversíveis e fáceis de resolver (ou seja, fáceis de inverter). Então resolver o sistema:\n",
    "$$ A\\mathbf{x} =b $$ \n",
    "é equivalente a resolver a equação:\n",
    "$$ N\\mathbf{x} = b+ P\\mathbf{x} \\text{ ou } \\mathbf{x}=N^{-1}(b+ P\\mathbf{x})$$\n",
    "A recursão\n",
    "$$ N\\mathbf{x}^{(k+1)} = b+ P\\mathbf{x}^{(k)}$$\n",
    "é a que converge para a solução caso $ ||N^{-1}P|| < 1$ onde a norma de matrizes é induzida de uma norma em $\\mathbb{R}^n$ (*veja no apêndice*).\n",
    "Destacamos dois casos\n",
    "\n",
    "*  Método de Jacobi: $ N = \\text{diag}(a_{11}, \\cdots, a_{nn})$\n",
    "\n",
    "*  Método de Gauss-Seidel : $ N = (n_{ij}) \\text{ com } n_{ij} = a_{ij} \\text{ se } i\\geq j \\\\ n_{ij}=0 \\text{ se } i\\lt j$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Método de Gauss-Seidel\n",
    "Vamos ver com mais detalhes só o método de Gauss-Seidel. Vamos denotar por $x_i^{(k)}$ a $i$-ésima coordenada do vetor estimativa, no passo $k$. Então podemos escrever:\n",
    "$$ \\begin{eqnarray}\n",
    "x_1^{(k+1)} & = &\\frac{1}{a_{11}}\\left( b_1 - \\sum_{j=2}^n a_{1j}x_j^{(k)}\\right) \\\\\n",
    "x_2^{(k+1)} & = &\\frac{1}{a_{22}}\\left( b_2 -  a_{21}x_1^{(k+1)} -\\sum_{j=3}^n a_{2j}x_j^{(k)}\\right) \\\\\n",
    "  & \\vdots & \\\\\n",
    "x_l^{(k+1)} & = &  \\frac{1}{a_{ll}}\\left( b_l -  \\sum_{j=1}^{l-1}a_{lj}x_j^{(k+1)} -\\sum_{j=l+1}^n a_{lj}x_j^{(k)}\\right)\\\\\n",
    "  & \\vdots & \\\\\n",
    " x_n^{(k+1)} & = &  \\frac{1}{a_{nn}}\\left( b_n -  \\sum_{j=1}^{n-1}a_{lj}x_j^{(k+1)}\\right)\n",
    "\\end{eqnarray}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(matrix([[ 0.03039834],\n",
      "        [ 0.32285115],\n",
      "        [ 0.38784066]]), 16, 3.5604346043801354e-08)\n"
     ]
    }
   ],
   "source": [
    "from pylab import *\n",
    "\n",
    "# dou a matriz quadrada de coeficientes e vetor independente\n",
    "\n",
    "A=matrix('4 2 0.6; 2 5.5 3; 1 0 2.5')\n",
    "b=matrix('1. ; 3.; 1.')\n",
    "\n",
    "def iterGS(x, A, b):\n",
    "    xold=x.copy()\n",
    "    L=list(x)\n",
    "    for i in range(len(A)):\n",
    "        soma = 0.\n",
    "        for j in range(len(A)):\n",
    "            if (j != i): soma += -A[i,j]*L[j]\n",
    "        L[i] = (b[i] + soma)/A[i,i]\n",
    "    return L\n",
    "\n",
    "def GaussSeidel(x0,A,b,eps,n):\n",
    "    IterGauss=0\n",
    "    dx=10.\n",
    "    while ((IterGauss<=n)&(dx>eps)):\n",
    "        xold = x0.copy()\n",
    "        L=iterGS(x0, A, b)\n",
    "        for k in range(len(L)) : L[k]=float(L[k])\n",
    "        x0=matrix(L).T\n",
    "        IterGauss += 1\n",
    "        dx = norm(x0-xold)\n",
    "    return x0 , IterGauss, dx\n",
    "\n",
    "########## Teste com um exemplo ##############\n",
    "\n",
    "x0=matrix(' 0; 0; 0 ')\n",
    "\n",
    "sol = GaussSeidel(x0,A,b,0.0000001, 100)\n",
    "\n",
    "print (sol)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'https://gist.github.com/9784091'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%pastebin '03-SistemasLineares-III.ipynb'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}