A idéia dos métodos iterativos é produzir uma sequência de estimativas vetoriais $ \mathbf{x}^{(k)} $, que idealmente convirja para a solução $\mathbf{x}$ e cada passo é calculado a partir da estimativa anterior. Costumam ser bem mais rápidos que oe métodos diretos no caso de matrizes muito grandes. Para construir recursivamente a sequência que dá origem ao método iterativo fazemos uma decomposição da matriz $A$ de coeficientes do sistema como $$ A = N - P$$ Onde $N$ é de uma classe de sistemas inversíveis e fáceis de resolver (ou seja, fáceis de inverter). Então resolver o sistema: $$ A\mathbf{x} =b $$ é equivalente a resolver a equação: $$ N\mathbf{x} = b+ P\mathbf{x} \text{ ou } \mathbf{x}=N^{-1}(b+ P\mathbf{x})$$ A recursão $$ N\mathbf{x}^{(k+1)} = b+ P\mathbf{x}^{(k)}$$ é a que converge para a solução caso $ ||N^{-1}P|| < 1$ onde a norma de matrizes é induzida de uma norma em $\mathbb{R}^n$ (veja no apêndice). Destacamos dois casos
Método de Jacobi: $ N = \text{diag}(a_{11}, \cdots, a_{nn})$
Método de Gauss-Seidel : $ N = (n_{ij}) \text{ com } n_{ij} = a_{ij} \text{ se } i\geq j \\ n_{ij}=0 \text{ se } i\lt j$
Vamos ver com mais detalhes só o método de Gauss-Seidel. Vamos denotar por $x_i^{(k)}$ a $i$-ésima coordenada do vetor estimativa, no passo $k$. Então podemos escrever: $$ \begin{eqnarray} x_1^{(k+1)} & = &\frac{1}{a_{11}}\left( b_1 - \sum_{j=2}^n a_{1j}x_j^{(k)}\right) \\ x_2^{(k+1)} & = &\frac{1}{a_{22}}\left( b_2 - a_{21}x_1^{(k+1)} -\sum_{j=3}^n a_{2j}x_j^{(k)}\right) \\ & \vdots & \\ x_l^{(k+1)} & = & \frac{1}{a_{ll}}\left( b_l - \sum_{j=1}^{l-1}a_{lj}x_j^{(k+1)} -\sum_{j=l+1}^n a_{lj}x_j^{(k)}\right)\\ & \vdots & \\ x_n^{(k+1)} & = & \frac{1}{a_{nn}}\left( b_n - \sum_{j=1}^{n-1}a_{lj}x_j^{(k+1)}\right) \end{eqnarray}$$
from pylab import *
# dou a matriz quadrada de coeficientes e vetor independente
A=matrix('4 2 0.6; 2 5.5 3; 1 0 2.5')
b=matrix('1. ; 3.; 1.')
def iterGS(x, A, b):
xold=x.copy()
L=list(x)
for i in range(len(A)):
soma = 0.
for j in range(len(A)):
if (j != i): soma += -A[i,j]*L[j]
L[i] = (b[i] + soma)/A[i,i]
return L
def GaussSeidel(x0,A,b,eps,n):
IterGauss=0
dx=10.
while ((IterGauss<=n)&(dx>eps)):
xold = x0.copy()
L=iterGS(x0, A, b)
for k in range(len(L)) : L[k]=float(L[k])
x0=matrix(L).T
IterGauss += 1
dx = norm(x0-xold)
return x0 , IterGauss, dx
########## Teste com um exemplo ##############
x0=matrix(' 0; 0; 0 ')
sol = GaussSeidel(x0,A,b,0.0000001, 100)
print (sol)
(matrix([[ 0.03039834], [ 0.32285115], [ 0.38784066]]), 16, 3.5604346043801354e-08)
%pastebin