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

Más allá del muro de explosión del esquema explícito, la estabilidad cuesta una pasada O(N)

Una vez lancé una simulación a las 11 de la noche y al volver a la oficina por la mañana abrí el archivo de resultados lleno de NaN. Acababa de duplicar la resolución de la malla y olvidé reducir el paso de tiempo. El precio de olvidar la condición CFL (Courant–Friedrichs–Lewy) fue una hora de tiempo de cluster.

La difusión explícita tiene una ley de escala brutal: si reduces el paso de la malla a la mitad, debes reducir el paso de tiempo a una cuarta parte. Un problema de 1024 puntos cuesta cientos de miles de pasos por la única razón de no explotar. Este artículo trata el atajo de una sola línea — el método implícito — y el algoritmo que en 1D casi elimina su coste: el solucionador tridiagonal de Thomas. Al final podrás empujar el paso de tiempo tú mismo y ver cómo el explícito explota mientras el implícito se mantiene tranquilo.

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$ .

El problema es que esa condición está atada al cuadrado de $\Delta x$ . Pasar de $1024 \to 2048$ celdas recorta $\Delta t$ a la cuarta parte. Es una tiranía de orden de magnitud.

Mover el futuro al lado derecho#

El truco es casi vergonzosamente simple. Sustituye $q^n$ del lado derecho por $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 en el lado derecho, ya no se puede resolver celda a celda. Todas las incógnitas quedan acopladas y aparece una ecuación matricial:

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

A cambio, la recompensa es enorme. Repitiendo el análisis de von Neumann, el factor de amplificación pasa a ser $1/(1 + 4r)$ , siempre menor que 1 para cualquier $r$ positivo. Es decir, el esquema es incondicionalmente estable. Por más que se aumente $\Delta t$ , no diverge.

El precio es la precisión. Un $\Delta t$ desmesuradamente grande sigue siendo estable, pero deja de ser preciso. Aun así, "no va a explotar" es una garantía valiosa para quien tiene que ejecutar el código.

Matrices tridiagonales y el barrido único de Thomas#

Escribamos la matriz $A$ . En 1D, cada celda $i$ se acopla solo con $i-1$ e $i+1$ , así que

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

Esta es una matriz tridiagonal: solo tres diagonales son distintas de cero. La eliminación gaussiana general es $\mathcal{O}(N^3)$ , pero para una tridiagonal se reduce a $\mathcal{O}(N)$ . Este algoritmo se llama de Thomas — Llewellyn Thomas lo dejó por escrito en una nota interna de IBM en 1949.

El núcleo de Thomas son dos fases.

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 paso son cuatro o cinco operaciones por celda, y la memoria se reduce a tres vectores. La matriz completa nunca se materializa.

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. Con el mismo paso de tiempo, el implícito sigue difundiendo suavemente.

Empuja tú mismo el paso de tiempo#

A continuación se puede manipular la simulación. El deslizador cambia $r = D\,\Delta t/\Delta x^2$ entre 0.05 y 5.

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

Por debajo de 0.5 las dos curvas casi se superponen. En cuanto se cruza 0.5, la curva naranja (explícita) empieza a vibrar con un patrón en diente de sierra; pasada la unidad sale disparada de la pantalla. La curva cian (implícita) no se rompe a ningún $r$ . Eso sí, también queda claro el contrapeso: con $r$ muy alto la resolución cae y el pulso se difunde demasiado en un solo paso.

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

La estabilidad incondicional no es comida gratis. Vienen tres costes.

1. Precisión vs estabilidad — La estabilidad solo garantiza que la amplitud no explota; no garantiza que la respuesta sea correcta. Euler hacia atrás es solo de primer orden en tiempo, $\mathcal{O}(\Delta t)$ . Para mantener precisión con $\Delta t$ grande, mejor cambiar a Crank–Nicolson, que es $\mathcal{O}(\Delta t^2)$ y sigue siendo incondicionalmente estable.

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 — iteración de Picard o Newton sobre cada barrido 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 lejos de la diagonal. Entonces se separa el problema con ADI (Implícito en Direcciones Alternas) en barridos 1D tridiagonales sucesivos, o se recurre a solucionadores iterativos — gradiente conjugado, BiCGSTAB, multimalla.

Ese tercer punto reaparece tal cual 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 estructura dispersa. Quien ya resolvió difusión implícita usa la misma caja de herramientas para el paso de presión.

Tres ideas para recordar#

El paso de tiempo de la difusión explícita queda atado a $\Delta x^2$ , lo que se vuelve brutal al refinar la malla.
El método implícito compra estabilidad incondicional moviendo el valor futuro al lado derecho — el coste es una resolución matricial por paso.
En 1D la matriz es tridiagonal y el algoritmo de Thomas termina en $\mathcal{O}(N)$ — incluso un millón de celdas se resuelve en un solo barrido rápido.