Refinar la malla 8 veces no movió el 2.4674 — dos constantes que Fourier dejó en 1807

Refinar la malla reduce el error L2, pero el amortiguamiento del modo diente de sierra sigue siendo pi^2/4 veces más lento. Hay errores que el refinamiento no toca.

Refiné la malla 8 veces y un número se negó a moverse#

Resolví conducción de calor 1D con 64 celdas. Después 128, 256, 512. El error L2 bajó tal como se esperaba.

Un número impreso al lado no se movió. 2.4674. Las cuatro mallas coincidieron hasta la cuarta cifra decimal.

El espaciado $\Delta x$ se redujo ocho veces. ¿Por qué un error se quedó quieto? Este artículo persigue ese número. La respuesta es $\pi^2/4$ y el punto de partida es París, 1807.

En 1807 los revisores objetaron la convergencia, no la física#

Joseph Fourier presentó en 1807 una memoria sobre conducción de calor en sólidos a la Academia de París. Su afirmación: la distribución de temperatura dentro de un sólido se escribe como suma de senos y cosenos.

Los revisores fueron tres: Lagrange, Laplace, Legendre. La memoria fue rechazada. El motivo no era físico sino matemático: nadie había demostrado que la serie convergiera.

Lagrange intentó aportar él mismo la demostración. No lo consiguió y murió en 1813. La prueba llegó mucho después, de la mano de Dirichlet, alumno del propio Fourier.

Los revisores no estaban siendo quisquillosos. Lo que vieron estaba realmente ahí. Al aplicar una serie de senos a una función con salto, algo sobrevive por muchos términos que se añadan.

En un pulso cuadrado de ancho 0.5, el deslizador N es el número de términos que se conservan. Al subirlo, el pico se pega al salto y la altura casi no se mueve.

terms N = 41Fejer averagezoom near jump

max overshoot = 0.000%  · peak at x = 0.0000 (jump at x = 0.25)  · distance to jump = -0.2500

Al empujar N, el pico se pega a la discontinuidad. La altura se queda cerca del 8.95%. Al activar Fejer average el pico desaparece y el escalón se difumina.

Apilar senos hasta formar una onda cuadrada deja un 8.9490%#

Coloquemos un pulso cuadrado de ancho 0.5 en un intervalo periódico. Sus coeficientes de Fourier salen en forma cerrada.

f(x) = \frac{1}{2} + \sum_{m=1}^{\infty} \frac{2}{\pi m}\sin\!\left(\frac{m\pi}{2}\right)\cos\big(2\pi m (x - 0.5)\big)

Aquí $m$ es el índice del término y $x$ la posición. Solo sobreviven los $m$ impares. Al truncar en $N$ términos, la suma parcial se dispara por encima de 1 justo al lado del salto.

Es el fenómeno de Gibbs. El valor límite es 8.9490%, una constante independiente del número de términos. No es un error que se compre con precisión.

Por qué la ecuación del calor prefiere las series de senos#

Fourier no eligió senos por gusto. En la ecuación del calor

\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}

$\alpha$ es la difusividad térmica. Al sustituir $T = e^{ikx}$ , la segunda derivada devuelve $-k^2 e^{ikx}$ . Es decir, $e^{ikx}$ es autofunción del laplaciano con autovalor $-k^2$ .

Por eso cada modo decae solo, como $e^{-\alpha k^2 t}$ , sin hablar con sus vecinos. Una ecuación en derivadas parciales se descompone en una ecuación ordinaria por modo.

Ese $k^2$ dentro del exponente fija la velocidad de amortiguamiento. Si la longitud de onda se reduce a la mitad, el amortiguamiento es cuatro veces más rápido. Las ondas cortas se borran primero. Lo que hace la difusión cabe en ese único sitio.

Qué devuelve la malla en lugar de $-k^2$ #

Ahora el mismo cálculo sobre una malla. Se reemplaza la segunda derivada por una diferencia central de segundo orden y se sustituye $T_j = e^{ikj\Delta x}$ .

\frac{T_{j+1} - 2T_j + T_{j-1}}{\Delta x^2} = -\frac{4}{\Delta x^2}\sin^2\!\left(\frac{k\Delta x}{2}\right) e^{ikj\Delta x}

Aquí $j$ es el índice de nodo. El operador discreto también tiene $e^{ikj\Delta x}$ como autofunción, lo que significa que la descomposición de Fourier sobrevive intacta a la discretización.

Lo que cambia es el autovalor. No $-k^2$ sino $-k'^2$ , con

k'^2 = \frac{4}{\Delta x^2}\sin^2\!\left(\frac{k\Delta x}{2}\right)

Este $k'$ se llama número de onda modificado (modified wavenumber). Es el esquema confesando qué $k$ cree estar recibiendo.

Con $\theta = k\Delta x$ pequeño, la curva efectiva se pega a la parábola $\theta^2$ . Hacia la derecha se separan. La misma onda muere más tarde en la malla que en el exponente exacto.

Dónde aparece $\pi^2/4$ en el límite de Nyquist#

La longitud de onda más corta que una malla puede llevar es $2\Delta x$ : sube una celda, baja la siguiente, un diente de sierra. Allí $\theta = k\Delta x = \pi$ .

Al sustituir en la fórmula anterior:

\frac{k'^2}{k^2} = \frac{4\sin^2(\pi/2)}{\pi^2} = \frac{4}{\pi^2} = 0.4053

El recíproco es $\pi^2/4 = 2.4674$ , el número que no se movía.

El detalle importante es que $\Delta x$ se canceló. La razón no depende del espaciado. Por mucho que se refine, el diente de sierra de esa malla siempre muere 2.4674 veces más lento.

Subir el orden ayuda poco. Aquí está con esténciles reales.

import numpy as np
 
STENCILS = {
    2: [1.0, -2.0, 1.0],
    4: [-1/12, 4/3, -5/2, 4/3, -1/12],
    6: [1/90, -3/20, 3/2, -49/18, 3/2, -3/20, 1/90],
}
 
def modified_k2(order, theta):
    """Cuadrado del numero de onda efectivo de la 2a derivada centrada (dx=1). El exacto es theta**2."""
    c = STENCILS[order]
    h = len(c) // 2
    sym = sum(c[h + j] * np.cos(j * theta) for j in range(-h, h + 1))
    return -sym
 
theta = np.array([np.pi/8, np.pi/4, np.pi/2, 3*np.pi/4, np.pi])
print("theta/pi   exact k^2    order2    order4    order6    (ratio)")
for th in theta:
    ex = th**2
    row = [modified_k2(o, th) for o in (2, 4, 6)]
    print(f"{th/np.pi:6.3f}  {ex:10.5f}  " + "  ".join(f"{v:7.4f}" for v in row)
          + "   " + " ".join(f"{v/ex:6.4f}" for v in row))
 
print()
for o in (2, 4, 6):
    r = modified_k2(o, np.pi) / np.pi**2
    print(f"order {o}: Nyquist ratio = {r:.4f}   decay is {1/r:.3f}x slower")

theta/pi   exact k^2    order2    order4    order6    (ratio)
 0.125     0.15421   0.1522   0.1542   0.1542   0.9872 0.9997 1.0000
 0.250     0.61685   0.5858   0.6144   0.6166   0.9496 0.9960 0.9996
 0.500     2.46740   2.0000   2.3333   2.4222   0.8106 0.9457 0.9817
 0.750     5.55165   3.4142   4.3856   4.8278   0.6150 0.7900 0.8696
 1.000     9.86960   4.0000   5.3333   6.0444   0.4053 0.5404 0.6124
 
order 2: Nyquist ratio = 0.4053   decay is 2.467x slower
order 4: Nyquist ratio = 0.5404   decay is 1.851x slower
order 6: Nyquist ratio = 0.6124   decay is 1.633x slower

El dial recorre $\theta$ y el orden de esa tabla. Subir el orden pega el lado de baja frecuencia a $\theta^2$ . La razón del diente de sierra apenas se mueve.

theta/pi = 1.00

k'^2 / k^2 = 0.4053  · grid decays 2.467x slower  · a*t/dx^2 = 0.000

Al llevar $\theta$ hasta $\pi$ , la copia de la malla decae más tarde que la onda exacta. En sexto orden la razón de Nyquist sigue cerca de 0.6124, un factor 1.633 demasiado lento.

Las ondas que abarcan cuatro celdas o más ( $\theta/\pi \le 0.25$ ) quedan casi perfectas al subir el orden. El diente de sierra no. Incluso en sexto orden marca 0.6124.

La misma onda cuadrada sobre cuatro mallas#

Ahora toca resolverlo de verdad. Para que no se mezcle el error de integración temporal, el tiempo se integra de forma analítica y queda solo la discretización espacial.

import numpy as np
 
L, ALPHA, TEND = 1.0, 1.0, 2.0e-3
 
def tophat_coeffs(m):
    """Coeficientes complejos de Fourier del pulso 0.25 <= x < 0.75 en el periodo [0,1)."""
    m = np.asarray(m, dtype=float)
    out = np.zeros(m.shape, dtype=complex)
    out[m == 0] = 0.5
    mm = m[m != 0]
    out[m != 0] = (np.exp(-2j*np.pi*mm*0.75) - np.exp(-2j*np.pi*mm*0.25)) / (-2j*np.pi*mm)
    return out
 
def gibbs_overshoot(nterm):
    """Sobreoscilacion maxima (%) de la suma parcial de 2N+1 terminos, muestreada cerca del salto."""
    x = np.linspace(0.25, 0.25 + 4.0/nterm, 40001)
    m = np.arange(-nterm, nterm + 1)
    c = tophat_coeffs(m)
    s = (c[None, :] * np.exp(2j*np.pi*m[None, :]*x[:, None])).sum(axis=1).real
    return (s.max() - 1.0) * 100.0
 
def heat_error(ngrid, smooth):
    """Ecuacion del calor integrada exactamente en el tiempo; solo cambia la discretizacion espacial."""
    dx = L / ngrid
    x = np.arange(ngrid) * dx
    m = np.fft.fftfreq(ngrid, d=1.0/ngrid)
    k = 2*np.pi*m
    kmod2 = (4.0/dx**2) * np.sin(k*dx/2.0)**2          # k^2 que entrega la malla
    mref = np.arange(-8192, 8193)
    if smooth:
        t0 = np.cos(2*np.pi*x) + 0.3*np.cos(6*np.pi*x)
        ref = (np.cos(2*np.pi*x)*np.exp(-ALPHA*(2*np.pi)**2*TEND)
               + 0.3*np.cos(6*np.pi*x)*np.exp(-ALPHA*(6*np.pi)**2*TEND))
    else:
        t0 = ((x >= 0.25) & (x < 0.75)).astype(float)
        c = tophat_coeffs(mref)*np.exp(-ALPHA*(2*np.pi*mref)**2*TEND)
        ref = (c[None, :]*np.exp(2j*np.pi*mref[None, :]*x[:, None])).sum(axis=1).real
    num = np.fft.ifft(np.fft.fft(t0)*np.exp(-ALPHA*kmod2*TEND)).real
    nyq = ngrid//2
    return np.sqrt(np.mean((num-ref)**2)), 1.0/(ALPHA*k[nyq]**2), 1.0/(ALPHA*kmod2[nyq])
 
print("N_terms   Gibbs overshoot[%]")
for n in (16, 32, 64, 128, 256, 512):
    print(f"{n:6d}    {gibbs_overshoot(n):9.5f}")
 
for smooth in (True, False):
    print(f"\ninitial condition: {'smooth cosine' if smooth else 'square wave (discontinuous)'}")
    print("grid      L2 error     order     sawtooth decay time exact/grid    ratio")
    prev = None
    for ng in (64, 128, 256, 512):
        l2, t_ex, t_gr = heat_error(ng, smooth)
        order = "  -  " if prev is None else f"{np.log2(prev/l2):5.2f}"
        print(f"{ng:4d}   {l2:.4e}    {order}     {t_ex:.3e} / {t_gr:.3e}   {t_gr/t_ex:.4f}")
        prev = l2

N_terms   Gibbs overshoot[%]
    16      9.01421
    32      8.96527
    64      8.95306
   128      8.95000
   256      8.94924
   512      8.94905
 
initial condition: smooth cosine
grid      L2 error     order     sawtooth decay time exact/grid    ratio
  64   5.3683e-04      -       2.474e-05 / 6.104e-05   2.4674
 128   1.3424e-04     2.00     6.184e-06 / 1.526e-05   2.4674
 256   3.3562e-05     2.00     1.546e-06 / 3.815e-06   2.4674
 512   8.3907e-06     2.00     3.865e-07 / 9.537e-07   2.4674
 
initial condition: square wave (discontinuous)
grid      L2 error     order     sawtooth decay time exact/grid    ratio
  64   2.3383e-02      -       2.474e-05 / 6.104e-05   2.4674
 128   1.1673e-02     1.00     6.184e-06 / 1.526e-05   2.4674
 256   5.8342e-03     1.00     1.546e-06 / 3.815e-06   2.4674
 512   2.9168e-03     1.00     3.865e-07 / 9.537e-07   2.4674

Errores que bajan con el refinamiento y errores que no#

Con la condición inicial suave el orden de convergencia es exactamente 2.00. El esquema es de segundo orden, tal como promete.

Al darle al mismo esquema una onda cuadrada, cae a 1.00. Una discontinuidad hace que los coeficientes de Fourier decaigan solo como $1/m$ , y esa cola lenta regresa en forma de aliasing. El orden de convergencia lo fija la suavidad de la solución, no el orden del esquema. Apagar la oscilación junto al escalón es el mismo trabajo que armar un esquema TVD de segundo orden con un limitador de flujo.

La razón de la última columna vale 2.4674 en las ocho filas. El error L2 cayó ocho veces; este valor es idéntico hasta el decimal.

No conviene mezclar esos dos hechos. L2 baja porque la amplitud del modo diente de sierra se hace menor, no porque el esquema haya mejorado en tratar dientes de sierra. Refinar una malla no resuelve mejor el problema: lo traslada.

En 64 celdas el coseno suave tiene un L2 de 5.3683e-04. Ocho veces más fino, 512 celdas, queda en 8.3907e-06. La razón es cerca de 64, y el logaritmo sigue diciendo orden 2.00. La onda cuadrada en las mismas mallas baja de 2.3383e-02 a 2.9168e-03, un factor ocho, orden 1.00. El tiempo de amortiguamiento del diente de sierra en 64 celdas es 2.474e-05 exacto frente a 6.104e-05 en la malla. En 512 celdas el par es 3.865e-07 y 9.537e-07. La razón se queda en 2.4674 hasta el cuarto decimal. El $k^2$ exacto en Nyquist es 9.86960. El esténcil de segundo orden devuelve 4.0000, y $4/\pi^2 = 0.4053$ . Una vez que $\Delta x$ se cancela, el refinamiento ya no tiene qué cambiar.

Gibbs se comporta igual. Con 16 términos el sobrepaso es 9.01421%. Con 512 términos es 8.94905%. La diferencia ronda 0.065 puntos porcentuales. En la fila $\theta/\pi = 0.125$ , el orden 2 ya marca 0.9872, así que dieciséis celdas por longitud de onda quedan dentro del 1%. El sexto orden en Nyquist devuelve 6.0444, razón 0.6124, amortiguamiento 1.633 veces demasiado lento. Comprar orden compra cifras en baja frecuencia. No compra el diente de sierra.

Al dividir el L2 de la onda cuadrada en 64 celdas, 2.3383e-02, por el L2 suave, 5.3683e-04, salen unas 44 veces. Misma malla, mismo esquema, solo cambia la condición inicial. En 512 celdas el par es 2.9168e-03 contra 8.3907e-06, unas 348 veces. El refinamiento empeora el caso discontinuo en términos relativos, porque los órdenes 1 y 2 se separan por una potencia de la razón de mallas. 512/64 = 8, así que la razón de errores debe abrirse unas 8 veces. 44 por 8 es 352, y la tabla dice 348.

Por qué sexto orden sigue marcando 0.6124, y para qué sirven los filtros#

Subir el orden solo rinde en números de onda bajos. En la fila $\theta/\pi = 0.125$ el sexto orden da 1.0000; en Nyquist da 0.6124. Ningún esténcil finito trata bien el diente de sierra, porque derivar una onda de $2\Delta x$ es genuinamente ambiguo sobre la malla.

Por eso la práctica se divide en dos caminos. Uno es impedir que el diente de sierra llegue a formarse. Las mallas desplazadas y la interpolación de Rhie–Chow hacen ese trabajo. Apretar solo las celdas malas deja la razón del diente de sierra de esa celda en $\pi^2/4$ . Un quadtree que hace crecer la malla por su cuenta tampoco borra esa razón.

El otro es borrar explícitamente el diente de sierra una vez que existe. Los filtros espectrales, la viscosidad artificial y los pesos de WENO pertenecen a esta familia. El Fejer average de la viz anterior es el miembro más simple, y muestra el precio en la misma imagen: el pico de Gibbs se va y el escalón se difumina.

Ninguna vía es gratis. Lo que se elige no es si el error existe, sino dónde se coloca.

El hueco que Dirichlet rellenó en 1829#

La objeción de Lagrange no era una formalidad. Dudaba de lo que hace una serie cuando se encuentra con un salto, y el 8.9490% estaba realmente ahí.

Dirichlet demostró la convergencia condicional en 1829 y rellenó el hueco. La serie converge en cada punto, solo que no de manera uniforme. El pico de Gibbs nunca fue un contraejemplo: era el aspecto que tiene ese "no uniformemente".

Todo estudio de convergencia de malla repite esa pregunta. ¿Es este error del tipo que baja al refinar, o del tipo en que $\Delta x$ se cancela y $\pi^2/4$ permanece? Una sola fila de orden de convergencia responde únicamente a lo primero.

Si un 2.00 limpio alguna vez dio tranquilidad, vale la pena graficar el número de onda modificado una vez. Lo que un esquema confiesa sobre el diente de sierra no está escrito en ningún otro sitio.