Métodos iterativos: Jacobi y Gauss-Seidel

La idea

Resolver Ax = b por eliminación gaussiana cuesta del orden de n³ operaciones y destruye la estructura dispersa de la matriz. Cuando n es grande (mallas de elementos finitos, redes eléctricas, PageRank) se prefieren los métodos iterativos: partir de un vector inicial y mejorarlo repetidamente hasta que el residuo r = b − Ax sea pequeño. El método de Jacobi despeja cada incógnita de su propia ecuación usando los valores de la iteración anterior: xᵢ⁽ᵏ⁺¹⁾ = (bᵢ − Σⱼ≠ᵢ aᵢⱼ xⱼ⁽ᵏ⁾)/aᵢᵢ. Todas las componentes se actualizan a la vez, lo que lo hace trivialmente paralelo.

El método de Gauss-Seidel hace lo mismo pero usa inmediatamente los valores ya actualizados en la misma pasada, lo que suele acelerar la convergencia (en muchos casos la duplica) y no necesita guardar dos vectores. Ambos son casos de la iteración x⁽ᵏ⁺¹⁾ = T x⁽ᵏ⁾ + c y convergen desde cualquier punto inicial si y solo si el radio espectral ρ(T), el mayor valor propio en módulo de la matriz de iteración, es menor que 1; el error se reduce aproximadamente por un factor ρ(T) en cada paso.

Calcular ρ(T) es caro, por eso se usan condiciones suficientes fáciles de comprobar. La más útil es la dominancia diagonal estricta: si en cada fila |aᵢᵢ| > Σⱼ≠ᵢ |aᵢⱼ|, ambos métodos convergen. Gauss-Seidel converge además para toda matriz simétrica definida positiva. A veces basta reordenar las ecuaciones para conseguir la dominancia; si no se cumple, la iteración puede oscilar o divergir. El método SOR añade un factor de relajación ω que, bien elegido, acelera Gauss-Seidel de forma espectacular, y los métodos de Krylov (gradiente conjugado, GMRES) son la evolución moderna de esta idea.

Fórmulas clave

• Jacobixᵢ⁽ᵏ⁺¹⁾ = (bᵢ − Σ_{j≠i} aᵢⱼ xⱼ⁽ᵏ⁾) / aᵢᵢ
• Gauss-Seidelxᵢ⁽ᵏ⁺¹⁾ = (bᵢ − Σ_{j<i} aᵢⱼ xⱼ⁽ᵏ⁺¹⁾ − Σ_{j>i} aᵢⱼ xⱼ⁽ᵏ⁾) / aᵢᵢ
• Forma matricialx⁽ᵏ⁺¹⁾ = T x⁽ᵏ⁾ + c, converge ⇔ ρ(T) < 1
• Dominancia diagonal estricta|aᵢᵢ| > Σ_{j≠i} |aᵢⱼ| para toda fila i ⇒ convergen Jacobi y Gauss-Seidel
• Residuo‖r⁽ᵏ⁾‖ = ‖b − A x⁽ᵏ⁾‖₂

Dónde se aplica

Los simuladores de fluidos de las películas y los videojuegos resuelven en cada fotograma una ecuación de Poisson para la presión con Gauss-Seidel o con sus descendientes (multimalla), sobre mallas de millones de celdas donde la eliminación gaussiana sería impensable. El cálculo de PageRank es una iteración tipo Jacobi sobre la matriz de la web. En ingeniería civil y térmica, la distribución de temperaturas o de tensiones en una pieza discretizada da sistemas dispersos enormes que se resuelven iterativamente, y los solucionadores de los programas de circuitos y de tomografía médica (reconstrucción ART) usan exactamente estas ideas.

Con programación (Python)

import math

A = [[4, -1, 1], [4, -8, 1], [-2, 1, 5]]
b = [7, -21, 15]          # solución exacta: (2, 4, 3)
n = len(A)

def residuo(x):
    return math.sqrt(sum((b[i] - sum(A[i][j] * x[j] for j in range(n))) ** 2 for i in range(n)))

def jacobi(x):
    return [(b[i] - sum(A[i][j] * x[j] for j in range(n) if j != i)) / A[i][i] for i in range(n)]

def gauss_seidel(x):
    y = x[:]
    for i in range(n):
        y[i] = (b[i] - sum(A[i][j] * y[j] for j in range(n) if j != i)) / A[i][i]
    return y

for i in range(n):
    fuera = sum(abs(A[i][j]) for j in range(n) if j != i)
    print(f"fila {i+1}: |a_ii| = {abs(A[i][i])} > {fuera}  ->", "dominante" if abs(A[i][i]) > fuera else "NO dominante")

xJ = xG = [0.0, 0.0, 0.0]
for k in range(1, 13):
    xJ, xG = jacobi(xJ), gauss_seidel(xG)
    if k <= 4 or k == 12:
        print(f"k={k:2d}  Jacobi {[round(v, 5) for v in xJ]} r={residuo(xJ):.1e}   G-S {[round(v, 5) for v in xG]} r={residuo(xG):.1e}")

La versión en JavaScript, editable y ejecutable en el navegador, está en la página interactiva del tema.

Temas relacionados

Sistemas lineales y eliminación de Gauss · Valores y vectores propios · Iteración de punto fijo · La matriz inversa