2.2.1 Sistemas Lineares — Métodos Iterativos — Gauss-Siedel
1. O método de Gauss-Siedel
A resolução de sistemas lineares através de métodos iterativos foram desenvolvidos como uma forma específica de busca de raízes a partir de derivações dos métodos que buscam raízes em funções não-lineares. Como em todos métodos numéricos iterativos, a abordagem para solução consiste em escolher um valor inicial (estimado) para convergir à uma solução a partir deste. Como sistemas lineares utilizam um conjunto de equações tais que as raízes devam satisfazê-las simultaneamente, os métodos iterativos são úteis neste contexto.
Como trata-se de um sistema linear de equações, o método de Gauss-Siedel resolve problemas na forma
Primeiramente, suponhamos um sistema de três equações, gerando uma matriz de coeficientes . Se a diagonal é formada apenas por elementos não-nulos, podemos isolar na primeira equação, na segunda e na terceira. Assim, obteremos
x_2 = \dfrac{b_2 – a_{21}x_1 – a_{23}x_3}{a_{22},
x_3 = \dfrac{b_3 – a_{31}x_1 – a_{32}x_2}{a_{33} $$
Pelas equações acima notamos que para estimar , inicialmente precisamos de e . A ideia do método consiste em escolher um valor inicial para estas variáveis e, em seguida, fazer sucessivas substituições até que os valores de estejam se aproximando de uma solução. Por esta razão, o método de Gauss-Siedel é conhecido como método das substituições sucessivas.
Para generalizar a ideia do parágrafo anterior, denotaremos o método através de duas variáveis e , onde é a componente do vetor de incógnitas e é a iteração. Desta maneira, uma variável qualquer é expressa como . Considerando o sistema de equações e incógnitas, começamos por calcular a primeira variável
\left( \sum_{j=2}^{n} a_{1j}x_{i}^{(j)} -b_1 \right) $$
Para próxima componente , utilizamos a segunda equação, já substituindo pelo valor calculado anteriormente
\left( a_{21}x_{i+1}^{(1)} +
\sum_{j=3}^{n} a_{2j}x_{i}^{(j)} – b_2 \right) $$
a partir destas duas últimas equações, podemos facilmente observar a relação de recorrência, a qual expressamos como a fórmula geral do método de Gauss-Siedel:
\left( \sum_{j=1}^{r-1} a_{rj}x_{i+1}^{(j)} +
\sum_{j=r+1}^{n} a_{rj}x_{i}^{(j)} – b_r \right) $$
onde .
2. Convergência do Método de Gauss-Siedel
Assim como todo método de busca de raízes iterativo, o algoritmo de Gauss-Siedel é uma técnica baseada no método do ponto fixo, fazendo com que o método de Gauss-Siedel seja suscetível à duas limitações fundamentais:
- casos onde não há convergência
- casos onde a convergência é muito lenta
Para ambas situações, observamos que os valores do sistema que condicionam os casos problemáticos. Para verificar quais casos são estes, deduzimos os critérios de convergência a partir de duas equações não-lineares e tal que
\left| \dfrac{\partial u}{\partial x_1} \right| +
\left| \dfrac{\partial u}{\partial x_2} \right| < 1 [/latex]
e
</center>
\left| \dfrac{\partial v}{\partial x_1} \right| +
\left| \dfrac{\partial v}{\partial x_2} \right| < 1 [/latex]
Supondo [latex]u v $$ como sendo duas funções lineares pertencente à um sistema com apenas duas variáveis, temos que tais funções são expressas como
</center>
u(x_1, x_2) = \frac{b_1}{a_{11} – \frac{a_{12}{a_{11} x_2
$$
e
v(x_1, x_2) = \frac{b_2}{a_{22} – \frac{a_{21}{a_{22} x_1
$$
Para sistemas lineares, as derivadas parciais eliminam todas as componentes que não contém a variável a qual estamos diferenciando, ou seja
$$
\frac{\partial u}{\partial x_1} = 0 \text{~,~}
\frac{\partial u}{\partial x_2} = – \frac{a_{12}{a_{11}
$$
e
$$
\frac{\partial v}{\partial x_1} = – \frac{a_{21}{a_{22}
\text{~,~} \frac{\partial v}{\partial x_2} = 0
$$
Utilizando os critérios das funções e , podemos reescrever as condições como
$$
\left| \frac{a_{12}{a_{11} \right| < 1 \text{~e~} \left| \frac{a_{21}{a_{22} \right| < 1 [/latex]
Isso significa que para haver convergência, o sistema precisa ter os maiores coeficientes concentrados na diagonal. Ou seja,
</center>
[latex]
|a_{ii}| = \sum_{j=1 and j \ne i}^{n} |a_{ij}|
$$
3. Aceleração de Processos Iterativos Estacionários
A taxa de convergência de um processo iterativo estacionário depende do raio espectral do vetor de elementos independentes. Assim, qualquer qualquer modificação no vetor implica em uma modificação na taxa de convergência. Isto é, a redução do raio espectral de irá aumentar a taxa de convergência. Considerando que queremos acelerar a convergência às raízes no processo iterativo, adotaremos um método conhecido como relaxação sucessiva.
Primeiramente, rearranjamos o sistema linear para uma forma onde todos elementos da diagonal são unitários, ou seja, . Isso pode ser feito dividindo-se toda a linha pelo elemento da diagonal pertencente à ela.
Pelo método de Gauss-Siedel, sabemos que
\sum_{j=r+1}^{n} a_{rj} x_{i}^{(j)} + b_r
= x_{i}^{(r)}
– \sum_{j=1}^{r-1} a_{rj} x_{i+1}^{(j)} –
\sum_{j=r+1}^{n} a_{rj} x_{i}^{(j)} + b_r $$
para . Assim, podemos melhorar a formula acima reescrevendo-a como
x_{i+1}^{(r)} = x_{i}^{(r)} – \omega \left(
\sum_{j=1}^{r-1} a_{rj} x_{i+1}^{(j)} –
\sum_{j=r+1}^{n} a_{rj} x_{i}^{(j)} + b_r \right)
$$
onde é o fator de relaxação do processo iterativo.
Para facilitar a visualização deste processo, imagine que a matriz de coeficientes é decomposta e escrita na forma
onde é uma matriz diagonal, é triangular inferior e é triangular superior. Assim, podemos escrever o método de Gauss-Siedel como
onde e .
Tal processo pode ser interpretado como uma seleção baseada em uma média ponderada dos resultados da iteração anterior e da iteração corrente:
x_{i+1}^{(j)} = \omega x_{i+1}^{(j)} + (1 – \omega) x_{i}^{(j)}
$$
onde deve ser escolhido como um valor entre e .
3.1. Interpretações para </p>
Para , temos que e o resultado não é modificado. Contudo, se estiver no intervalo entre e , o resultado é uma média ponderada do resultado atual e do anterior. Neste último caso, temos uma sub-relaxação, que é usada para forçar o sistema a convergir, mesmo que ele seja não-convergente.
Para no intervalo entre e , supomos que o condicionamento do sistema é certamente convergente, mas converge muito lentamente. Neste caso, atua como um fator de peso que favorece a convergência entre duas iterações. Assim, dizemos que esta abordagem é uma super-relaxação, também conhecida como sucessive overrelaxation — SOR.
3.2. Análise de SOR </p>
A escolha de requer um certo cuidado, uma vez que este valor é fortemente ligado aos dados do sistema, sendo dependente do problema e determinado de maneira empírica. Para exemplificar a utilização de um fator de relaxação, a seguir apresentamos o comportamento do método de Gauss-Siedel para um sistema em relação à diferentes valores de .
O primeiro procedimento para análise consiste em resolver um sistema linear na forma diversas vezes, variando apenas a variável de relaxação no intervalo com um passo de e verificando o comportamento das iterações necessárias para resolução do problema. Tal procedimento nos permite gerar o gráfico da Figura sor1.
Mostra o número de iterações necessárias para resolver o mesmo sistema linear sem variação na precisão do resultado
Como podemos observar da Figura sor1, o ponto ótimo para resolução desse sistema consiste em utilizar um valor para o fator de relaxação próximo ao ponto . Sendo assim, repetimos o experimento para o subintervalo e observamos o comportamento ilustrado na
Figura sor2.
Mostra o número de iterações necessárias para resolver o mesmo sistema linear sem variação na precisão do resultado no intervalo próximo ao mínimo global
Desta última figura podemos concluir que o valor que minimiza o tempo de resolução do sistema linear consiste em um , sendo notável a diferença no desempenho proporcionada por por este valor: menos de iterações ante as iterações do pior caso, significando uma melhoria de quase .
A Figura sor2 também nos fornece informações para avaliação da necessidade de escolha empírica do valor do parâmetro de relaxação. Se notarmos o comportamento do gráfico após o valor ótimo, notamos que o crescimento é contínuo, mas com uma pequena oscilação. Este tipo de comportamento gera diversos pontos de mínimo, o que torna o problema de encontrar o valor ótimo um problema difícil de se automatizar com eficiência. Geralmente, problemas como estes são resolvidos com heurísticas — tais como algoritmos genéticos — e costumam ser muito mais caros de se resolver do que problemas de resolução de sistemas lineares, o que torna injustificável a utilização da busca pelo valor ótimo do fator de relaxação.
Portanto, para definir um adequado para algum problema, é preferível que façamos alguns testes com instâncias menores do problema geral para conseguirmos obter um valor próximo ao ideal.
4. Implementação
def gauss_siedel(A, b, errto=10E-10, maxiter=10E5, SOR=1.77, getiters=False):
"""
Return a list with the solution of linear system.
gauss_siedel(A, b)
INPUT:
* A: coefficients matrix
* b: a list, relative a vector of independent terms
return:
if getiters=False:
a list with unknown vector 'x' of system A*x = b.
else:
the tuple (list x, integer iters)
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/
Sep 2010
"""
n = len(A)
errno = 0.0
iters = 0
sentry = False
x = [1.0 / float(n) for i in xrange(n)]
for i in xrange(n):
diagonal = float(A[i][i])
for j in xrange(n):
A[i][j] /= diagonal
b[i] /= diagonal
for i in xrange(n):
summ = b[i]
for j in xrange(n):
if i != j:
summ -= A[i][j] * x[i]
x[i] = summ
while not sentry and iters < = maxiter:
sentry = True
for i in xrange(n):
xold = x[i]
summ = b[i]
for j in xrange(n):
if i != j:
summ -= A[i][j] * x[j]
x[i] = SOR * summ + (1.0 - SOR) * xold
if sentry and x[i] != 0:
errno = abs((x[i] - xold) / x[i])
if errno > errto:
sentry = False
iters += 1
if getiters:
return (x, iters)
else:
return x
Nesta implementação fazemos duas considerações importantes para a melhoria do desempenho:
- a primeira consiste em dividir cada linha da matriz de coeficientes e dos termos independentes pelo elemento respectivo da diagonal. Conforme discutimos na seção de análise de convergência, isto reduz o número total de operações no algoritmo.
- a segunda consideração está na definição da variável de controle sentry. Esta variável verifica o erro da aproximação e controla a execução da função. Somente se alguma das equações possuir um erro de aproximação maior do que o erro tolerado é que a variável de controle permite a execução. O uso da variável sentry nos permite economizar alguns laços de repetição e evitar cálculos desnecessários. </p>
Esta função está disponível em http://www.sawp.com.br/code/linear/gausssiedel.py, assim como um exemplo de como utilizá-la.
5. Copyright
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/.