Métodos iterativos de Jacobi y Gauss-Seidel

01 · Otra forma de resolver

De una factorización exacta a una sucesión de aproximaciones

Gauss y LU son métodos directos: en aritmética exacta producen la solución después de una cantidad finita de operaciones. Los métodos iterativos comienzan con una estimación x⁽⁰⁾ y construyen:

El objetivo es que la sucesión se acerque a la solución de Ax=b. Podemos detenerla cuando el residuo o el cambio entre iteraciones alcance la tolerancia requerida.

02 · Cuándo convienen

Sistemas grandes y matrices dispersas

Una matriz es dispersa cuando la mayoría de sus coeficientes son cero. Este patrón aparece al discretizar temperaturas, difusión, redes, estructuras, imágenes y ecuaciones diferenciales.

No siempre convergen, por lo que la estructura de A, el orden de las ecuaciones y el criterio de parada son fundamentales.

03 · Construcción

Despejar cada incógnita

La fila i del sistema es:

Si aᵢᵢ≠0, podemos despejar:

Al evaluar repetidamente estas expresiones obtenemos un proceso de punto fijo. Jacobi y Gauss-Seidel se diferencian en qué versión de cada componente utilizan durante una iteración.

04 · Actualización simultánea

Método de Jacobi

Jacobi calcula todos los componentes nuevos usando únicamente los valores de la iteración anterior:

Necesita dos vectores: uno conserva x⁽ᵏ⁾ y otro recibe x⁽ᵏ⁺¹⁾. Esto evita mezclar iteraciones.

Como cada componente nuevo depende del mismo vector anterior, muchos cálculos pueden realizarse en paralelo.

05 · Ejemplo manual

Primeras iteraciones de Jacobi

Consideremos el sistema, cuya solución es (1,2,−1):

Despejamos:

Desde x⁽⁰⁾=(0,0,0):

En la segunda iteración, los tres cálculos usan exclusivamente (0.6,2,−1), aunque un componente nuevo ya haya sido obtenido.

06 · Actualización inmediata

Método de Gauss-Seidel

Gauss-Seidel usa inmediatamente los componentes recién calculados. En la fila i, utiliza valores de la iteración nueva para j<i y valores anteriores para j>i:

Con el mismo sistema y x⁽⁰⁾=(0,0,0), la primera iteración es:

1. x⁽¹⁾=0.6.
2. y⁽¹⁾=(22+0.6+0)/11≈2.054545, usando el nuevo x.
3. z⁽¹⁾=(−10−2(0.6)+2.054545)/10≈−0.914545, usando los nuevos x e y.

Esta reutilización suele acelerar la convergencia, pero introduce una dependencia secuencial entre componentes.

07 · Diferencia esencial

Jacobi frente a Gauss-Seidel

Que Gauss-Seidel sea frecuentemente más rápido no significa que siempre converja cuando Jacobi converge, ni que cualquier orden de las filas produzca el mismo comportamiento.

08 · Matrices de iteración

Separar A en diagonal, parte inferior y superior

Escribimos A=D+L+U, donde D contiene la diagonal, L la parte estrictamente inferior y U la superior.

Ambos tienen la forma x⁽ᵏ⁺¹⁾=Tx⁽ᵏ⁾+c. La convergencia desde cualquier valor inicial ocurre exactamente cuando el radio espectral de la matriz de iteración satisface:

El radio espectral es la mayor magnitud de los valores propios de T. En la práctica suelen utilizarse condiciones suficientes más fáciles de comprobar.

09 · Condiciones útiles

Cuándo podemos esperar convergencia

Dominancia diagonal estricta

Si para cada fila se cumple:

entonces Jacobi y Gauss-Seidel convergen. Es una condición suficiente, no necesaria: algunas matrices no dominantes también convergen.

Matriz simétrica definida positiva

Si A es simétrica y xᵀAx>0 para todo x≠0, Gauss-Seidel converge. Para Jacobi, esta propiedad por sí sola no constituye la misma garantía general.

Orden de ecuaciones

Una permutación puede colocar coeficientes convenientes en la diagonal y mejorar la dominancia. Sin embargo, reordenar no siempre convierte un sistema divergente en convergente.

10 · Diagnóstico

Residuo, cambio y presupuesto

Una implementación robusta combina más de una señal:

• El residuo comprueba la ecuación original.
• El cambio detecta estabilización o estancamiento.
• El máximo de iteraciones limita el costo.
• La finitud evita continuar después de producir NaN o infinito.
• Un crecimiento sostenido del residuo permite diagnosticar divergencia.

11 · JavaScript

Un solucionador para ambos métodos

function resolverIterativo(A, b, opciones = {}) {
  const n = A.length;
  const metodo = opciones.metodo ?? "jacobi";
  const tolerancia = opciones.tolerancia ?? 1e-10;
  const maxIter = opciones.maxIter ?? 1000;
  let x = opciones.inicial ? [...opciones.inicial] : Array(n).fill(0);

  if (n === 0 || b.length !== n || A.some(fila => fila.length !== n)) {
    throw new Error("Dimensiones incompatibles");
  }
  if (!['jacobi', 'gauss-seidel'].includes(metodo)) {
    throw new Error("Método desconocido");
  }
  if (A.some((fila, i) => !Number.isFinite(fila[i]) || fila[i] === 0)) {
    throw new Error("La diagonal debe ser finita y no nula");
  }

  const normaB = Math.max(1, ...b.map(Math.abs));
  const residuo = vector => Math.max(...A.map((fila, i) =>
    Math.abs(b[i] - fila.reduce((suma, aij, j) =>
      suma + aij * vector[j], 0)))) / normaB;

  for (let iteracion = 1; iteracion <= maxIter; iteracion++) {
    const anterior = [...x];
    const siguiente = metodo === 'jacobi' ? Array(n).fill(0) : x;

    for (let i = 0; i < n; i++) {
      let suma = 0;
      for (let j = 0; j < n; j++) {
        if (j === i) continue;
        const xj = metodo === 'gauss-seidel' && j < i
          ? siguiente[j] : anterior[j];
        suma += A[i][j] * xj;
      }
      siguiente[i] = (b[i] - suma) / A[i][i];
    }
    x = siguiente;
    const cambio = Math.max(...x.map((valor, i) =>
      Math.abs(valor - anterior[i]))) / Math.max(1, ...x.map(Math.abs));
    const r = residuo(x);

    if (!x.every(Number.isFinite)) {
      return { solucion: x, iteraciones: iteracion,
        convergio: false, motivo: "valor_no_finito" };
    }
    if (r <= tolerancia && cambio <= tolerancia) {
      return { solucion: x, iteraciones: iteracion,
        residuo: r, cambio, convergio: true, motivo: "tolerancia" };
    }
  }
  return { solucion: x, iteraciones: maxIter, residuo: residuo(x),
    convergio: false, motivo: "max_iter" };
}

const A = [[10, -1, 2], [-1, 11, -1], [2, -1, 10]];
const b = [6, 22, -10];
console.log(resolverIterativo(A, b, { metodo: 'gauss-seidel' }));

En una aplicación de producción también conviene validar todos los coeficientes, registrar el historial y detectar crecimiento persistente del residuo.

12 · Implementación eficiente

No recorrer ceros innecesarios

Una matriz densa de n×n requiere recorrer n² coeficientes por iteración. Si solo existen nnz valores no nulos, una estructura dispersa reduce el trabajo a aproximadamente:

En lugar de almacenar cada fila completa, podemos guardar pares [columna,valor]:

function productoDisperso(filas, x) {
  return filas.map(fila =>
    fila.reduce((suma, [columna, valor]) =>
      suma + valor * x[columna], 0)
  );
}

const filas = [
  [[0, 10], [1, -1], [2, 2]],
  [[0, -1], [1, 11], [2, -1]],
  [[0, 2], [1, -1], [2, 10]]
];
console.log(productoDisperso(filas, [1, 2, -1])); // [6, 22, -10]

Formatos como CSR almacenan todos los valores, índices de columna y comienzos de fila en arreglos compactos. Son habituales en bibliotecas numéricas.

13 · Ajustar la actualización

Relajación y método SOR

Podemos combinar el valor anterior con la propuesta de Gauss-Seidel mediante un parámetro ω:

0<ω<1

Subrelajación: pasos más prudentes, útiles para amortiguar oscilaciones.

ω=1

Recupera exactamente el método de Gauss-Seidel.

1<ω<2

Sobrerrelajación: puede acelerar ciertos problemas; una mala elección puede empeorar o destruir la convergencia.

Elegir ω

Depende de la matriz. No existe un valor universalmente óptimo.

SOR muestra una idea general de los métodos iterativos: además de la fórmula, importa cómo se controla y precondiciona el proceso.

14 · Exploración interactiva

Laboratorio: comparar la caída del residuo

Elegí un sistema, un valor inicial y la cantidad de iteraciones. El gráfico muestra el residuo relativo de Jacobi y Gauss-Seidel en escala logarítmica.

Preguntas para explorar

1. En el sistema dominante, compará las pendientes de ambas curvas.
2. Cambiá el valor inicial. ¿Modifica la solución límite o principalmente la trayectoria?
3. Observá que la dominancia diagonal es suficiente, pero no necesaria, usando el sistema simétrico.
4. En los casos desfavorable y divergente, buscá crecimiento u oscilación del residuo.

15 · Comprensión

Confusiones frecuentes

16 · Práctica guiada

Ejercicios

Ejercicio 1: una iteración de Jacobi

Desde (0,0,0), calculá la primera iteración para el sistema del ejemplo.

Ejercicio 2: Gauss-Seidel

¿Por qué el primer valor de y es aproximadamente 2.054545 y no 2?

Ejercicio 3: dominancia diagonal

¿Es estrictamente dominante la matriz [[4,−1,0],[−1,4,−1],[0,−1,3]]?

Ejercicio 4: criterio de parada

Si el cambio relativo es 10⁻¹² pero el residuo relativo es 10⁻³, ¿deberíamos aceptar la solución con tolerancia 10⁻⁸?

Ejercicio 5: dominancia en código

Implementá una función que compruebe dominancia diagonal estricta por filas.

function esDiagonalDominante(A) {
  return A.every((fila, i) => {
    const fueraDiagonal = fila.reduce((suma, valor, j) =>
      suma + (j === i ? 0 : Math.abs(valor)), 0);
    return Math.abs(fila[i]) > fueraDiagonal;
  });
}

console.log(esDiagonalDominante([
  [4, -1, 0],
  [-1, 4, -1],
  [0, -1, 3]
])); // true

17 · Síntesis

Ideas para recordar

En el próximo tema iniciaremos el bloque de interpolación y aproximación de datos: construiremos funciones a partir de muestras y distinguiremos interpolar de ajustar.