Lo que Jacobi no mata en 10 000 barridos, un V-ciclo lo mata en 10 — el verdadero secreto del multigrid

Por qué 10 000 barridos de Jacobi se niegan a terminar#

El ingeniero B resolvía un Poisson de presión en una malla 1024×1024 con Jacobi puntual. A los doscientos segundos el residuo seguía en el mismo orden. El CFL era pequeño. La malla era ortogonal. El código estaba bien. Hay errores que Jacobi no mata.

El mismo error cambia de frecuencia cuando cambia la malla. El V-ciclo es esa frase hecha algoritmo.

Al mover Pre/Post entre 1 y 5, en la misma malla de 64 celdas y el mismo término fuente nulo, solo el residuo del V-ciclo se empina, y el suavizado extra deja de pagar enseguida.

La traza naranja es Jacobi. Las oscilaciones bajan y el montículo ancho se queda. El log-residuo cian cae casi en línea recta, y el naranja se aplana en la meseta de baja frecuencia.

Frecuencias del error — quién muere rápido y quién vive para siempre#

El sistema lineal es $A u = f$. El error de la estimación $u^{(k)}$ es

un vector. El residuo es $r^{(k)} = f - A u^{(k)} = A e^{(k)}$.

En la propagación del Jacobi ponderado, un modo de Fourier $\sin(j \pi x)$ decae por barrido como

donde $\omega$ es el peso de relajación, $j$ el número de modo y $h$ el espaciado. Las magnitudes del factor:

Las ondas cortas bajan a un tercio en cada barrido. Las largas casi no se mueven. Lo que Jacobi no elimina es la colina suave tendida sobre toda la malla.

En una malla más gruesa, esa baja frecuencia se ve alta#

Ahí el multigrid se separa de un suavizador plano. La onda más larga en $N$ celdas es $\lambda \sim N h$. Al pasarla a $N/2$ celdas ocupa la mitad de celdas y se ve más corta. En concreto, el modo más bajo de la malla fina es un modo medio en la gruesa, y ahí Jacobi sí mata los modos medios.

Un nivel más grueso acorta otra vez la onda. Tras $\log_2 N$ niveles cada modo es alta frecuencia en algún sitio. Unos pocos barridos de Jacobi por nivel, y en teoría muere todo error. Dejar los modos bajos del lado de los iterativos es el tema de la convergencia de Jacobi a BiCGSTAB.

El V-ciclo — el camino que baja y vuelve a subir#

Esos pasos dibujan una V.

1. Pre-suavizado — $\nu_1$ barridos de Jacobi en la malla fina.
2. Residuo — $r_h = f_h - A_h u_h$.
3. Restricción — $r_{2h} = R\, r_h$. El residuo baja.
4. Recursión — se resuelve $A_{2h} e_{2h} = r_{2h}$ con otro V-ciclo, o de forma directa si la malla ya es pequeña.
5. Prolongación — $e_h = P\, e_{2h}$. Sube la corrección.
6. Corrección — $u_h \leftarrow u_h + e_h$.
7. Post-suavizado — $\nu_2$ barridos de Jacobi en la malla fina.

Al bajar solo viaja el residuo. Al subir, solo la corrección. Un ciclo cuesta $O(N)$: un barrido fino más una cola geométrica. Un ciclo baja cerca de un dígito de residuo.

Restricción y prolongación — el tipo de cambio entre dos mallas#

Los dos operadores de transferencia tienen que ir en pareja. En una malla 1D centrada en vértices, la pareja simple es esta.

Restricción de ponderación completa: el punto grueso $i$ es el promedio ponderado de los puntos finos $2i-1, 2i, 2i+1$.

Prolongación lineal: se copia en los nodos alineados con la malla gruesa y se promedian los dos vecinos en medio.

Son traspuestas salvo una constante de razón de malla ($P = c R^T$). Esa pareja entra en la garantía de convergencia. Una pareja asimétrica puede divergir.

Un V-ciclo de dos mallas en NumPy#

Se resuelve $-u''(x) = f(x)$ con $u(0)=u(1)=0$. Se parte de 64 celdas.

import numpy as np
 
class GridLevel:
    """Estado de un nivel de malla: solución, lado derecho, espaciado."""
    def __init__(self, n):
        self.n = n              # número de celdas
        self.h = 1.0 / n
        self.u = np.zeros(n + 1)
        self.f = np.zeros(n + 1)
 
def jacobi_relax(u, f, h, iters):
    """Jacobi ponderado para Poisson 1D. Dirichlet cero en ambos extremos."""
    n = len(u) - 1
    nxt = np.zeros_like(u)
    for _ in range(iters):
        nxt[1:n] = 0.5 * (u[0:n-1] + u[2:n+1] + h * h * f[1:n])
        u[1:n] = nxt[1:n]
    return u
 
def residual_norm(u, f, h):
    n = len(u) - 1
    r = np.zeros_like(u)
    inv = 1.0 / (h * h)
    r[1:n] = f[1:n] - (-u[0:n-1] + 2 * u[1:n] - u[2:n+1]) * inv
    return float(np.sqrt(np.mean(r[1:n] ** 2)))
 
def restrict_full_weighting(rh):
    n = len(rh) - 1
    n2 = n // 2
    r2 = np.zeros(n2 + 1)
    r2[1:n2] = 0.25 * rh[1:n-1:2] + 0.5 * rh[2:n:2] + 0.25 * rh[3:n+1:2]
    return r2
 
def prolong_linear(e2h):
    n2 = len(e2h) - 1
    n = 2 * n2
    eh = np.zeros(n + 1)
    eh[2:n:2] = e2h[1:n2]
    eh[1:n+1:2] = 0.5 * (e2h[0:n2] + e2h[1:n2+1])
    return eh
 
def v_cycle(u, f, h, n_pre=2, n_post=2):
    n = len(u) - 1
    if n <= 4:
        return jacobi_relax(u, f, h, 50)
    jacobi_relax(u, f, h, n_pre)
    inv = 1.0 / (h * h)
    r = np.zeros_like(u)
    r[1:n] = f[1:n] - (-u[0:n-1] + 2 * u[1:n] - u[2:n+1]) * inv
    r2 = restrict_full_weighting(r)
    e2 = np.zeros_like(r2)
    v_cycle(e2, r2, 2 * h, n_pre, n_post)
    e = prolong_linear(e2)
    u += e
    jacobi_relax(u, f, h, n_post)
    return u
 
# --- ejecución ---
n = 64
h = 1.0 / n
x = np.linspace(0, 1, n + 1)
f = (np.pi ** 2) * np.sin(np.pi * x)   # solución exacta u = sin(pi x)
 
u_jac = np.zeros(n + 1)
u_mg = np.zeros(n + 1)
 
for k in range(20):
    jacobi_relax(u_jac, f, h, 4)
    v_cycle(u_mg, f, h, n_pre=2, n_post=2)
    print(f"step {k:3d} | jacobi r={residual_norm(u_jac, f, h):.3e} "
          f"v-cycle r={residual_norm(u_mg, f, h):.3e}")

Un step impreso gasta 4 barridos de Jacobi, y el V-ciclo gasta 4 barridos finos (2 pre + 2 post). Trabajo parecido, residuos distintos. Tras 20 pasos el residuo del V-ciclo suele quedar bajo $10^{-10}$ y el de Jacobi se queda cerca de $10^{-2}$.

La tabla cierra en $\omega = 2/3$. El modo más rápido, $jh = 1$, tiene $\cos(\pi) = -1$, así que el factor es $1 - 2\omega = -1/3$. El valor absoluto coincide con el $1/3$ de la tabla. Diez barridos solo de ese modo dejan $(1/3)^{10} \approx 1.7 \times 10^{-5}$. El modo medio es $\approx 0.67$, y $0.67^{20} \approx 3.3 \times 10^{-4}$. El factor más lento es $1 - O(h^2)$. Con $N = 64$ y $h = 1/64$, $h^2 = 1/4096 \approx 2.4 \times 10^{-4}$. Ese es el tamaño de un corte sobre la onda larga.

La salida encaja con la tabla. Cuatro barridos de Jacobi por paso, veinte pasos, ochenta barridos, y el residuo sigue cerca de $10^{-2}$. El modo rápido ya murió con factor $1/3$. El $10^{-2}$ que queda es el modo bajo de la tabla. Los mismos veinte V-ciclos, 2 pre y 2 post en la malla fina, suelen quedar bajo $10^{-10}$. Ocho órdenes separan $10^{-2}$ y $10^{-10}$.

El nivel más grueso son 50 barridos de Jacobi cuando $n \le 4$. Sesenta y cuatro celdas se restringen como 64, 32, 16, 8, 4: cuatro transferencias. La ponderación completa es $1/4$, $1/2$, $1/4$. El trabajo grueso suma $N/2 + N/4 + \cdots$ y no pasa de $N$. Diez V-ciclos repiten esa suma diez veces. Diez mil barridos de Jacobi son otro recuento, antes de mirar la constante.

El mismo factor, puesto en el modo más lento, da cerca de $0.9992$ con $N = 64$, $jh = 1/N$ y $\omega = 2/3$. Los ochenta barridos del listado multiplican hasta cerca de $0.94$. Un dígito cuesta unos $2900$ barridos. Con $N = 1024$ el factor queda cerca de $0.999997$ y un dígito cuesta unos $7.3 \times 10^{5}$ barridos. Diez mil barridos dejan ese modo cerca de $0.97$. La corrida de 1024×1024 que se quedó en un solo orden durante 200 segundos es el mismo atasco. La fórmula es Jacobi ponderado en 1D. La malla del ingeniero era 2D. Un V-ciclo que baja un dígito por ciclo hace en unos diez ciclos lo que esos barridos no hacen con el modo lento.

El $0.9992$ no es un residuo medido. Sale de la fórmula de la tabla en el modo más lento. Ochenta barridos lo mueven cerca de un seis por ciento. El residuo de Jacobi impreso, cerca de $10^{-2}$ tras veinte pasos, es lo que queda cuando un factor tan próximo a 1 se queda solo, después de que el modo rápido — amortiguado a $1/3$ — ya desapareció. En $N = 1024$ el mismo factor deja diez mil barridos cerca de $0.97$, de modo que el orden de magnitud no se mueve.

Por qué no siempre es más rápido — dos trampas#

El V-ciclo no es magia. En la práctica se rompe en dos sitios.

Coeficientes abruptos. En un Poisson anisótropo o con saltos, Jacobi no es un buen suavizador. En $-\nabla \cdot (k(x) \nabla u) = f$, si $k$ salta por órdenes de magnitud el V-ciclo estándar deja el residuo en algún dígito. La salida es line-Jacobi, una malla alineada, o multigrid algebraico, que arma la jerarquía desde la matriz.

Mallas no anidadas. En mallas no estructuradas o AMR el paso de fino a grueso no es partir un índice. El multigrid por aglomeración agrupa celdas finas en celdas gruesas, restringe con pesos de volumen y prolonga con interpolación centrada en la celda. GAMG de OpenFOAM sigue ese patrón.

La tasa es independiente de la malla solo si el suavizador, el operador y la razón de engrosamiento coinciden. Si falta uno, el número de iteraciones crece con $N$. Ahí entra un W-ciclo o un F-ciclo.

Dónde dejar el V-ciclo estándar#

Si $k$ salta por órdenes de magnitud, el Jacobi plano no es un suavizador y el residuo se pega a algún dígito. El cambio es line-Jacobi o AMG. La versión para mallas estiradas está en un precondicionador linelet.

Una malla no estructurada o AMR no fabrica el nivel grueso partiendo un índice. GAMG de OpenFOAM aglomera celdas. Si el número de iteraciones empieza a seguir a $N$, el siguiente ensayo es un W-ciclo o un F-ciclo. Las 64 celdas, los 50 barridos de Jacobi en $n \le 4$ y la brecha de 20 pasos entre $10^{-10}$ y $10^{-2}$ se leyeron solo en un Poisson 1D alineado.