03-SistemasLineares-I.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Resolução de Sistemas Lineares\n",
    "\n",
    "Vamos considerar aqui os métodos para solução de sistemas lineares com $n$ equações e $n$ incógnitas da forma\n",
    "$$ \\begin{array}{ccccc}\n",
    "a_{11}x_1 & + a_{12}x_2 & + \\cdots & + a_{1n}x_n& = b_1 \\\\\n",
    "\\vdots & & & &  \\vdots \\\\\n",
    "a_{n1}x_1& + a_{n2}x_2& + \\cdots & + a_{nn}x_n & = b_n\n",
    "\\end{array}$$\n",
    "ou na forma matricial:\n",
    "$$ \\begin{bmatrix}\n",
    "a_{11} & a_{12} & \\cdots & a_{1n} \\\\\n",
    "\\vdots & \\vdots & & \\vdots \\\\\n",
    "a_{n1} & a_{n2} &\\cdots & a_{nn}\n",
    "\\end{bmatrix}\\begin{pmatrix} x_1 \\\\ \\vdots \\\\ x_n \\end{pmatrix} = \\begin{pmatrix} b_1 \\\\ \\vdots \\\\ b_n \\end{pmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sistemas fáceis de resolver\n",
    "Em primeiro lugar caracterizamos os sistemas fáceis de resolver. Um sistema é **fácil de resolver** quando a matriz de coeficientes $A$ é triangular superior, isto é quando tivermos:\n",
    "$$ a_{ij}=0 \\text{ para } i > j. $$ Neste caso resolvemos o sistema começando da última equação e chegando até a primeira incógnita por substituição. Claro que o sistema só será possivel e determinado quando cada elemento $a_{ii}$ da diagonal for diferente de zero. E podemos aplicar a fórmula:\n",
    "$$ x_k = \\frac{1}{a_{kk}}\\left(b_k - \\sum_{k < j \\leq n} a_{kj}x_j\\right) $$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "# para a maioria dos exemplos vamos usar o numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def BackStep(A,b):\n",
    "    \"\"\" Se a matriz for triangular superior resolve Ax=b \n",
    "    A é uma matriz triangular.\n",
    "    b é uma lista de números reais de mesma dimensão de A. \"\"\"\n",
    "    #i Vamos testar se A é TS e se a a diagonal não tem zero\n",
    "    # Esta parte consome muito tempo e não deve ser usada depois\n",
    "    n = len(A)\n",
    "    for i in range(n):\n",
    "        teste1 = A[i,i]==0\n",
    "        if i==n:\n",
    "            teste2 = False\n",
    "        else: teste2 = False in (A[i+1:n,i]==0)\n",
    "        if teste1 or teste2:\n",
    "            raise ValueError(\" Matriz A não pode ser parametro do BackStep \")\n",
    "    # Esta foi só a parte de teste.\n",
    "    x=np.zeros(n)\n",
    "    k=n-1\n",
    "    while k>=0:\n",
    "        x[k] = (1/A[k,k])*(b[k]-sum([A[k,j]*x[j] for j in range(k+1,n)]))\n",
    "        k-=1\n",
    "    return x\n",
    "                           "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transformar um sistema qualquer num sistema fácil de resolver\n",
    "O método direto consiste em transformar um sistema $Ax=b$ num sistema $A_1x=b_1$ sem alterar as soluções do sistema e tal que o segundo sistema seja fácil de resolver.\n",
    "\n",
    "### Operações de linhas elementares\n",
    "São as operações básicas que não alteram o conjunto solução e que iremos utilizar.\n",
    "\n",
    "*  $E_{ij}$ Significa trocar as posições das equações $j$ e $i$. Na matriz de coeficientes e  no vetor $b$ corresponde a trocar estas linhas.\n",
    "\n",
    "* $E_{ij,\\alpha}$: isto corresponde a fazer a substituição $L_i = L_i + \\alpha L_j$, ou seja, somar a equação $i$ com um múltiplo da equação $j$\n",
    "\n",
    "* $E_{i,\\beta}$ para $\\beta \\neq 0$ correspode a multiplicação da $i$-ésima equação por $\\beta$ ($=E_{ii,{\\beta-1}}$).\n",
    "\n",
    "### Algoritmo da Eliminação de Gauss\n",
    "\n",
    "Para o $k$ de $1$ até $n-1$ executamos:\n",
    "\n",
    "*Passo $k$:* Caso $a_{kk} \\neq 0$, faça de $ i=k+1$ até $ i=n $: $ m_{ik} = a_{ik}/a_{kk} $ e aplique ao  sistema a operação elementar $ E_{ik,m_{ik}}$.\n",
    "\n",
    "Ao fim do passo $n-1$ o sistema estará triangularizado.\n",
    "Para evitar a possibilidade de não se triangularizar o sistema, mesmo ele sendo possível e determinado, executamos a **pivotação** que é o seguinte procedimento:\n",
    "antes da execução do passo $k$ do MEG. Buscamos, na $k-$ésima coluna, a partir da $k-$ésima linha o índice $i$\n",
    "tal que $|a_{ik}|$ seja o maior possível, e executamos $E_{ik}$. Com isto garantimos a triangularização de qualquer sistema que tenha solução única e evitamos perdas de algarismos significativos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "# A is type array from numpy\n",
    "# remember: A[0.0] is the first element\n",
    "# Elementary row operations:\n",
    "\n",
    "def  troca_linha(A,i,j):\n",
    "    '''Troca as linhas i e j da matriz A'''\n",
    "    buffer = A[i].copy()\n",
    "    A[i] = A[j]\n",
    "    A[j] = buffer\n",
    "    return A\n",
    "\n",
    "def mult_linha(A,i,alfa):\n",
    "    '''Multiplica por alfa a linha i'''\n",
    "    A[i] = alfa*A[i]\n",
    "    return A\n",
    "\n",
    "def subs_linha(A,i,k,alfa):\n",
    "    '''soma a linha i com um multiplo da linha k'''\n",
    "    A[i] = A[i] + alfa*A[k]\n",
    "    return A\n",
    "\n",
    "# step k of gauss elimination:\n",
    "def gauss_step(A,k):\n",
    "    '''O k-esimo passo na eliminacao de gauss'''\n",
    "    L=range(len(A))\n",
    "    for j in L[k+1:]:\n",
    "        m=-A[j,k]/A[k,k]\n",
    "        A=subs_linha(A,j,k,m)\n",
    "        A[j,k]=-m # Aproveitamos a matriz A pra guardar o multiplicador\n",
    "    return A\n",
    "\n",
    "# triangle form, without care and pivoting\n",
    "def triangle_form(A):\n",
    "    '''Retorna a matriz escalonada,\n",
    "    com os multiplicadores na parte triangular inferior'''\n",
    "    L=len(A)\n",
    "    for k in range(L-1):\n",
    "        p=pivo(A,k)\n",
    "        troca_linha(A,k,p)\n",
    "        gauss_step(A,k)\n",
    "     \n",
    "    return A\n",
    "   \n",
    "def pivo(A,n):\n",
    "    ''' retorna o indice de pivotacao da matriz A na coluna n,\n",
    "    a partir da linha n'''\n",
    "    tamanho = shape(A)\n",
    "    l=tamanho[0] # numero de linhas de A\n",
    "    c= tamanho[1] # numero de colunas\n",
    "    if (l<=n) or (c<=n):\n",
    "        raise ValueError(\"n deve ser menor que dimensao de A\")\n",
    "    p=n # inicio o pivo como n, e vou aumentando\n",
    "    for k in range(n,l):\n",
    "        if abs(A[k,n]) > abs(A[p,n]):\n",
    "            p=k\n",
    "    return p\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1 2 3]\n",
      " [0 1 2]\n",
      " [0 0 3]]\n",
      "[-2.  0.  1.]\n"
     ]
    }
   ],
   "source": [
    "#Teste\n",
    "# do Backstep\n",
    "A=np.matrix([[1,2,3],[0,1,2],[0,0,3]])\n",
    "print(A)\n",
    "b=[1,2,3]\n",
    "x=BackStep(A,b)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# teste da função troca_linha\n",
    "A1=troca_linha(A,0,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[0, 1, 2],\n",
       "        [1, 2, 3],\n",
       "        [0, 0, 3]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[0, 1, 2],\n",
       "        [1, 2, 3],\n",
       "        [0, 0, 3]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A2 = troca_linha(A,1,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[0, 1, 2],\n",
       "        [0, 0, 3],\n",
       "        [1, 2, 3]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[0, 1, 2],\n",
       "        [0, 0, 3],\n",
       "        [1, 2, 3]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A3 = mult_linha(A,1,1/3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[0, 1, 2],\n",
       "        [0, 0, 1],\n",
       "        [1, 2, 3]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# teste da função de pivotação\n",
    "pivo(A,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[1, 2, 3],\n",
       "        [0, 1, 2],\n",
       "        [0, 0, 1]])"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#teste da função de triangularização\n",
    "triangle_form(A)"
   ]
  },
  {
   "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
}