Física Computacional - FSC-5705

só um divisor

Sistema de equações Lineares - Método de Gauss

Se denomina de sistema linear de equações ao conjunto $m$ equações, com $n$ incógnitas do tipo

\( \begin{array}{c c c c c c c c c c c c} E_1: & a_{11}x_1 &+& a_{12}x_2 &+& a_{13}x_3 &+& \cdots &+& a_{1n}x_n &=& b_1\\ E_2: & a_{21}x_1 &+& a_{22}x_2 &+& a_{23}x_3 &+& \cdots &+& a_{2n}x_n &=& b_2\\ E_3: & a_{31}x_1 &+& a_{32}x_2 &+& a_{33}x_3 &+& \cdots &+& a_{3n}x_n &=& b_3\\ & \vdots & & \vdots & & \vdots & & \vdots & & \vdots &=& \vdots\\ E_n: & a_{n1}x_1 &+& a_{n2}x_2 &+& a_{n3}x_3 &+& \cdots &+& a_{nn}x_n &=& b_n \end{array} \)

Observe que, ainda que o nome do sistema de equações indique que é linear, devemos entender que os $X_n$ representam uma variável, dessa forma podemos ter algo como $x_n=x^n$, isto é, as equações podem ser não lineares, mas o sistema pode ser tratado com sendo não linear.

A fim de resolver sistemas como esse existem diversos métodos, nestas notas nos concentraremos no método conhecido como método de Gauss. O método de Gauss está baseado em manipulações algebraicas simples realizadas ao sistema de equação a fim de deixar esse sistema na sua forma triangular ou reduzida, estas manipulações são:

  1. A equação $E_i$ pode ser multiplicada por qualquer constante $\lambda$ diferente de zero, e a equação resultante pode ser utilizada em lugar de $E_i$. Essa operação é indicada como $(\lambda E_i)\rightarrow (E_i)$
  2. A equação $E_j$ pode ser multiplicada por qualquer constante $\lambda$ e adicionada à equação $E_j$, e a equação resultante pode ser utilizada em lugar da equação $E_i$. Essa operação é indicada como $(E_i + \lambda E_j)\rightarrow (E_i)$
  3. As equações $E_i$ e $E_j$ podem ser transpostas em ordem. Essa operação é indicada por $(E_i) \rightarrow (E_j)$
Com esas simples regrinhas podemos levar o nosso sistema à forma triangular e assim obter mais simplesmente a resposta ao sistema desejado. Vejamos o seguinte exemplo

Exemplo

$\displaystyle{ \begin{array}{c c c c c c c c c} E_1: & x_1 & + & x_2 & & & + & 3x_4 & = & 4, \\ E_2: & 2x_1 & + & x_2 & - & x_3 & + & x_4 & = & 1, \\ E_3: & 3x_1 & - & x_2 & - & x_3 & + & 2x_4 & = & -3, \\ E_4: & -x_1 & + & 2x_2 & + & 3x_3 & - & x_4 & = & 4, \end{array} }$


$(E_2 - 2 E_1)\rightarrow (E_2)$, $(E_3 - 3 E_1)\rightarrow (E_3)$ e $(E_4 + E_1)\rightarrow (E_4)$

$\displaystyle{ \begin{array}{c c c c c c c c c} E_1: & x_1 & + & x_2 & & & + & 3x_4 & = & 4, \\ E_2: & & - & x_2 & - & x_3 & - & 5x_4 & = & -7, \\ E_3: & & - & 4x_2 & - & x_3 & - & 7x_4 & = & -15, \\ E_4: & & & 3x_2 & + & 3x_3 & + & 2x_4 & = & 8, \end{array} }$


$(E_3 - 4 E_2)\rightarrow (E_3)$ e $(E_4 + 3 E_2)\rightarrow (E_4)$

$\displaystyle{ \begin{array}{c c c c c c c c c} E_1: & x_1 & + & x_2 & & & + & 3x_4 & = & 4, \\ E_2: & & - & x_2 & - & x_3 & - & 5x_4 & = & -7, \\ E_3: & & - & & & 3x_3 & - & 13x_4 & = & 13, \\ E_4: & & & & & & - & 13x_4 & = & -13, \end{array} }$

De onde:
$x_4 = 1$,
$x_3 = 1/3(13 - 13x_4)=0$,
$x_2 = -(7+5x_4+x_3)=2$,
$x_1 = 4-3x_4-x_2=-1$

Resulta conveniente escrever o sistema na forma da matriz expandida:

$\displaystyle{ \left[ \begin{array}{c c c c c c c} a_{11} & a_{12} & a_{13} & \cdots & a_{1n} & & a_{1,n+1} \\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n} & & a_{2,n+1}\\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n} & & a_{3,n+1}\\ \vdots & \vdots & \vdots & \vdots & \vdots & & \vdots\\ a_{n1} & a_{n2} & a_{n3} & \cdots & a_{nn} & & a_{n,n+1} \end{array} \right] }$

Onde os $a_{i,n+1}$ são os valores de $b_i$.

A partir dessa matriz expandida podemos desenvolver um procedimento genérico que permita encontrar a solução de um sistema linear. Admitindo-se que $a_{ii}\neq 0$, devemos realizar sucessivas aplicações de $(E_j-(a_{ji}/a_{ii})E_i) \rightarrow (E_j)$, para cada $j = i+1, i+2, \ldots , n$, a fim de obter uma matriz triangular como

$\displaystyle{ \left[ \begin{array}{c c c c c c} a_{11} & a_{12} & \cdots & a_{1n} & & a_{1,n+1} \\ 0 & a_{22} & \cdots & a_{2n} & & a_{2,n+1}\\ \vdots & \ddots & \ddots & & & \vdots\\ \vdots & & \ddots & \ddots & & \vdots\\ 0 & \cdots & 0 & a_{nn} & & a_{n,n+1} \end{array} \right] }$

Esse resultado implica em que temos um sistema triangular

$\displaystyle{ \begin{array}{c c c c c c c c} a_{11}x_1 & a_{12}x_2 & + &\cdots & + & a_{1n}x_n & = & a_{1,n+1} \\ & a_{22}x_2 & + &\cdots & + & a_{2n}x_n & = & a_{2,n+1}\\ & & \ddots & & & \vdots & & \vdots\\ & & &\ddots & & \vdots & & \vdots \\ & & & &\ddots & \vdots & & \vdots \\ & & & & & a_{nn}x_n & = & a_{n,n+1} \end{array} }$

De onde se pode calcular diretamente

$\displaystyle{ x_n = \frac{a_{n,n+1}}{a_{nn}} }$

E, seguidamente

$\displaystyle{ x_{n-1} = \frac{a_{n-1,n+1} - a_{n-1,n}x_n}{a_{n-1,n-1}} }$

Assim, continuando o processo temos

$\displaystyle{ x_{i} = \frac{a_{i,n+1} - \sum _{j=i+1}^{n}a_{ij}x_j}{a_{ii}} }$

Essa procedimento funciona bem desde que nenhum dos $a_{ii}$ for igual a zero. Afortunamento existe um complemento ao procedimento que permite evitar esse problema. A ideia é trocar as linhas de forma a colocar o zeros do lado esquerda da diagonal principal até atingir a forma triangular desejada na matriz. A o método completo se lhe chama de eliminação de Gauss com pivotamento parcial

Algoritmo de eliminação de Gauss com substituição retroativa

Para resolver o sistema linear $n\times n$

$\displaystyle{ \begin{array}{c c c c c c c c c c c} a_{11}x_1 &+& a_{12}x_2 &+& a_{13}x_3 &+& \cdots &+& a_{1n}x_n &=& b_1\\ a_{21}x_1 &+& a_{22}x_2 &+& a_{23}x_3 &+& \cdots &+& a_{2n}x_n &=& b_2\\ \vdots & & \vdots & & \vdots & & \vdots & & \vdots &=& \vdots\\ a_{n1}x_1 &+& a_{n2}x_2 &+& a_{n3}x_3 &+& \cdots &+& a_{nn}x_n &=& b_n \end{array} }$

A fim de resolver um sistema como o anterior devemos construir duas subroutinas. A primeira subroutina levara a matriz estendida para sua forma diagonal superior mediante a aplicação das propriedades das matrizes acima descritas, a esse processo se lhe chama de eliminação progressiva pois os elementos por baixo do pivô são eliminados (suponhamos que a linha pivô fosse a linha $k$, consequentemente o elemento pivô será o elemento $a_{kk}$, a eliminação é realizada fazendo $E_{j} - \lambda \cdot E_k \rightarrow E_{j}$, onde $j \leq k$.

Algoritmo de eliminação

ENTRADA: número de equações $n$, matriz de coeficientes $a_{ij}$, e vetor $b_j$

Para $k=1,\; n-1$ faça
Para $i=k+1,\; n$ faça
$\lambda = \frac{a_{i,k}}{a_{k,k}}$
Para $j=k,\; n$ faça
$a_{i,j} = a_{i,j} - \lambda \cdot a_{k,j}$
fim do faça
$b_i = b_i - \lambda \cdot b_k$
fim do faça
fim do faça

Uma vez realizada a eliminação progressiva, estamos habilitados a calcular as solução para as diversas variáveis, para isso realizamos as substituições regressiva necessárias. Chamamos está etapa de substituição regressiva já que como vimos no exemplo anterior primeiro encontramos a solução para a variável de maior ordem $x_n$, com essa solução é possível encontrar a outra solução $x_{n-1}$ e com essas duas encontramos a próxima $x_{n-2}$ e assim sucessivamente. Como vimos, a solução para a $i_{essima}$ variável é

$\displaystyle{ x_{i} = \frac{a_{i,n+1} - \sum _{j=i+1}^{n}a_{ij}x_j}{a_{ii}} }$

Exemplo

Converta em traingular superior o seguinte sistema de equações utilizando o algoritmo acima: $\displaystyle{ \begin{array}{c c c c c c c} 3x_1 & - & 0.1x_2 & - & 0.2x_3 & = & 7.85\\ 0.1x_1 & + & 7.0x_2 & - & 0.3x_3 & = & -19.3\\ 0.3x_1 & - & 0.2x_2 & + & 10.0x_3 & = & 71.4 \end{array} }$
Expressando o pseudocódigos em python obtemos:

          def eliminacao (a, b):
            n = len(a)
            for k in xrange(0, n-1, 1):
              for i in xrange(k+1, n, 1):
                lambda = a[i*n+k] / a[k*n+k]
                for j in xrange(k, n, 1):
                  a[i*n+j] = a[i*n+j] - lambda * a[k*n+j]
                b[i] = b[i] - lambda * b[k]
            return(b, a)
        
[usuario@pclabfis: ]# python resolveGaussJordan.py principal.py
3.000000 -0.100000 -0.200000 7.850000
-0.000000 7.003334 -0.293333 -19.561666
0.000000 0.000000 10.012042 70.084290

Algoritmo de substituição regressiva

ENTRADA: matriz triangular de coeficientes $a_{ij}$ mais vetor $b_j$

$x_n = \frac{b_n}{a_{n,n}}$
Para $i=n-1,\; 1,\; -1$ faça
$soma = b_i$
Para $j=i+1,\; n$ faça
$soma = soma - a_{i,j}\cdot x_j$
fim do faça
$x_i = \frac{soma}{a_{i,i}}$
fim do faça

Os algoritmos anteriores permitem resolver alguns sistemas de equações lineares. Contudo, a grande maioria de sistemas não podem ser resolvidos por essas equações, um exemplo é o sistema a seguir

$\displaystyle{ \begin{array}{c c c c c c c c c c c} &-& x_2 &+& 6x_3 &-& x_4 &=& 5\\ -2x_1 &+& 2x_2 &+& 4x_3 &-& x_4 &=& 4\\ -x_1 &+& 3x_2 &-& 6x_3 &+& x_4 &=& -4\\ -4x_1 &+& 8x_2 &-& 4x_3 &-& x_4 &=& 0 \end{array} }$

Observe que o elemento $a_{1,1}$ (o pivô da linha $1$) é zero, isso implica que seria impossível calcular $\lambda = a_{i,k} /a_{k,k}$. Por outro lado, como os computadores lidam com poucos algarismos significativos em geral não é recomendável tentar resolver sistema em que o elemento pivô for muito pequeno (que este na ordem de grandeza da maquina - simples precisão $\sim 10^5$) pois isso pode incorrer em erros de arredondamento consideráveis (devido à substituição regressiva, que propaga o erro). A fim de resolver estes inconvenientes, não levados em consideração no algoritmo anterior, devemos utilizar a técnica do pivotamento parcial isto é, intercambiar a linha problemática com alguma outra linha que não apresente problema.

O algoritmo que vamos desenvolver a continuação deve ser utilizado na etapa de eliminação, logo após a escolha da $k_{esima}$ linha de pivô. A ideia é encontrar a linha $i \leq k$ que tenha o maior valor de $a_{i,k}$, encontrada essa linha deveremos trocar de posição a linha $k$ com essa linha $i$. Algo que deve ser verificado antes de efetuar a troca é que a linha em questão não seja menor que uma certa tolerância ($tol$) a qual consideraremos como sendo o verdadeiro zero numérico.

Algoritmo de Pivotamento Parcial

ENTRADA: $k$ linha pivô, matriz com coeficientes $a_{i,j}$ e $b_i$ onde $0 \leq i \leq n$ e $j \leq k$, $tol$, $erro$

imax = k
Para $i=k+1,\; n$ faça
Se $|a_{i,k}| > |a_{k,k}|$ então
imax = i
fim do Se
Se $|a_{imax,k}| < tol$ então
erro
PARA o programa
fim do Se
fim do faça
Se $imax \neq k$ então
Para $i=1,\; n$ faça
$tempA = a_{imax,i}$
$a_{imax,i} = a_{k,i}$
$a_{k,i} = tempA$
fim do faça
$tempB = b_{imax}$
$b_{imax} = b_k$
$b_k = tempB$
fim do Se

Os programas

Todos os pseudocódigos será explicitado em forma de programa a fim de você ver como se traduz do pseudocódigo para python. O principal ponto a ser notado é o referente a os índices, o python inicia em zero suas matrizes e os pseudocódigos iniciam em um, deve prestar atenção a esse fato.

Triangularização

      
      #!/usr/bin/env
      #-*- coding: utf-8 -*-

      import pivotamento

      def Superior(A, b):
        Ln = len(b)
        cn = Ln
        for Li in range(0, Ln-1, 1):
          A[:], b[:] = pivotamento.pivo(Li, b, A)
          for Lj in range(Li+1, Ln, 1):
            lam   =  A[Lj][Li] /A[Li][Li]
            b[Lj] = b[Lj] - b[Li] * lam
            for ci in range(Li, cn, 1):
              A[Lj][ci] =  A[Lj][ci] - A[Li][ci] * lam
        return (b, A)
      

Perceba que utilizamos as variáveis $L$ e $c$ para denotarmos as linhas e colunas, assim $Ln$ diz o número total de linhas entanto que $Li$ e $Lj$ dizem respeito às linhas $i$ e $j$ entanto que $ci$ diz respeito à coluna $i$.

Pivotamento:

      #!/usr/bin/env
      #-*- coding: utf-8 -*-

      import numpy as np

      def pivo(k, b, A):

        print 'Entrou: pivotando a linha', k
        print A
        n = len(b)
        jmax = k
        for j in xrange(k+1, n, 1):
          if ( np.abs(A[j][k]) > np.abs(A[jmax][k]) ):
            jmax = j

        TempA       = np.zeros_like(b)
        TempA[:]    = A[k][:]
        A[k][:]     = A[jmax][:]
        A[jmax][:]  = TempA[:]
        
        TempB  = b[k]
        b[k]   = b[jmax]
        b[jmax] = TempB
        
        print 'linha', k, 'trocou com', jmax
        print A, '\n'

        return A, b
      

Aqui tem um ponto importante a ser considerado, na linha 17 está escrito $TempA[:]$, se $[:]$ não for colocado o resultado será a criação de um link e não uma cópia com desejamos.

Retrossubstituição:

      #!/usr/bin/env
      #-*- coding: utf-8 -*-

      import numpy as np

      def retrosub(b, A):
        n = len(b)
        x = np.zeros_like(b)
        x[n-1] = b[n-1] / A[n-1][n-1]
        
        for i in xrange(n-2, -1, -1):
          soma = b[i]
          for j in xrange(i+1, n, 1):
            soma = soma - A[i][j]*x[j]
          x[i] = soma / A[i][i]
          
        return x
      

Neste programa deve se prestar atenção à forma como o índice $i$ da linha 11 se executa, por exemplo $xrange(5, 2, -1) = 5,\, 4,\, 3$, ou seja, incluí o inicio dos dados, por tanto, como as matrizes e vetores só chegam até $n-1$, no caso de series decrescentes devemos substituir o $n$ por $n-1$

Programa principal

      #!/usr/bin/env
      #-*- coding: utf-8 -*-

      import numpy as np
      import triangular
      import substituicao

      arq = open('matriz.dat', 'r')
      n = int(arq.readline().rstrip())

      A = np.zeros((n,n), dtype='float')
      b = np.zeros(n, dtype='float')

      for i in xrange(n):
        linha = arq.readline().rstrip().split()
        for j in xrange(n):
          A[i][j] = linha[j]


      linha = arq.readline().rstrip().split()
      for j in xrange(n):
        b[j] = linha[j]

      b, A = triangular.Superior(A, b)

      x = substituicao.retroativa(b, A)

      print '\nO resultado eh' 
      print 'x1 =', x[0]
      print 'x2 =', x[1]
      print 'x3 =', x[2]
      

Na versão do script principal, mostrada acima, a formatação dos dados de entrada se realiza, primeiro se coloca a dimensão da matriz que foi formada com os coeficientes das equações lineares, a continuação é colocada a matriz e a última linha é é colocado e vetor $b$ como sendo a última linha da nova matriz, por exemplo, o sistema

$\displaystyle{ \begin{array}{c c c c c c c c c} E_1: & x_1 & + & x_2 & & & + & 3x_4 & = & 4, \\ E_2: & 2x_1 & + & x_2 & - & x_3 & + & x_4 & = & 1, \\ E_3: & 3x_1 & - & x_2 & - & x_3 & + & 2x_4 & = & -3, \\ E_4: & -x_1 & + & 2x_2 & + & 3x_3 & - & x_4 & = & 4, \end{array} }$

é escrito como

\[ \begin{array}{rrrr} 4 & & & \\ 1 & 1 & 0 & 3 \\ 2 & 1 & -1 & 1 \\ 3 & -1 & -1 & 2 \\ -1 & 2 & 3 & -1 \\ 4 & 1 & -3 & 4 \\ \end{array} \nonumber \]

Mas há um jeito mais simples de se fazer o programa principal, utilizando loadtxt

      #!/usr/bin/env python
      #-*- coding: utf-8 -*-

      import numpy as np
      import trinagular
      from retrosubstituicao import *

      temp = np.loadtxt('matriz.dat')
      A = temp[:,0:-1]
      b = temp[:,-1]

      b, A = trinagular.Superior(A, b)
      x = retrosub(b, A)

      print '\nO resultado eh' 
      print 'x1 =', x[0]
      print 'x2 =', x[1]
      print 'x3 =', x[2]
      

Nesta abordagem os dados são dados na forma da matriz aumentada. Na linha 8 lemos toda a matriz aumentada; na linha 9 lemos todas as linhas (: significa todas) mas somente são lidas as colunas desde zero até a coluna $-1$ que significa a antepenúltima coluna (o $-1$ significa ultima coluna mas quando colocamos por exemplo $a[0:4]$ significa os elementos $a[0],\,a[1],\,a[2]$, ou seja, não incluí a fronteira). Dessa forma o arquivo com os dados é

\[ \begin{array}{rrrrr} 1 & 1 & 0 & 3 & 4 \\ 2 & 1 & -1 & 1 & 1 \\ 3 & -1 & -1 & 2 & -3 \\ -1 & 2 & 3 & -1 & 4\\ \end{array} \nonumber \]

executando o script (lembre que todos as funções devem estar no mesmo arquivo da função principal):

[usuario@pclabfis: ]# python principal.py
Entrou: pivotando a linha 0
[[ 1. 1. 0. 3.]
 [ 2. 1. -1. 1.]
 [ 3. -1. -1. 2.]
 [-1. 2. 3. -1.]]
linha 0 trocou com 2
[[ 3. -1. -1. 2.]
 [ 2. 1. -1. 1.]
 [ 1. 1. 0. 3.]
 [-1. 2. 3. -1.]]

Entrou: pivotando a linha 1
[[ 3. -1. -1. 2. ]
 [ 0. 1.66666667 -0.33333333 -0.33333333]
 [ 0. 1.33333333 0.33333333 2.33333333]
 [ 0. 1.66666667 2.66666667 -0.33333333]]
linha 1 trocou com 3
[[ 3. -1. -1. 2. ]
 [ 0. 1.66666667 2.66666667 -0.33333333]
 [ 0. 1.33333333 0.33333333 2.33333333]
 [ 0. 1.66666667 -0.33333333 -0.33333333]]

Entrou: pivotando a linha 2
[[ 3.00000000e+00 -1.00000000e+00 -1.00000000e+00 2.00000000e+00]
 [ 0.00000000e+00 1.66666667e+00 2.66666667e+00 -3.33333333e-01]
 [ 0.00000000e+00 0.00000000e+00 -1.80000000e+00 2.60000000e+00]
 [ 0.00000000e+00 0.00000000e+00 -3.00000000e+00 5.55111512e-17]]
linha 2 trocou com 3
[[ 3.00000000e+00 -1.00000000e+00 -1.00000000e+00 2.00000000e+00]
 [ 0.00000000e+00 1.66666667e+00 2.66666667e+00 -3.33333333e-01]
 [ 0.00000000e+00 0.00000000e+00 -3.00000000e+00 5.55111512e-17]
 [ 0.00000000e+00 0.00000000e+00 -1.80000000e+00 2.60000000e+00]]

O resultado eh
x1 = 1.0
x2 = 1.66533453694e-16
x3 = 2.0

Função intrínseca do numpy

O objetivo de criar nossa própria função é poder entender e interpretar algoritmos que frequentemente serão visto ao longo de nossa careira como físico, isso permitirá que seu futuro você veja uma equação matemática e consiga traduzir a alguma linguagem sem maiores problemas (passando a se dedicar assim à comprensã da estrutura matemática da equação e não em como colocar no computador essa equação). Obviamente o numpy possui várias funçõe que são própria para tratar com sistemas lineares, entre elas uma função que resolve o sistema linear:

      #!/usr/bin/env python
      #-*- coding: utf-8 -*-

      import numpy as np

      temp = np.loadtxt('matriz.dat')
      A = temp[:,0:-1]
      b = temp[:,-1]

      x = np.linalg.solve(A,b)

      print '\nO resultado eh' 
      print 'x1 =', x[0]
      print 'x2 =', x[1]
      print 'x3 =', x[2]
      
[usuario@pclabfis: ]# python sistLinear.py
O resultado eh
x1 = 1.0
x2 = 0.0
x3 = 2.0

Tarefa

Implemente o algoritmo de pivotamento parcial e resolva os seguintes sistema de equações lineares