Cuando FTCS tiembla, el upwind se entumece — esquema advectivo de primer orden, viscosidad artificial y el trueque de CFL

A comienzos de los años sesenta, gente junto a von Neumann metió diferencias centradas en la advección lineal unidimensional. El pulso cuadrado debía derivar en calma. El borde de salida creció en oscilaciones y salió de la malla en pocos pasos.

La causa no es el orden. Traer el valor de aguas abajo con el mismo peso rompe la causalidad. Un upwind de una celda frena esa huida y mata las esquinas. El límite de ese trueque es el número CFL.

Al mover la marca CFL de 0.05 a 1.5, sobre el mismo pulso cuadrado, FTCS sigue temblando y el upwind cae sobre la línea de trazos solo cerca de 1.

Con la marca en 0.5 el pulso upwind se abre más. Pasado 1.05, el upwind y Lax–Wendroff se van juntos del marco. El stencil quedó más corto que la característica.

La diferencia más natural cae primero#

La advección lineal en una dimensión es una sola línea.

$q$ es la cantidad transportada y $u>0$ una velocidad constante. La solución es la forma inicial trasladada a velocidad $u$.

Euler hacia adelante en el tiempo y una diferencia centrada de segundo orden en el espacio dan esta actualización.

$c \equiv u\,\Delta t / \Delta x$ es el número de Courant, las celdas que la onda cruza en un paso. El esquema es de segundo orden en espacio y de primer orden en tiempo, así que el LTE (error local de truncamiento) es $\mathcal{O}(\Delta t,\Delta x^2)$. Por más que se reduzca $c$, el esquema se dispara.

Al sustituir $q^n_i = g^n e^{i k x_i}$, el factor de amplificación es

y $|g|^2 = 1 + c^2 \sin^2(k\Delta x) > 1$. Todo modo crece en cada paso. Basta un modo con $|g|$ mayor que 1 para romper el campo.

Con el $c=0.8$ del código, el modo con $\sin=1$ mide $\sqrt{1+c^2}=\sqrt{1.64}\approx 1.28$. Es 1.28 por paso, y 400 pasos son esa potencia. El $1.78^{400}$ aproximado de más abajo no es este modo. La fórmula da 1.28.

La información solo fluye desde un lado#

Si $u>0$, $q_i$ en $t_{n+1}$ solo puede mirar aguas arriba, en $x<x_i$. La diferencia centrada toma el valor de aguas abajo $q_{i+1}^n$ con el mismo peso. Lee un lado al que el flujo no ha llegado.

El upwind fija esa dirección.

El factor de amplificación es

y $|g|^2 = 1 - c(1-c)\cdot 4\sin^2(k\Delta x/2)$. Para $0 \le c \le 1$, $|g|\le 1$.

El coeficiente $c(1-c)$ llega a 0.25 en $c=0.5$. En $c=0.8$ vale 0.16. En la onda más corta, $|g|^2 = 1 - 4\cdot 0.16 = 0.36$, o sea $|g|=0.6$. En $c=0.5$ el $|g|$ de ese modo cae a 0. En $c=1$ la resta desaparece y todo modo tiene $|g|=1$. El emborronado máximo está en 0.5, no en el borde del intervalo estable.

Upwind = centrada + una viscosidad oculta#

Al partir la pendiente upwind de primer orden sale la identidad

El término de delante es la diferencia centrada. El de detrás es una segunda derivada discreta. Así, el upwind integra en realidad

El objetivo era la advección. El algoritmo metió una viscosidad $\nu_{\text{num}} = u\Delta x/2$. Viscosidad artificial, también llamada difusión numérica.

Esa viscosidad mata la huida. Las esquinas del pulso cuadrado se abren hasta un espesor $\sqrt{4\nu_{\text{num}} t}$. Un $\Delta x$ menor recorta la viscosidad, pero llegar al mismo tiempo pide más pasos, y el emborronado se queda.

El coeficiente posterior $u\Delta x(1-c)/2$ no coincide con esta fórmula por culpa de $c$. El $u\Delta x/2$ de arriba es el extremo $c\to 0$. Con $c=0.8$, $(1-c)=0.2$, así que el coeficiente que queda es la quinta parte, $u\Delta x/10$. Solo la forma con $(1-c)$ se anula en $c=1$.

CFL — cuando la malla no alcanza la característica#

El valor de $q_i$ en $t_{n+1}$ es el que estaba en $x_i - u\Delta t$ en el tiempo $t_n$. Si ese pie cae entre $x_{i-1}$ y $x_i$, dos valores de malla lo interpolan. En cuanto $u\Delta t > \Delta x$, el pie queda más allá de $x_{i-2}$. El stencil no mira hasta ahí.

Esa geometría es la condición de Courant–Friedrichs–Lewy.

El dominio numérico de dependencia (el stencil) tiene que contener al físico (la característica). CFL es necesaria para todo esquema explícito y nunca suficiente. FTCS puede cumplir CFL y aun así detonar. La explosión es causalidad rota, no un stencil corto.

Tres esquemas en cincuenta líneas de NumPy#

La misma malla, el mismo pulso cuadrado, el mismo CFL.

import numpy as np
 
N = 200
dx = 1.0 / N
u = 1.0
cfl = 0.8
dt = cfl * dx / u
n_steps = 400
 
x = (np.arange(N) + 0.5) * dx
q0 = ((x > 0.2) & (x < 0.4)).astype(float)  # pulso cuadrado
 
def ftcs_advect(q, c):
    return q - 0.5 * c * (np.roll(q, -1) - np.roll(q, 1))
 
def upwind_advect(q, c):
    return q - c * (q - np.roll(q, 1))  # u > 0
 
def lax_wendroff_advect(q, c):
    flux_diff = 0.5 * c * (np.roll(q, -1) - np.roll(q, 1))
    visc = 0.5 * c * c * (np.roll(q, -1) - 2 * q + np.roll(q, 1))
    return q - flux_diff + visc
 
q_ftcs, q_up, q_lw = q0.copy(), q0.copy(), q0.copy()
for _ in range(n_steps):
    q_ftcs = ftcs_advect(q_ftcs, cfl)
    q_up = upwind_advect(q_up, cfl)
    q_lw = lax_wendroff_advect(q_lw, cfl)
 
# Solucion exacta: una vuelta completa volviendo al inicio
q_exact = q0
print("FTCS  L1 err:", np.abs(q_ftcs - q_exact).sum() * dx)
print("Upwind L1 err:", np.abs(q_up - q_exact).sum() * dx)
print("LW    L1 err:", np.abs(q_lw - q_exact).sum() * dx)

Con cfl=0.8 y 400 pasos, las cifras aproximadas que deja este texto son:

• FTCS: la amplitud crece cerca de $1.78^{400}$ y pasa a NaN o inf.
• Upwind: el cuadrado se engorda hacia una forma gaussiana y el error L1 queda cerca de 0.04.
• Lax–Wendroff: la amplitud se conserva, pero quedan rizos (Gibbs) tras las esquinas.

c=0.5 es el que más abre el upwind#

El máximo de $c(1-c)$ en 0.5 es el emborronado de la figura. En $c=1$, $(1-c)$ vale cero y el upwind de primer orden desplaza el cuadrado una celda, intacto. $c=0.8$ queda en medio, en 0.16. Estable no quiere decir menos romo. 0.8 solo es menos viscoso que 0.5.

El primer orden es muy romo, pero el segundo trae nuevos dolores de cabeza#

La viscosidad upwind viene bien junto a un choque. Aplana oscilaciones no físicas. El mismo término sigue encendido en zonas suaves, así que el error queda atado al tamaño de una segunda derivada. Los códigos de producción eligen una de dos salidas.

• Usar un esquema de orden mayor (Lax–Wendroff, MUSCL, WENO) y añadir un limitador de flujo o de pendiente que baje el orden solo cerca de los choques. El cálculo que sujeta con un limitador el temblor de Lax–Wendroff en un escalón es esa elección.
• Llamar un solucionador de Riemann (Godunov, HLLC, Roe) en cada cara e integrar la causalidad de forma directa.

Los dos son mecanismos para decidir dónde se acepta la viscosidad y dónde se rechaza. El upwind de primer orden la acepta en todas las celdas. Simple, seguro y el más difusivo.

Llamar una vuelta a 400 pasos es un fallo#

Con $N=200$, $u=1$ y $c=0.8$, $\Delta t = c/N = 0.004$. Una vuelta de un dominio de longitud 1 pide $1/\Delta t = 250$ pasos, es decir $N/c$. Cuatrocientos pasos recorren $400\times 0.004 = 1.6$. Un pulso de ancho 0.2, desplazado 0.6, pasa de $[0.2,0.4]$ a $[0.8,1.0]$. Esos dos cuadrados no se solapan. La L1 del par nítido es 0.4.

La L1 upwind de 0.04 citada arriba no puede ser esa comparación. Leer 0.04 como el emborronado de una vuelta limpia exige 250 pasos. El upwind de primer orden en una dimensión aguanta hasta $c=1$. Cuando lo que se transporta es una interfaz, ¿por qué el límite cae tanto? Donde la advección de interfaz se queda con un CFL de 0.05 muestra dónde se separan la longitud del stencil y la restricción de la interfaz.