2.1.1 Sistemas Lineares — Métodos Diretos — Eliminação de Gauss
1. Introdução
Em um sistema linear formado por um par de equações, o processo para obtenção se descobrir as duas incógnitas consiste em isolar uma delas e substituir na outra equação, eliminando assim, uma das incógnitas.
Para um sistema maior, podemos estender a abordagem para um número maior de equações através do algoritmo conhecido como “Eliminação de Gauss”, que basicamente manipula algumas variáveis, a fim de gerar uma eliminação ao longo da execução, ao tempo que vai substituindo regressivamente nas equações originais, a fim de obter a solução final.
Sendo assim, o sistema geral que iremos estudar é formado por um conjunto arbitrário de equações tais que:
$$ a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n = b_1 $$
$$ a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n = b_2 $$
$$ a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n = b_n $$
Entre os métodos computacionais temos que o processo de diagonalização de matrizes é uma técnica-chave para solução de sistemas lineares. Sendo assim, a primeira fase da eliminação de Gauss está em reduzir a matriz dos coeficientes para uma matriz triangular-superior, eliminando da primeira variável em todas as equações diferentes de .
Para isso, vamos multiplicar a Equação 1 por . Com isso, obtemos
$$ a_{21}x_1 + \dfrac{a_{21}{a_{11}a_{12}x_2 + \cdots + \dfrac{a_{21}{a_{11}a_{1n}x_n = \dfrac{a_{21}{a_{11}b_{1} $$
Agora, subtraímos a Equação 4 da Equação 2, o que nos permite obter
$$ \left(a_{22} – \dfrac{a_{21}{a_{11}a_{12}\right)x_2 + \cdots + \left(a_{2n} – \dfrac{a_{21}{a_{11}a_{1n}\right)x_n = b_2 – \dfrac{a_{21}{a_{11}b_1 $$
substituindo cada termo dentro dos parênteses por um novo , onde , temos que a Equação 5 pode ser escrita como
$$ a’_{22}x_2 + \cdots + a’_{2n}x_n = b’_2 $$
onde .
Assim como fizemos com a segunda linha do sistema, repetimos o processo para as demais equações. Ou seja, multiplicamos pela Equação 1 e subtraímos o resultado da equação 3, e continuamos o processo até a última equação.
No processo de eliminação, o termo é conhecido como “pivô”, e a multiplicação de cada linha pelo elemento é denominada “normalização”. Portanto, podemos notar que o pivô nulo é um problema no processo de normalização, já que provoca uma divisão por zero. Por isso, em implementações reais, devemos utilizar algumas técnicas de “re-normalização” para eliminarmos o pivô nulo do processo.
Após a primeira eliminação, a qual estamos descrevendo até aqui, o sistema linear modificado ficará:
$$ a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n = b_1 $$
$$ \hspace{35pt} a’_{22}x_2 + \cdots + a’_{2n}x_n = b’_2 $$
$$ \hspace{35pt} a’_{n2}x_2 + \cdots + a’_{nn}x_n = b’_n $$
agora, repetimos o mesmo processo de eliminação, mas tomando como coeficiente de normalização como sendo , aplicando-o à todas equações a partir da terceira. Este processo deve continuar usando os pivôs obtidos nas equações remanescentes. Ou seja, repetimos o pivotamento e normalização descritas para até a equação , eliminando até o termo da equação. Com isso, obteremos o sistema triangular a seguir:
$$ a_{11}x_1 + a_{12}x_2 + a_{13}x_3 + \cdots + a_{1n}x_n = b_1 $$
$$ \hspace{35pt} a’_{22}x_2 + a’_{23}x_3 + \cdots + a’_{2n}x_n = b’_2 $$
$$ \hspace{70pt} a”_{33}x_3 + \cdots + a”_{3n}x_n = b”_3 $$
$$ \hspace{110pt} a^{(n-1)}_{nn}x_n = b^{(n-1)}_n $$
A partir do sistema triangular, a solução é facilmente obtida, a começar por :
$$ x_n = \dfrac{b^{n-1}_n}{a^{n-1}_{nn} $$
substituindo esse resultado na equação anterior, isto é, na equação , cuja dependência está apenas nas variáveis e , podemos determinar também . Esse processo de ir substituindo as incógnitas já encontradas nas equações anteriores chamaremos de backtrack. A partir dele, encontramos todos pela da fórmula:
$$ x_i = \dfrac{ b^{(i-1)}_{i} – \sum^{n}_{j=i+1} a^{(i-1)}_{ij}x_j }{a^{i-1}_{ii} $$
para
2. Implementação </p>
A implementação da forma mais simples do processo de eliminação progressiva, seguida da substituição regressiva segue abaixo.
def gauss(A, b):
"""
Return a list with the solution of linear system.
gauss(A, b)
INPUT:
* A: a matrix of coefficients
* b: a list, relative a vector of independent terms
return: a list with all unknown variables of system.
Author: Pedro Garcia [sawp@sawp.com.br]
see: http://www.sawp.com.br
License: Creative Commons
http://creativecommons.org/licenses/by-nc-nd/2.5/br/
Mar 2010
"""
def progressive_elimination(A, b):
n = len(b)
for k in xrange(n-1):
for i in xrange(k+1, n):
factor = A[i][k]/A[k][k]
for j in xrange(k+1, n):
A[i][j] = A[i][j] - factor * A[k][j]
b[i] = b[i] - factor * b[k]
return (A, b)
def regressive_replace(A, b):
n = len(b)
x = [0.0 for i in xrange(n)]
x[n-1] = b[n-1]/A[n-1][n-1]
for i in xrange(n-1, -1, -1):
suma = b[i]
for j in xrange(i+1, n):
suma = suma - A[i][j] * x[j]
x[i] = suma/A[i][i]
return x
(A, b) = progressive_elimination(A, b)
return regressive_replace(A, b)
Embora este código funciona bem para uma série de equações ele produzirá um erro cada vez que o pivô for um elemento nulo, pois isso implicará em divisão por zero.
Como os pivôs são formados apenas pelos elementos da diagonal principal, antes do programa executar a eliminação progressiva, ele pode verificar a existência de nulidade dos pivôs. Caso haja, basta que a linha nula seja trocada por outra, já que isso não afeta o resultado final. Esta função pode ser, por exemplo:
def remove_null_pivot(A, b):
n = len(b)
checked = 0
i = 0
while i < n and checked < n:
if A[i][i] == 0.0:
if i + 1 == n-1:
A[i] = A[i-1]
b[i] = b[i-1]
else:
A[i] = A[i+1]
b[i] = b[i+1]
i = 0
checked += 1
return (A, b)
A implementação destas funções, bem como um exemplo de como utilizá-la está disponível em http://www.sawp.com.br/code/linear/gauss.py.
3. Refinando a Implementação </p>
Ainda que muitos sistemas de equações possam ser resolvidos com a última implementação apresentada, há algumas peculiaridades de casos em que ela não é eficiente. Para ser usado em um programa robusto, devemos fazer com que o algoritmo preveja as seguintes situações de erro:
- Divisão por zero: quando o pivô é nulo, o programa mostrado anteriormente acabará gerando uma divisão por zero, o que implicará em uma possível falha de todo o sistema. Para evitar este tipo de problema, um tratamento sobre os dados deve ser feito. Esse tratamento consiste em uma técnica chamada pivotamento.
-
Roundoff Error: “Erros de arrendondamento.” Um problema muito comum no tratamento de números reais por computadores. O problema de arrendondamento é um fator inerente às máquinas reais. Isso ocorre porque os computadores possuem um número finito de algarismos significativos, implicando em limites na expressão de números reais através da notação de ponto flutuante. O problema de arrendondamento pode ser propagando em um grande número de iterações ou, como é o caso da Eliminação de Gauss, que é um método de solução direta, na manipulação de um grande número de equações. Como erros de arrendondamento induzem mudanças nos coeficientes dessas equações, essas mudanças podem levar a erros grandes nas soluções de sistemas sem balanceamento.
Um fator importante para evitar a propagação de erros, a fim de se obter um resultado com certa precisão desejada, consiste em utilizar técnicas de arrendondamento. Ralson e Rabinowitz[rabinowitz] fazem uma excelente discussão formal de tratamento de erros de arrendondamento, que vale a pena ser consultada.
- Mal condicionamento: Sistemas “mal condicionados”, como o próprio nome sugere, são aqueles que por sua própria formação estão sob condições ruins. Isto é, eles estão desbalanceados, de forma que a realização de operações entre seus coeficientes geram resultados ruins, normalmente associados à propagação de erros de arrendondamento. A primeira implicação de problemas mal condicionados é que pequenas mudanças nos seus coeficientes implicam em grande divergência na solução correta. Assim, haveria uma quantidade grande de respostas que serviriam como aproximação da solução do sistema, o que é uma coisa indesejável, uma vez que sistemas lineares possuem solução única. </p>
-
Sendo assim, com o intuito de contornar esses erros, intrínsecos à precisão computacional, as técnicas para tornar o nosso programa mais robusto são:
- Aumentar o número de algarismos significativos: Com mais casas disponíveis para manipulação dos dados, os problemas do mal condicionamento e da propagação do erro são minimizados ao longo de uma série muito grande de cálculos. Por isso, o ideal é sempre usarmos a maior palavra disponível para manipulação de ponto flutuante disponível para nosso sistema (em geral, sempre que precisarmos de valores de até casas de precisão (equivalente ao tipo float em C), devemos prefirir usar o tipo double. As arquiteturas mais atuais possuem tecnologias avançadas de tratamento de ponto flutuante, tais como SSE, que contém registradores de até bits, o que permite uma precisão muito superior, e que deve ser explorada).
-
Pivoteamento: Para solucionar o problema do pivô nulo (ou muito menor que o elemento da normalização), podemos balancear o sistema, a fim de gerar uma pré-normalização, a fim de evitar uma disparidade muito grande entre os operandos. Assim, antes que cada linha seja normalizada, o programa deve determinar o maior coeficiente disponível na coluna logo abaixo do pivô. Em seguida, as linhas devem ser trocadas para que o maior elemento dado seja o pivô. Após este re-pivotamento parcial, buscamos o maior elemento nas colunas, a fim de encontrar também o maior elemento dentre as colunas, refinando assim o melhor pivô a ser usado.
O pivotamento, além de resolver o problema do pivô nulo, também auxilia na diminuição da propagação do erro de arrendondamento, sendo uma solução eficiente e altamente recomendada para implementações reais.
Um pivotamento pode ser realizado pela seguinte função:
def pivot(A, b, k): n = len(b) p = k bigger = abs(A[k][k]) for i in xrange(k+1, n): sentinel = abs(A[i][k]) if (sentinel > bigger): bigger = sentinel p = i if p != k: for i in xrange(k, n): (A[p][i], A[k][i]) = (A[k][i], A[p][i]) (b[p], b[k]) = (b[k], b[p])Este função é projetada para ser chamada antes de cada nova normalização do sistema. Isso é, como cada eliminação progressiva é feita em cada linha, pivôs serão gerados, como já discutimos anteriormente. Sendo assim, para cada um deles, devemos chamar a função pivot para permitir o balanceamento do sistema.
Da função pivot, temos que a linha captura inicialmente a posição do primeiro pivô, enquanto que a linha coloca na variável bigger, o seu valor. O laço de repetição definido da linha à é que realiza a busca na matriz dos coeficientes. Caso haja um valor maior que o pivô inicial, seu valor e posição são atualizados nas variáveis bigger e p. As linhas que seguem servem para verificar se a posição do maior pivô foi alterado. Caso isso tenha ocorrido, o sistema é rearranjado para a forma pivotada, conforme descrevemos anteriormente.
-
4. Implementação da Eliminação de Gauss </p>
A implementação completa da eliminação de Gauss, considerando o pivotamento segue abaixo:
def gauss_robust(A, b):
"""
Return a list with the solution of linear system.
gauss_robust(A, b)
INPUT:
* A: a matrix of coefficients
* b: a list, relative a vector of independent terms
return: a list with all unknown variables of system.
Author: Pedro Garcia [sawp@sawp.com.br]
see: http://www.sawp.com.br
License: Creative Commons
http://creativecommons.org/licenses/by-nc-nd/2.5/br/
Mar 2010
"""
def pivot(A, b, k):
n = len(b)
p = k
bigger = abs(A[k][k])
for i in xrange(k+1, n):
sentinel = abs(A[i][k])
if (sentinel > bigger):
bigger = sentinel
p = i
if p != k:
for i in xrange(k, n):
(A[p][i], A[k][i]) = (A[k][i], A[p][i])
(b[p], b[k]) = (b[k], b[p])
return (A, b)
def progressive_elimination(A, b):
n = len(b)
for k in xrange(n-1):
(A, b) = pivot(A, b, k)
for i in xrange(k+1, n):
factor = A[i][k]/A[k][k]
for j in xrange(k+1, n):
A[i][j] = A[i][j] - factor * A[k][j]
b[i] = b[i] - factor * b[k]
return (A, b)
def regressive_replace(A, b):
n = len(b)
x = [0.0 for i in xrange(n)]
x[n-1] = b[n-1]/A[n-1][n-1]
for i in xrange(n-1, -1, -1):
suma = b[i]
for j in xrange(i+1, n):
suma = suma - A[i][j] * x[j]
x[i] = suma/A[i][i]
return x
(A, b) = progressive_elimination(A, b)
return regressive_replace(A, b)
Note que pouca alteração foi necessária em relação à forma mais ingênua do algoritmo. Este código está disponível em http://www.sawp.com.br/code/linear/gauss_robust.py, bem como um exemplo de como utilizá-lo.
5. Copyright </p>
Este documento é disponível sob a licença Creative Commons. As regras dos direitos de cópia deste conteúdo estão acessíveis em http://creativecommons.org/licenses/by-nc-nd/2.5/br/.