2.2.2 Sistemas Lineares — Métodos Iterativos — Jacobi
1. O método de Jacobi
Seja o sistema linear tal que a matriz dos coeficientes possa ser decomposta na seguinte forma:
onde é uma matriz diagonal e é a matriz que possui diagonal nula. Desta forma, o sistema é decomposto como
$$
Dx + Rx = b
$$
onde podemos isolar o vetor formado pela diagonal .
O passo iterativo do método consiste em sucessivas iterações para aproximar a solução através da fórmula:
$$x^{(k+1)} = D^{-1} (b – Rx^{(k)}) $$
ou em uma forma mais simples
$$ x_{i}^{(k+1)} = \dfrac{\left( b_i – \sum_{j \neq i} a_{ij}x_{j}^{(k)} \right)}{a_{ii} $$
2. SOR Adaptativo
Como apresentamos no post sobre o método de Gauss-Siedel, podemos inserir uma variável na formula iterativa para controlar o comportamento da convergência. Tal recurso é chamado de SOR (Sucessive OverRelaxation).
Uma propriedade interessante da variável de relaxação é que ela permite transformar um processo não-convergente em outro convergente ou acelerar a taxa de convergência. Quando escolhemos um valor entre , observamos o primeiro caso, enquanto que valores entre favorecem o segundo. Assim, podemos dizer que, respeitando estes intervalos, quanto maior o valor para variável de relaxação, maior será a velocidade de convergência, caso o sistema seja predominantemente diagonal.
Na nossa implementação, utilizamos um valor dinâmico para variável de relaxação. Utilizamos um valor inicial que favorece a convergência (entre 0 e 1). Caso o sistema se mantenha com uma taxa de convergência constante entre duas iterações consecutivas, aumentamos o valor da variável com o objetivo de aumentar esta taxa. Se nesta configuração detectarmos divergência, diminuímos a variável de relaxação para forçar o solução do sistema.
3. Implementação
def jacobi(A, b, errto=10E-10, maxiter=10E5, SOR=1.5, getiters=False):
"""
Return a list with the solution of linear system.
jacobi(A, b)
INPUT:
* A: coefficients matrix
* b: a list, 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/
Oct 2010
"""
A = [map(float, i) for i in A]
n = len(A)
errnoold = 0.0
errno = 0.0
iters = 0
sentry = False
xold = [1.0 / float(n + errto) for i in xrange(n)]
xnew = [1.0 / float(n + errto) 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] * xold[i]
xold[i] = summ
while not sentry and iters < = maxiter:
sentry = True
for i in xrange(n):
summ = b[i]
for j in xrange(n):
if i != j:
summ -= A[i][j] * xold[j]
xnew[i] = SOR * summ + (1.0 - SOR) * xold[i]
if sentry and xnew[i] != 0:
errno = abs((xnew[i] - xold[i]) / xnew[i])
if errno > errnoold:
newSOR = SOR - 0.1
if newSOR >= errto:
SOR = newSOR
else:
newSOR = SOR + 0.1
if newSOR < = 1.9:
SOR = newSOR
if errno > errto:
sentry = False
xold[:] = xnew[:]
errnoold = errno
iters += 1
if getiters:
return (xnew, iters)
else:
return xnew
Esta função está disponível em http://www.sawp.com.br/code/linear/jacobi.py, assim como um exemplo de como utilizá-la.
4. 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/.