4.2 Gauss-Seidel – Algoritmo

El método se describe con un sistema de 3 incógnitas y 3 ecuaciones.
Note el uso de índices ajustado a la descripción de python para la primera fila=0 y primera columna=0:

a00x0 + a01x1+ a02x2 = b0
a10x0 + a11x1+ a12x2 = b1
a20x0 + a21x1+ a22x2 = b2

La forma matricial del sistema original de ecuaciones es: [A]{X} = {B}

Se usa los indices i para fila y j para columna, asi un elemento de la matriz es aij .
Para obtener los valores de cada xi, se despeja uno por cada por cada fila i:

x0 = (b0        - a01x1 - a02x2) /a00
x1 = (b1 - a10x0        - a12x2) /a11
x2 = (b2 - a20x0 - a21x1         ) /a22

Observe el patrón de cada fórmula usada para determinar x[i]:

x_i = \bigg(b_i - \sum_{j=0, j\neq i}^n A_{i,j}X_j\bigg) \frac{1}{A_{ii}}
la sumatoria para cada término de A[i,j] en la fila i, 
excepto para el término de la diagonal A[i,i]

Si se tiene conocimiento del problema planteado y se puede «intuir o suponer» una solución para el vector X. Por ejemplo, iniciando con el vector cero, es posible calcular un nuevo vector X usando las ecuaciones para cada X[i] encontradas.

X = [0, 0, 0]

Con cada nuevo valor se calcula el vector diferencias entre el vector X y cada nuevo valor calculado X[i] .

El error que llama la atención es al mayor valor de las diferencias, y se toma como condición para repetir la evaluación de cada vector.

     nuevo = [ 0, 0,  0]
         X = [x0, x1, x2]
diferencia = [|0 - x0|, |0 - x1|, |0 - x2|]
errado = max(diferencia)

Se observa los resultados de errado para cada iteración, relacionados con la convergencia. Si luego de «muchas» iteraciones se encuentra que (errado>tolera),  se detiene el proceso.

Con lo descrito, se muestra una forma de implementar el algoritmo como una función y un ejemplo de prueba con datos del problema conocido.

# Algoritmo Gauss-Seidel,
# matrices, métodos iterativos
# Referencia: Chapra 11.2, p.310, pdf.334
#      Rodriguez 5.2 p.162
# ingresar iteramax si requiere más iteraciones

import numpy as np

def gauss_seidel(A,B,tolera,X,iteramax=100):
    tamano = np.shape(A)
    n = tamano[0]
    m = tamano[1]
    diferencia = np.ones(n, dtype=float)
    errado = np.max(diferencia)
    
    itera = 0
    while not(errado<=tolera or itera>iteramax):
        for i in range(0,n,1):
            nuevo = B[i]
            for j in range(0,m,1):
                if (i!=j): # excepto diagonal de A
                    nuevo = nuevo-A[i,j]*X[j]
            nuevo = nuevo/A[i,i]
            diferencia[i] = np.abs(nuevo-X[i])
            X[i] = nuevo
        errado = np.max(diferencia)
        itera = itera + 1
    # Vector en columna
    X = np.transpose([X])
    # No converge
    if (itera>iteramax):
        X=0
    return(X)

# Programa de prueba #######
# INGRESO
A = np.array([[3,-0.1,-0.2],
              [0.1,7,-0.3],
              [0.3,-0.2,10]])

B = np.array([7.85,-19.3,71.4])

tolera = 0.00001

# PROCEDIMIENTO
n = len(B)
X = np.zeros(n, dtype=float)
respuesta = gauss_seidel(A,B,tolera,X)
verifica = np.dot(A,respuesta)

# SALIDA
print('respuesta de A.X=B : ')
print(respuesta)
print('verificar A.X: ')
print(verifica)

que da como resultado:

respuesta de A.X=B : 
[[ 3. ]
 [-2.5]
 [ 7. ]]
verificar A.X: 
[[  7.84999999]
 [-19.3       ]
 [ 71.4       ]]
>>> 

que es la respuesta del problema obtenida durante el desarrollo del ejemplo con otros métodos.