Cuando las características chocan, nace una onda de choque — Burgers y el plano (x,t)

En 1948, el físico holandés J. M. Burgers lanzó una ecuación de una sola línea como modelo de turbulencia. Una curva suave que, sin fuerza externa, se rompe en un acantilado vertical en tiempo finito — el drama más corto de las EDP hiperbólicas no lineales. Este artículo lee a mano el instante de formación del choque a partir de la geometría de las características (las curvas espacio-temporales por donde viaja la información), muestra cómo una solución débil (definida por la forma integral allí donde la diferenciabilidad falla) corta la región multivaluada con una sola línea, y reproduce la misma imagen en 60 líneas de Python. Al final podrás manipular el plano (x,t) y ver chocar las líneas en directo.

Dos líneas se encuentran en un mismo punto#

Las ecuaciones hiperbólicas tratan de propagación de señales. Mientras la información fluya por características a velocidades fijas, todo es tranquilo. En la advección lineal $\partial_t q + u\,\partial_x q = 0$ con $u$ constante, las características son rectas paralelas y el perfil inicial se desplaza lateralmente. En esa imagen no cabe ninguna onda de choque.

Todo cambia cuando aparece la no linealidad. Si $u$ depende de $q$ mismo, los valores grandes corren rápido y los pequeños despacio. En el instante en que un valor rápido alcanza a uno lento, dos valores llegan al mismo punto $(x,t)$. La función se vuelve multivaluada, físicamente imposible — ese sinsentido es justo el nacimiento de un choque.

Una ecuación que se transporta a sí misma#

La forma conservativa de la ecuación de Burgers es:

$q$ es el escalar conservado y $f(q) = q^2/2$ es el flujo. La regla de la cadena da la forma no conservativa:

La velocidad de advección es $u(x,t) = q(x,t)$. La cantidad se transporta a sí misma. Esa única no linealidad es el origen de todo el drama. El jacobiano $\partial f/\partial q = q$ es positivo donde la señal va a la derecha y negativo donde va a la izquierda. Un perfil inicial que cambia de signo en el espacio — la onda sinusoidal es el caso típico — obliga a los puntos rápidos a perseguir a los lentos.

El plano (x,t) que dibujan las características#

Una característica es una curva que cumple $dx/dt = q$. Como sobre ella la derivada material $Dq/Dt$ es nula, $q$ queda congelado en su valor inicial a lo largo de la curva. La característica resulta ser una recta en el plano (x,t) con pendiente $1/q_0(x_0)$; el valor inicial $q_0(x_0)$ es la velocidad de esa recta.

Para la condición inicial $q_0(x_0) = 0.55 + A\sin(2\pi k x_0)$:

• Las rectas que parten de los picos ($q_0 > 0.55$) se inclinan fuertemente a la derecha (rápidas).
• Las que parten de los valles ($q_0 < 0.55$) se inclinan poco (lentas).

La línea del pico acabará pisando la punta de la línea del valle inmediatamente delante. El punto donde dos líneas se cruzan es el primer choque.

El tiempo de ruptura t* — el primer cruce#

Dos características vecinas, partiendo de $x_0$ y $x_0 + \Delta x$, se encuentran en el tiempo $t$ cuando:

Tomando el límite $\Delta x \to 0$:

Sólo las pendientes negativas producen un encuentro. El choque más temprano ocurre donde $q_0'(x_0)$ es más negativa:

Cuanto mayores la amplitud o el número de onda de $q_0$, antes se rompe. Para el seno anterior, $q_0'(x_0) = 2\pi k A\cos(2\pi k x_0)$, con mínimo $-2\pi k A$:

Con $A = 0.35$, $k = 1$, sale $t^* \approx 0.455$. Antes la curva se desfasa suavemente; al cruzar $t^*$ se levanta una pared vertical.

Solución débil y Rankine–Hugoniot — el cuchillo que corta#

Más allá de $t^*$, la lectura ingenua escupe tres valores de $q$ en el mismo $x$ — imposible. La salida es la solución débil: dejar la forma diferencial y volver a la ley de conservación integral.

Esto tiene sentido aun cuando $q$ no sea derivable — es decir, aunque sea discontinuo. Si dos estados planos $q_L$ y $q_R$ están unidos por un choque que avanza a velocidad $s$:

Es la condición de Rankine–Hugoniot. La velocidad del choque es la media aritmética de los dos estados. Esta condición corta la zona multivaluada con una línea y exige que las áreas a ambos lados sean iguales — la regla de áreas iguales (equal-area rule).

Hace falta además una condición de entropía. El choque sólo es físico si $q_L > q_R$. Si $q_L < q_R$, la respuesta es una rarefacción en abanico (rarefaction). Si se ignora esta condición, queda como respuesta una expansión-choque matemáticamente válida pero no física.

Python: del trazado de características al choque de Godunov#

Trazamos las características y atrapamos el choque débil con un esquema de Godunov de primer orden.

import numpy as np
 
def burgers_breaking_time(q0_prime_min):
    # q0_prime_min: minimo de q0'(x) en x (debe ser negativo para tener choque)
    if q0_prime_min >= 0:
        return float('inf')
    return -1.0 / q0_prime_min
 
def trace_characteristics(q0_func, x0_arr, t_arr):
    # x[i, j] = x0_arr[i] + q0(x0_arr[i]) * t_arr[j]
    return x0_arr[:, None] + q0_func(x0_arr)[:, None] * t_arr[None, :]
 
def rankine_hugoniot_speed(qL, qR):
    # velocidad del choque s = (qL + qR) / 2
    return 0.5 * (qL + qR)
 
def godunov_burgers_step(q, dx, dt):
    # forma conservativa: q^{n+1} = q^n - (dt/dx) * (F_{i+1/2} - F_{i-1/2})
    qL = q                          # estado izquierdo de la interfaz derecha
    qR = np.roll(q, -1)             # estado derecho de la interfaz derecha
    s = 0.5 * (qL + qR)
    f_shock = np.where(s >= 0, 0.5 * qL**2, 0.5 * qR**2)
    f_rare  = np.where(qL >= 0, 0.5 * qL**2,
              np.where(qR <= 0, 0.5 * qR**2, 0.0))
    F_iph = np.where(qL > qR, f_shock, f_rare)   # flujo en la interfaz derecha
    F_imh = np.roll(F_iph, 1)                    # interfaz izquierda = desplazar
    return q - (dt / dx) * (F_iph - F_imh)
 
# Parametros de la simulacion --------------------------------
N, L = 256, 1.0
A, k = 0.35, 1
x = (np.arange(N) + 0.5) * (L / N)
q0 = lambda xx: 0.55 + A * np.sin(2 * np.pi * k * xx)
q  = q0(x)
 
t_star = burgers_breaking_time(-2 * np.pi * k * A)
print(f"breaking time t* = {t_star:.3f}")
 
t_end, t = 1.2, 0.0
while t < t_end:
    dt = 0.4 * (L / N) / np.max(np.abs(q))    # CFL = 0.4
    q  = godunov_burgers_step(q, L / N, dt)
    t += dt
print(f"shock height (max - min) at t = {t:.2f}: {q.max() - q.min():.3f}")

Aparece breaking time t* = 0.455 y una altura de choque cercana a 0.7. El esquema de Godunov corta automáticamente la región multivaluada en una sola línea — el cuchillo es la línea np.where(qL > qR, f_shock, f_rare). Esa línea resuelve el problema de Riemann exacto en cada interfaz y garantiza la solución débil.

Mira chocar las líneas#

En la simulación de abajo cambia la amplitud y el número de onda.

Si subes la amplitud $A$ a 0.45, $t^*$ se reduce a unos 0.35 y la línea roja punteada sube. Con número de onda $k = 2$, la ruptura ocurre en dos lugares a la vez; una vez formados, los dos choques se fusionan en uno mayor. En el panel superior, la curva cian y la curva punteada tenue (perfil inicial) encierran áreas iguales a cada lado del choque. Es la regla de áreas iguales hecha visible.

Dos motivos por los que falla la captura numérica del choque#

(1) Violación de CFL. Si rompes $\Delta t \le \Delta x / \max|q|$, el dominio numérico adelanta al físico. Si la amplitud parece crecer con el tiempo, sospecha primero del $dt$. El código de arriba recalcula $dt$ en cada paso con np.max(np.abs(q)), lo cual importa aún más justo después de formarse el choque, porque $q$ puede picar localmente.

(2) Entropía mal tratada. Sólo con Rankine–Hugoniot caben choques de expansión no físicos. Solvers linealizados sencillos como Roe producen un "entropy glitch" en rarefacciones transónicas (el abanico cruza el cero). Godunov, HLLE y los Riemann solvers exactos imponen la condición de entropía automáticamente; con Roe hace falta una corrección de entropía explícita estilo Harten.

Cuando un código LES compresible muestra posiciones de choque oscilantes o presiones negativas, casi siempre es uno de estos dos.

La próxima vez que empiece un atasco#

Los coches confluyen en un carril y reducen velocidad. El de delante va lento, el de detrás rápido. En cuanto el rápido toca el parachoques del lento, nace la cabeza del atasco. Es el punto donde se cruzan dos trayectorias en el plano (x,t). Un atasco es prácticamente un choque de Burgers y, hasta que algo lo resuelve, se propaga aguas arriba a la velocidad de Rankine–Hugoniot.

Resumen en una línea. Las EDP hiperbólicas no lineales generan discontinuidades en tiempo finito incluso desde datos iniciales suaves. Aceptar ese hecho y recuperar el sentido por el rodeo de la solución débil es el punto de partida del análisis numérico de los choques.