Difusión implícita y algoritmo de Thomas — comprar estabilidad con un solo barrido tridiagonal

Partir la malla a la mitad divide el paso explícito entre cuatro. En una dimensión la salida es un barrido tridiagonal, no un paso menor.

El cálculo dejado a las 11 de la noche amaneció en NaN. La malla acababa de partirse a la mitad. El paso de tiempo, no. Olvidar la línea de la CFL (Courant–Friedrichs–Lewy) costó una hora de clúster.

Si el espaciado se reduce a la mitad, el paso explícito de la difusión queda en un cuarto. Una malla de 1024 puntos igual camina cientos de miles de pasos, solo para no explotar. La salida es corta. El futuro pasa al lado derecho y Thomas factoriza la tridiagonal.

El deslizador de r va de 0.05 a 5 y superpone el mismo pulso gaussiano en ambos esquemas, así que el punto de ruptura queda en una sola figura.

CFL number r = D·dt / dx²: 0.45Explicit: stable (r ≤ 0.5)Implicit: unconditionally stable

Por debajo de 0.5 las curvas casi coinciden. Al cruzar 0.5 la traza naranja (explícita) empieza a vibrar en diente de sierra; pasada la unidad sale del marco, y la traza cian (implícita) se queda.

La CFL ata las mallas finas#

La ecuación de difusión es

\frac{\partial q}{\partial t} = D \frac{\partial^2 q}{\partial x^2}

donde $q$ es concentración (o temperatura) y $D$ el coeficiente de difusión. El esquema explícito más simple (Euler hacia adelante con diferencias centrales) se escribe

q_i^{n+1} = q_i^{n} + r\,(q_{i+1}^{n} - 2 q_i^{n} + q_{i-1}^{n}),\quad r = \frac{D\,\Delta t}{\Delta x^2}

Para que sea estable hace falta $r \le 1/2$ . Análisis de von Neumann en una línea: la longitud de onda más corta (longitud $= 2\Delta x$ ) tiene factor de amplificación $|1 - 4r|$ , y forzarlo $\le 1$ obliga a $r \le 1/2$ .

Pasar de $1024 \to 2048$ deja $\Delta t$ en $1/4$ . Un corte más, hasta $4096$ , deja $\Delta t$ en $1/16$ respecto del paso de 1024 puntos, así que el número de pasos se multiplica por 16. Ocho veces más puntos significan 64 veces más pasos para el mismo instante. La factura es ese cuadrado, no el número de celdas.

Con $r = 0.6$ el factor de la onda más corta es $|1 - 4r| = 1.4$ . Sesenta pasos son ese factor elevado a 60. El $1.33\times 10^{12}$ impreso también arrastra la parte del pulso inicial que cayó sobre ese modo. Con $r = 0.4$ el factor es $0.6$ y el pico explícito se queda en $4.21\times 10^{-1}$ .

Mover el futuro al lado derecho#

Una sola sustitución. $q^n$ del lado derecho pasa a ser $q^{n+1}$ .

q_i^{n+1} = q_i^{n} + r\,(q_{i+1}^{n+1} - 2 q_i^{n+1} + q_{i-1}^{n+1})

Con el valor futuro a la derecha ya no se resuelve celda a celda. Todas las incógnitas quedan acopladas y el paso es una ecuación matricial.

A\,\mathbf{q}^{n+1} = \mathbf{q}^{n}

El mismo análisis da un factor $1/(1 + 4r)$ , siempre menor que 1 si $r$ es positivo. Aumentar $\Delta t$ no hace divergir el esquema. Es incondicionalmente estable.

El precio es la precisión. Un $\Delta t$ enorme sigue estable y deja de ser preciso. Aun así, "no va a explotar" pesa cuando el trabajo queda corriendo de noche.

Con el mismo $r = 0.4$ y los mismos 60 pasos, los picos son $4.21\times 10^{-1}$ explícito y $0.4456$ implícito. Ambos estables, y ya discrepan en el tercer dígito. Mismo tiempo físico. La diferencia es error de truncamiento.

Matrices tridiagonales y el barrido único de Thomas#

Al escribir $A$ , en 1D la celda $i$ solo se acopla con $i-1$ e $i+1$ .

A = \begin{pmatrix} 1+r & -r & & & \\ -r & 1+2r & -r & & \\ & -r & 1+2r & -r & \\ & & \ddots & \ddots & \ddots \\ & & & -r & 1+r \end{pmatrix}

Es una matriz tridiagonal. Solo tres diagonales son distintas de cero. La eliminación gaussiana general es $\mathcal{O}(N^3)$ ; una tridiagonal es $\mathcal{O}(N)$ . El algoritmo es el de Thomas. Llewellyn Thomas lo dejó en una nota interna de IBM en 1949.

Thomas da dos pasadas.

Eliminación hacia adelante — barrido de arriba abajo, eliminando la subdiagonal y calculando una nueva diagonal $b'$ y un nuevo lado derecho $d'$ .
Sustitución hacia atrás — barrido de abajo arriba, sustituyendo los valores ya conocidos.

Cada pasada son cuatro o cinco operaciones por celda, y la memoria son tres vectores. La matriz completa no tiene que residir en memoria.

Solo la diagonal de frontera vale $1+r$ . Las filas interiores valen $1+2r$ . Flujo nulo iguala la celda fantasma al valor de frontera, así que el vecino exterior desaparece. Un $r$ sale del $2r$ diagonal, y el $-r$ que apuntaba a ese vecino se va con él. Por eso el código fija b[0] = 1 + r. Copiar la fila interior en la frontera hace que el total se escape.

Comparación en NumPy#

import numpy as np
 
def explicit_diffusion_step(q, r):
    """Euler hacia adelante + diferencias centrales."""
    next_q = q.copy()
    next_q[1:-1] = q[1:-1] + r * (q[2:] - 2 * q[1:-1] + q[:-2])
    return next_q
 
def thomas_solver(a, b, c, d):
    """Resuelve el sistema tridiagonal A x = d en O(N).
    a: subdiagonal (a[0] no se usa)
    b: diagonal principal
    c: superdiagonal (c[-1] no se usa)
    d: lado derecho
    """
    n = len(d)
    cp = np.zeros(n)
    dp = np.zeros(n)
    cp[0] = c[0] / b[0]
    dp[0] = d[0] / b[0]
    for i in range(1, n):
        m = b[i] - a[i] * cp[i - 1]
        cp[i] = c[i] / m
        dp[i] = (d[i] - a[i] * dp[i - 1]) / m
    x = np.zeros(n)
    x[-1] = dp[-1]
    for i in range(n - 2, -1, -1):
        x[i] = dp[i] - cp[i] * x[i + 1]
    return x
 
def implicit_diffusion_step(q, r):
    """Euler hacia atrás — Thomas resuelve (I - r·L) q^{n+1} = q^n."""
    n = len(q)
    a = np.full(n, -r)
    b = np.full(n, 1 + 2 * r)
    c = np.full(n, -r)
    b[0] = 1 + r        # frontera de flujo cero
    b[-1] = 1 + r
    return thomas_solver(a, b, c, q.copy())
 
# Comparar sobre un pulso gaussiano
N = 80
x = np.linspace(0, 1, N)
q0 = np.exp(-((x - 0.5) ** 2) / 0.005)
D, dx = 1.0, 1.0 / N
r_values = [0.4, 0.6, 2.0]
for r in r_values:
    dt = r * dx**2 / D
    qe, qi = q0.copy(), q0.copy()
    for _ in range(60):
        qe = explicit_diffusion_step(qe, r)
        qi = implicit_diffusion_step(qi, r)
    print(f"r={r:.2f} | explicit max={np.max(np.abs(qe)):.2e} | implicit max={np.max(np.abs(qi)):.4f}")
# r=0.40 | explicit max=4.21e-01 | implicit max=0.4456
# r=0.60 | explicit max=1.33e+12 | implicit max=0.3873
# r=2.00 | explicit max=2.78e+38 | implicit max=0.2517

Con $r=0.6$ el método explícito ha llegado a $10^{12}$ en 60 pasos. El implícito, con el mismo paso, sigue difundiendo, con un pico de $0.3873$ . Con $r=2$ el máximo explícito es $2.78\times 10^{38}$ y el pico implícito ha bajado a $0.2517$ .

Las tres filas impresas no son el mismo instante#

Con $\Delta t = r\,\Delta x^2 / D$ , $N = 80$ y $D = 1$ , sesenta pasos de tiempo físico escalan con $r$ . La fila $r = 2$ está cinco veces más adelante que la fila $r = 0.4$ . Leer $0.4456$ junto a $0.2517$ y decir que el implícito "murió más" es una comparación falsa. Igualar el tiempo exige igualar $r$ por el número de pasos. Doce pasos con $r = 2$ caen en el mismo instante que sesenta pasos con $r = 0.4$ . La salida impresa no hizo ese ajuste.

Por qué el implícito no siempre es la respuesta#

La estabilidad incondicional no es un pase libre. El coste se parte en tres.

1. Precisión vs estabilidad — La estabilidad solo dice que la amplitud no explota. No dice que la respuesta sea correcta. Euler hacia atrás es $\mathcal{O}(\Delta t)$ . Si el paso grande también tiene que ser preciso, conviene pasar a Crank–Nicolson ( $\mathcal{O}(\Delta t^2)$ , también incondicionalmente estable).

El factor de la onda más corta en Crank–Nicolson es $(1-2r)/(1+2r)$ . Con $r = 2$ vale $-0.6$ . La magnitud queda por debajo de 1, así que no explota, pero el signo es negativo y la fase se invierte en cada paso. Si $r$ sube más, el factor se pega a $-1$ .

2. No linealidad — Si $D$ depende de la solución ( $D = D(q)$ ), la matriz cambia en cada paso y hace falta un solucionador no lineal (Picard o Newton) encima de Thomas.

3. Más dimensiones — En 2D/3D la matriz deja de ser tridiagonal. Un stencil de cinco o siete puntos crea bandas laterales. El paso se separa con ADI (implícito en direcciones alternas) en dos o tres barridos tridiagonales 1D, o se pasa a los métodos iterativos de Jacobi a BiCGSTAB.

Ese tercer caso reaparece al pasar de flujo compresible a incompresible. La ecuación de Poisson para la presión $\nabla^2 P = -\rho \nabla\!\cdot(\mathbf{u}\!\cdot\!\nabla\mathbf{u})$ tiene la misma forma dispersa y grande. Quien ya resolvió difusión implícita usa las mismas herramientas en el paso de presión.

Anotar un r grande como respuesta correcta es un fallo#

El pico implícito $0.2517$ con $r = 2$ no es un registro de haberse acercado a la solución. El tiempo físico es cinco veces mayor, y el error temporal escala con $\Delta t$ . Guardar la curva sin anotar el tiempo físico de esos pasos vuelve falsa la comparación.

En una malla bidimensional este código tridiagonal no se reutiliza tal cual. Cuando un barrido que nunca forma la matriz parece más barato, ¿hasta dónde lleva la misma promesa de estabilidad LU-SGS, que avanza de forma implícita sin construir la matriz? Si la dirección del barrido está mal, ¿esa promesa se rompe con ella? Esa parte sigue abierta.