2.2.2 Sistemas Lineares — Métodos Iterativos — Jacobi

1. O método de Jacobi

Seja o sistema linear \(Ax = b \) tal que a matriz dos coeficientes \(A \) possa ser decomposta na seguinte forma:

\(A = D + R \)

onde \(D \) é uma matriz diagonal e \(R \) é a matriz que possui diagonal nula. Desta forma, o sistema é decomposto como

\(
Dx + Rx = b
\)

onde podemos isolar o vetor formado pela diagonal \(Dx = b – Rx \) .

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 \(]0,1[ \) , observamos o primeiro caso, enquanto que valores entre \(]1,2[ \) 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/.

References

[1] Anthony Ralston and Philip Rabinowitz, A First Course in Numerical Analysis (2nd ed.), McGraw-Hill and Dover, (2001), 423-424.

[2] N.B Franco, Cálculo Numérico, Pearson Prentice Hall, (2006), 69.