Discretización espacial, estabilidad y convergencia

Año 2 · Módulo 15: Dinámica espacial y sistemas continuos · Transferencia e integración

149.1. Del campo continuo a un cálculo que pueda ser auditado

Las ecuaciones en derivadas parciales describen magnitudes que cambian en el tiempo y de un lugar a otro. La concentración de una sustancia depende de la posición donde se mide; una temperatura tiene gradientes y una perturbación en un medio excitable puede recorrer el espacio. Una computadora, en cambio, almacena un conjunto finito de números. Resolver una ecuación espacial exige elegir qué significan esos números y qué propiedades de la ecuación sobrevivirán a la traducción. Determinar si una imagen calculada se parece a un fenómeno conocido es solo una parte de la tarea; también debemos establecer qué problema aproxima la simulación y cómo cambia la respuesta cuando mejoramos su resolución.

Consideremos la ecuación de difusión en un intervalo de longitud $L$

La dimensión de $D$ es longitud al cuadrado dividida por tiempo. Sin condiciones iniciales y de frontera, la ecuación no define un problema completo. Si fijamos $u(0,t) = u(L,t) = 0$, la magnitud puede escapar por los bordes, según su interpretación física. Si prescribimos flujo cero, la cantidad total se conserva. Con una condición periódica el dominio equivale a un anillo, en el que lo que sale por un extremo entra por el otro. Esta diferencia reaparece en los autovalores del operador espacial y puede decidir si un patrón crece o desaparece, asunto central del caso resuelto C150.

Dividimos el intervalo en puntos $x_{j} = j\Delta x$ (con $\Delta x = L/N$) y elegimos tiempos $t_{n} = n\Delta t$. La notación $U_{j}^{n}$ designa un número que intenta aproximar $u\left( x_{j},t_{n} \right)$. Este significado debe mantenerse explícito. En diferencias finitas, $U_{j}^{n}$ suele ser una aproximación puntual; en volúmenes finitos, puede representar el promedio dentro de una celda. Intercambiar ambas interpretaciones sin modificar las fórmulas produce errores en los balances y en los bordes.

Hay tres decisiones conectadas. La representación espacial establece cómo aproximamos derivadas o flujos; la integración temporal indica cómo avanza el estado discreto; el tratamiento de fronteras cierra el sistema algebraico. El costo de cada elección afecta a las otras. Una malla fina puede exigir pasos temporales extremadamente pequeños; un método implícito permite pasos mayores, aunque cada paso requiere resolver un sistema. Un esquema puede ser estable y disipar artificialmente los patrones que tratamos de estudiar. La convergencia al problema continuo, por su parte, solo tiene sentido cuando este está bien planteado y se aproxima de forma coherente.

149.2. Diferencias finitas, consistencia y el error que introduce una malla

La derivada segunda expresa la curvatura local del campo. En una malla uniforme, una aproximación habitual es

Podemos justificarla mediante expansiones de Taylor. Si $u$ tiene derivadas espaciales suficientemente regulares, alrededor de $x_{j}$ se cumple

Al sumar las expansiones se cancelan las potencias impares y obtenemos

suponiendo regularidad suficiente para la forma indicada del resto. El error de truncamiento espacial es, por tanto, de segundo orden en la malla. Esto quiere decir que, para una función suave, el defecto local de la aproximación a la derivada, evaluado sobre la solución exacta, disminuye aproximadamente cuatro veces cuando reducimos $\Delta x$ a la mitad. Eso no implica todavía que toda la solución numérica converja al mismo ritmo.

Para la derivada temporal podemos usar Euler hacia adelante,

cuya consistencia temporal es de primer orden. Sustituyendo las expresiones obtenemos un esquema con defecto de truncamiento $O\left( \Delta t + (\Delta x)^{2} \right)$. Esta estimación necesita que la solución exacta tenga las derivadas involucradas y que no aparezcan dificultades singulares en las fronteras. Si la condición inicial presenta una discontinuidad, el comportamiento de los primeros pasos puede ser muy diferente del que sugiere una expansión de Taylor ordinaria.

Una discretización es consistente cuando sus ecuaciones discretas se aproximan a las continuas al refinar la malla y el paso, dentro de un régimen apropiado. La consistencia expresa una correspondencia local entre ambas formulaciones, pero no asegura que el error acumulado sea pequeño. Un defecto diminuto puede amplificarse durante miles de iteraciones, del mismo modo que una perturbación física puede crecer cerca de un equilibrio inestable. Necesitamos una noción separada de estabilidad numérica.

También importa el objeto al que aplicamos la malla. Si la difusión es heterogénea, la ecuación de conservación correcta suele adoptar la forma $u_{t} = \partial_{x}\left( D(x)u_{x} \right)$. Sustituirla por $D(x)u_{xx}$ omite el término $D\prime(x)u_{x}$, cuando existe. El error persiste aunque aumentemos indefinidamente la resolución, porque estamos resolviendo otra ecuación. Algo semejante ocurre al sustituir una función de flujo no lineal por su aproximación lineal fuera del régimen para el que fue obtenida.

A veces se llama «resolución» solamente al número de puntos por unidad de longitud. Sin embargo, la información relevante es qué estructuras espaciales y temporales interesa distinguir. Una interfaz estrecha, una capa límite o un frente de reacción puede requerir muchos puntos dentro de una región reducida, mientras que un patrón suave de gran longitud de onda puede tolerar una malla más gruesa. Por eso conviene relacionar $\Delta x$ con una escala característica justificada por el modelo y conservar esa justificación junto con los resultados.

149.3. Volúmenes finitos y prioridad del balance sobre la imagen

Cuando la magnitud es masa, partículas, energía u otra cantidad conservada, es especialmente útil comenzar por una ley integral. Para una densidad $u$, un flujo $J$ y una producción local $R$, el balance en una celda $\left\lbrack x_{j - 1/2},x_{j + 1/2} \right\rbrack$ es

Definamos ${\bar{U}}_{j}$ como el promedio de celda. Al aproximar los flujos en las caras obtenemos

Esta forma tiene una ventaja estructural, porque al sumar todos los balances cada flujo interior aparece dos veces con signos opuestos. Si $R = 0$, la variación de la cantidad total discreta depende únicamente de los flujos a través de los extremos del dominio. El resultado es algebraicamente exacto, salvo errores de redondeo. Una aproximación que viole sistemáticamente ese balance puede generar o destruir materia incluso cuando reproduce una curva aparentemente plausible.

Para la difusión de Fick, $J = - D\, u_{x}$, un flujo intercelular de primer orden en el gradiente es

Si $D$ es constante y la malla regular, el balance recupera en las celdas interiores el conocido laplaciano de tres puntos. Cuando $D$ cambia bruscamente entre materiales, el valor interfacial $D_{j + 1/2}$ necesita atención. En una interfaz donde deben ser continuos el flujo y, bajo condiciones apropiadas, la variable potencial, el promedio armónico suele respetar mejor el transporte a través de capas en serie que el promedio aritmético. Un promedio equivocado puede ser especialmente dañino cuando las difusividades de dos regiones son muy distintas.

La definición de la frontera debe mantenerse en términos físicos. Para imponer aislamiento basta fijar $J_{1/2} = J_{N + 1/2} = 0$ en la formulación de volúmenes finitos, con lo que se conserva la suma $\Delta x\sum_{j}^{}{\bar{U}}_{j}$. En un esquema nodal de diferencias finitas también es posible introducir puntos fantasma para representar derivadas normales, pero hay que cuidar las ponderaciones de los nodos extremos. Sumarlos como si todos representaran celdas completas puede producir una pérdida aparente de masa que en realidad procede de una interpretación inadecuada del esquema.

La conservación discreta es una buena prueba automática. Podemos iniciar un campo no negativo y heterogéneo, cerrar las fronteras, desactivar reacciones y ejecutar el código durante muchos pasos. Registramos la masa inicial y final, su deriva acumulada y el efecto de disminuir $\Delta t$. Si la deriva supera lo esperable por redondeo o por los errores del resolvedor lineal, conviene revisar primero las caras y los signos de los flujos. En sistemas con reacción, la prueba consiste en contrastar la variación total con la integral de $R$.

149.4. Difusión explícita y restricción sobre el paso temporal

Combinemos Euler hacia adelante con el laplaciano centrado en una dimensión. El esquema de difusión es

Reordenando,

Así escrito, cada paso es una combinación de los valores vecinos y del valor local. Para $0 \leq r \leq 1/2$, los coeficientes son no negativos y suman uno. En presencia de condiciones de frontera compatibles, el esquema conserva cotas máximas y mínimas, propiedad esperable de la difusión pura. Si $r > 1/2$, el peso del punto central se vuelve negativo; pueden surgir oscilaciones y amplificación de las componentes de longitud de onda corta.

La restricción clásica es

Es una condición de estabilidad propia de este esquema particular en una dimensión, ajena a la física de la ecuación de difusión. Si reducimos la distancia espacial a la mitad, el paso temporal admisible cae a la cuarta parte. Para simular durante el mismo intervalo necesitamos, por tanto, alrededor de cuatro veces más pasos; al mismo tiempo se duplica el número de puntos. En una dimensión el trabajo total de un algoritmo explícito simple puede multiplicarse aproximadamente por ocho; en dos o tres dimensiones el crecimiento es todavía más rápido.

Supongamos, por ejemplo, $D = 0,1\ m^{2}/s$ y $\Delta x = 0,05\ m$. El límite de estabilidad es $\Delta t \leq 0,0125\ s$. Elegir $\Delta t = 0,01\ s$ da $r = 0,4$, dentro del rango estable. Si mantuviéramos ese tiempo al pasar a $\Delta x = 0,025\ m$, obtendríamos $r = 1,6$, claramente inadmisible. El aspecto suave de las primeras iteraciones no invalidaría el diagnóstico, porque la inestabilidad puede comenzar como ruido de redondeo y terminar dominando toda la solución.

En una malla rectangular bidimensional, con difusividad constante, los pesos del punto central y de sus cuatro vecinos llevan a la condición $r_{x} + r_{y} \leq 1/2$, donde $r_{x} = D\Delta t/(\Delta x)^{2}$ y $r_{y} = D\Delta t/(\Delta y)^{2}$. Para $d$ dimensiones con espaciado común la condición es $D\Delta t/(\Delta x)^{2} \leq 1/(2d)$. Cuando aparecen reacciones acopladas, advección, mallas irregulares o coeficientes variables, el análisis debe rehacerse.

149.5. Análisis de von Neumann y comportamiento de cada longitud de onda

El razonamiento basado en pesos positivos da una condición clara para la difusión. El análisis espectral muestra además que cada modo espacial discreto se amplifica o atenúa con un factor propio. En una malla periódica, consideremos una perturbación de la forma

La sustitución en el esquema explícito da

Una perturbación inicialmente pequeña permanece acotada paso a paso si $\left| G(\theta) \right| \leq 1$ para los modos relevantes. El caso más restrictivo aparece en la escala más corta que la malla representa, $\theta = \pi$, para la cual $G = 1 - 4r$. Exigir que su módulo no supere uno conduce a $r \leq 1/2$. El resultado coincide con el argumento de convexidad y muestra que la restricción la impone el modo de alta frecuencia.

Una malla no distingue todas las longitudes de onda. Una función coseno muestreada sobre puntos separados por $\Delta x$ puede resultar indistinguible de otra con un número de onda diferente. Este fenómeno, denominado aliasing, no desaparece al avanzar con mayor precisión temporal. Para señales sobre una malla uniforme la frecuencia de Nyquist es $k_{N} = \pi/\Delta x$, pero estar formalmente por debajo de ella no asegura una representación precisa de las derivadas. Cerca de ese límite, la versión discreta de la curvatura difiere sustancialmente de $- k^{2}$.

En efecto, el laplaciano centrado aplicado a un modo de número de onda $k$ da

Para longitudes de onda grandes respecto de la malla, ${\widetilde{k}}^{\, 2} = k^{2} - k^{4}(\Delta x)^{2}/12 + O\left( (\Delta x)^{4} \right)$. El operador discreto subestima ligeramente el cuadrado del número de onda; la diferencia aumenta en modos menos resueltos. Si el crecimiento de un patrón depende de que $k$ caiga dentro de una banda estrecha de inestabilidad, sustituir $k^{2}$ por ${\widetilde{k}}^{\, 2}$ puede desplazar artificialmente los umbrales.

En medios con varias especies la perturbación de Fourier es un vector y el factor de amplificación se convierte en una matriz. Para un modelo de reacción–difusión linealizado cerca de un equilibrio, $\partial_{t}\mathbf{\eta} = J\mathbf{\eta} + D_{m}\partial_{xx}\mathbf{\eta}$, el modo de número de onda $k$ evoluciona mediante la matriz $J - k^{2}D_{m}$, donde $D_{m}$ es la matriz de difusividades. El esquema espacial aproxima esa matriz por $J - {\widetilde{k}}^{\, 2}D_{m}$ y el integrador temporal añade su propia modificación a las tasas de crecimiento. En el régimen cercano a una bifurcación, pequeñas diferencias pueden cambiar qué modos parecen inestables.

149.6. Métodos implícitos, estabilidad ampliada y nuevas exigencias

La limitación cuadrática del paso temporal puede evitarse mediante una aproximación implícita. Para la difusión lineal, Euler hacia atrás evalúa el laplaciano en el nuevo instante, de modo que

Los valores desconocidos de $U^{n + 1}$ están acoplados. Debemos resolver un sistema lineal $\left( I - \Delta tDL_{h} \right)U^{n + 1} = U^{n}$, donde $L_{h}$ es el laplaciano discreto con las fronteras incorporadas. En una dimensión el sistema es tridiagonal y hay algoritmos eficientes que explotan esa estructura. En dominios grandes o geometrías complejas se necesitan resolvedores iterativos, precondicionadores y un control explícito del error algebraico.

Para un modo con autovalor espacial $- {\widetilde{k}}^{\, 2}$, Euler implícito produce

El módulo de este factor nunca supera uno cuando $D \geq 0$ y $\Delta t > 0$, así que la difusión lineal se mantiene estable sin imponer $\Delta t \leq (\Delta x)^{2}/(2D)$. Cuando $D\Delta t{\widetilde{k}}^{\, 2}$ es grande, las componentes de alta frecuencia se atenúan con fuerza. Esa amortiguación puede ser ventajosa si los datos iniciales contienen estructuras mal resueltas. El precio es un sesgo temporal, porque la atenuación del método no coincide exactamente con $\exp\left( - D{\widetilde{k}}^{\, 2}\Delta t \right)$.

Una alternativa habitual es Crank–Nicolson, que promedia el operador de difusión entre los instantes inicial y final. Para cada modo su factor es

El método es de segundo orden temporal para soluciones regulares y estable en el sentido lineal habitual de la difusión. Sin embargo, si el producto $D\Delta t{\widetilde{k}}^{\, 2}$ es muy grande, el factor tiende a $- 1$ y el modo alterna de signo sin recibir la intensa amortiguación que debería producir la difusión verdadera. Pueden aparecer oscilaciones temporales de modos pequeños, especialmente con datos poco suaves. De aquí se desprende que la estabilidad en norma y la preservación de la positividad son propiedades distintas. Un método que no hace crecer la norma de un error puede crear valores negativos cuando la solución física debería permanecer no negativa.

Una fuente de incertidumbre que suele omitirse es que el sistema implícito rara vez se resuelve con exactitud matemática. Se usa una tolerancia de iteración. El residuo $\parallel b - AU^{n + 1} \parallel$ mide la calidad algebraica; si la matriz está mal condicionada, un residuo pequeño no garantiza un error de solución pequeño. Una tolerancia demasiado laxa puede impedir observar la convergencia espacial; una excesivamente estricta desperdicia recursos sin mejorar las conclusiones. La tolerancia debe ser pequeña respecto del error de discretización que pretendemos medir. En la documentación del cálculo conviene separar errores de tiempo, espacio y resolución algebraica.

149.7. Advección y difusión artificial cuando el método cambia la física

La advección describe el transporte de una magnitud por una velocidad. En una dimensión, para una velocidad constante $c > 0$, escribimos

Si las fronteras y el dominio lo permiten, la solución es la traslación $u(x,t) = u_{0}(x - ct)$; el perfil conserva su forma y se desplaza. Una discretización numérica inadecuada puede convertir ese movimiento en una serie de oscilaciones o transformar un frente agudo en una transición excesivamente suave. Ambas alteraciones pueden confundirse con propiedades emergentes del modelo.

La diferencia centrada combinada con Euler explícito da

Al sustituir un modo de Fourier obtenemos $G = 1 - iC\sin\theta$, de modo que $|G|^{2} = 1 + C^{2}\sin^{2}\theta$. Para cualquier $C \neq 0$ existen modos con amplificación superior a uno. Por eso el esquema FTCS centrado es inestable para advección pura con cualquier $\Delta t > 0$ en la malla considerada; reducir el paso sin cambiar el método no elimina el mecanismo de amplificación.

Con $c > 0$, la diferencia a barlovento aproxima $u_{x}$ con los valores del punto y de su vecino izquierdo. El esquema resultante es

Si $0 \leq C \leq 1$, el nuevo valor es una combinación convexa y el método evita aquella inestabilidad. El vecino elegido depende de la dirección del transporte; para $c < 0$ corresponde usar el derecho. Este esquema, sin embargo, suaviza los perfiles incluso cuando la ecuación exacta no incluye difusión. El análisis de ecuación modificada para una solución suficientemente suave muestra que su error principal se comporta como un término difusivo con coeficiente aproximado

Así, si un investigador observa cómo se ensancha una mancha transportada, parte de ese ensanchamiento podría proceder del método. Con $C = 1$ el esquema se reduce a desplazar un punto de malla por paso y, en un dominio periódico, transporta exactamente los valores nodales; esa exactitud no se extiende a cualquier velocidad variable o geometría.

En la ecuación de advección–difusión, $u_{t} + cu_{x} = Du_{xx}$, aparece el número de Péclet de malla $Pe_{h} = |c|\Delta x/D$. Cuando es grande, el transporte domina sobre la difusión física a escala de una celda y el tratamiento de flujos convectivos se vuelve especialmente delicado. Métodos de orden superior pueden reducir la disipación numérica, pero necesitan limitadores u otras técnicas para evitar oscilaciones cerca de discontinuidades. La expresión «alta resolución» designa aquí una familia de compromisos matemáticos, sin garantía automática de fidelidad.

149.8. Modos de una frontera finita y selección de patrones por la malla

Los métodos numéricos interactúan con la geometría del dominio. Para la ecuación lineal de difusión, las condiciones periódicas en un intervalo de longitud $L$ permiten modos $\exp(i2\pi mx/L)$, con $m$ entero. En cambio, con derivada normal cero en ambos extremos aparecen modos coseno con $k_{n} = n\pi/L$, $n = 0,1,2,\ldots$. Si se imponen valores nulos de la perturbación, surgen modos seno con $k_{n} = n\pi/L$, $n \geq 1$. La diferencia afecta a las escalas espaciales que pueden existir, incluso antes de introducir una malla computacional.

Para una dinámica cuya tasa de crecimiento depende de $k$, es útil representar el conjunto $\lambda(k)$. Si la tasa es positiva en un intervalo estrecho de números de onda, el dominio puede admitir uno de esos modos, varios o ninguno. Las fronteras fijan una lista discreta de candidatos y la ecuación establece cómo crece cada uno. La malla finita introduce luego otra modificación, porque el laplaciano centrado sustituye $k^{2}$ por su contraparte discreta ${\widetilde{k}}^{\, 2}$. Si el máximo de $\lambda$ está cerca del borde de una banda inestable, una malla basta para alterar la conclusión.

En los modelos de formación de patrones, un diagrama de manchas que aparece en una simulación con fronteras reflectantes y desaparece con fronteras periódicas podría ser un efecto genuino del espectro del dominio. También podría deberse a un cambio inadvertido en el código, a una malla demasiado gruesa o a condiciones de frontera implementadas con signos incorrectos. Sin un análisis modal o una prueba de refinamiento no tenemos manera de separar estas explicaciones.

El cambio de frontera puede modificar además las cantidades conservadas y las restricciones topológicas. En una difusión pura, un dominio periódico o cerrado conserva el promedio; con bordes absorbentes se pierde masa. En un medio excitable, imponer una frontera puede absorber o reflejar perturbaciones según el modelo; ciertas trayectorias pueden quedar ancladas a defectos. La geometría afecta también al comportamiento de ondas y espirales descrito en C145–C148. Por ello, la condición de frontera es un componente del modelo y debe consignarse junto con sus ecuaciones.

149.9. Ejemplo resuelto con una onda de calor calculada y contrastada

Tomemos un intervalo de longitud $L = 1\ m$, con $0 \leq x \leq L$, difusividad $D = 0,1\ m^{2}/s$, fronteras con valor nulo y condición inicial $u(x,0) = \sin(\pi x/L)$. La solución exacta es

El perfil no cambia de forma; su amplitud disminuye de manera exponencial. Este caso es una prueba de código especialmente exigente porque podemos medir el error sin aproximaciones Monte Carlo ni incertidumbre observacional. Una simulación que falle aquí difícilmente merece confianza cuando se la aplique a un sistema no lineal sin solución cerrada.

Elegimos $N = 20$ (de modo que $\Delta x = 0,05$) y $\Delta t = 0,01\ s$. La razón $r = D\Delta t/(\Delta x)^{2}$ vale $0,4$, así que el esquema explícito cumple su condición de estabilidad. El modo inicial seno es también un autovector del laplaciano discreto con condiciones de Dirichlet nodales. Cada paso multiplica su amplitud por

Tras diez pasos, en $t = 0,1\ s$, la amplitud calculada es $G^{10}$ y la exacta es $e^{- 0,1\pi^{2} \cdot 0,1}$. Ambas magnitudes pueden calcularse directamente; su diferencia absoluta es un error comprobable de la implementación. El signo y el tamaño del error dependen de los términos temporales y espaciales de truncamiento; no debemos atribuirlo solo a la malla espacial.

Con los parámetros indicados se obtiene $G \approx 0,9901506725$. Diez iteraciones producen una amplitud $G^{10} \approx 0,9057594371$, mientras que la solución exacta en el mismo instante tiene amplitud $0,9060180558$. El error absoluto del factor de amplitud es aproximadamente $2,5862 \times 10^{- 4}$. Repetir el cálculo con $N = 40$ y $N = 80$, conservando $r = 0,4$, produce respectivamente errores cercanos a $6,4427 \times 10^{- 5}$ y $1,6093 \times 10^{- 5}$. Los cocientes de errores sucesivos están próximos a cuatro, de acuerdo con la convergencia cuadrática esperada al ligar el paso temporal al cuadrado del espacial. Estos valores se han recalculado directamente con la fórmula del factor de amplificación, sin ejecutar aún un algoritmo general de difusión.

Para construir un estudio de convergencia podemos repetir la prueba con $N = 40$ y $N = 80$, conservando $r = 0,4$. Ello exige disminuir $\Delta t$ por un factor cuatro cada vez que se duplica $N$ y elegir un tiempo final que sea múltiplo de todos los pasos. Al usar $\Delta t \propto (\Delta x)^{2}$, tanto el defecto espacial de segundo orden como el error global temporal de primer orden pasan a ser del orden de $(\Delta x)^{2}$. Se espera entonces, en el régimen asintótico y para esta solución suave, que el error se reduzca aproximadamente cuatro veces por refinamiento. Esta predicción debe confrontarse con resultados medidos.

En esta prueba fijamos $U_{0}^{n} = U_{N}^{n} = 0$ en todo instante, valor que coincide con la condición exacta. La ecuación de actualización se aplica a los nodos interiores $j = 1,\ldots,N - 1$. Si alguien actualizara los extremos mediante una fórmula periódica, la simulación podría seguir produciendo un campo suave, pero correspondería a otro problema. Conservar las condiciones de frontera en el código, en el informe y en las pruebas es parte de reproducir una solución matemática.

149.10. Protocolo de convergencia y criterios para aceptar una simulación

Una simulación de ecuaciones espaciales debe justificarse mediante una secuencia de pruebas que respondan a fallos distintos. El primer nivel verifica propiedades algebraicas elementales (tamaños de matrices, índices de vecinos, conservación cuando corresponde, positividad cuando el esquema la garantiza y ausencia de errores al aplicar las fronteras). El segundo compara el algoritmo con soluciones de referencia. El tercero explora cómo cambian los observables al modificar la resolución y la estrategia temporal.

Para aislar fuentes de error conviene hacer dos refinamientos complementarios. En el primero se fija una malla espacial suficientemente precisa y se disminuye solo el paso temporal; en el segundo se escoge un integrador temporal cuyo error sea despreciable frente al espacial y se refina solo la malla. Los esquemas explícitos pueden dificultar este aislamiento porque sus restricciones de estabilidad ligan ambos pasos, pero los esquemas implícitos o los casos con soluciones modales conocidas permiten organizar la comparación. El resultado deseable es una región donde el observable apenas cambia al continuar refinando y donde se entiende el orden de los términos residuales. Un único contraste entre dos mallas deja abierta la posibilidad de una coincidencia accidental. En problemas complejos estas pruebas requieren tiempo, pero sin ellas la interpretación científica depende de una elección de programación.

Para estudiar la convergencia, fijamos un observable $Q$ y calculamos $Q_{h}$, $Q_{h/2}$ y $Q_{h/4}$ con pasos temporales compatibles. Si el error dominante se comporta como $Q_{h} = Q_{*} + Ch^{p} + \cdots$, las diferencias entre resoluciones sucesivas permiten estimar un orden observado,

La fórmula presupone que las tres mallas están en el régimen asintótico, que el coeficiente dominante no se anula y que los denominadores no son demasiado pequeños. Si los errores temporal, espacial y algebraico se compensan, el cociente puede adoptar valores engañosos. Un aparente orden elevado en un único observable tampoco asegura que otras propiedades, como una probabilidad de escape o la posición de una interfaz, sean igual de precisas.

Cuando existe una solución exacta conviene usarla. Si no la hay, puede compararse con una solución de referencia calculada mediante otra discretización, una malla mucho más fina o un integrador independiente; en ese caso hay que declarar que la referencia también tiene error. La extrapolación de Richardson puede estimar y reducir el error dominante en condiciones favorables, pero pierde fiabilidad junto a discontinuidades, cambios de régimen o bifurcaciones. Cerca de un umbral crítico, el signo de una tasa de crecimiento diminuta puede ser más importante que la distancia global entre dos curvas.

En sistemas no lineales conviene observar también invariantes, balances y escalas características. Puede ocurrir que una solución sea estable en todas las simulaciones y, aun así, converja a un estado espurio por mala elección de frontera. También es posible que dos mallas produzcan imágenes similares pero difieran notablemente en el período de una onda, la masa total o la velocidad de un frente. Un informe responsable identifica qué observables están numéricamente resueltos, cuáles son sensibles y cuáles siguen siendo exploratorios.

Las simulaciones caóticas añaden otra dificultad. Dos trayectorias próximas pueden separarse rápidamente por sensibilidad a condiciones iniciales o a errores de redondeo. En horizontes largos la coincidencia punto a punto entre trayectorias deja de ser una exigencia razonable. En ese contexto se examinan propiedades estadísticas, espectros, distribuciones y balances, con un tiempo de observación suficiente y medidas que no estén dominadas por transitorios.

Para dejar un cálculo auditable hay que registrar el problema continuo, el dominio y sus unidades, las condiciones iniciales y de frontera, la forma de los flujos, la malla, el integrador temporal, sus parámetros, las tolerancias del resolvedor y las cantidades usadas para juzgar la convergencia. También se guardan las versiones del programa y, cuando interviene aleatoriedad, la estrategia de semillas. La figura final debería surgir de esos datos y permitir reconstruir sus valores. Poder repetir una imagen todavía no demuestra su fidelidad matemática, pero es la condición básica para examinarla.

En términos generales, una aproximación consistente y estable converge bajo las hipótesis adecuadas del problema lineal bien planteado al que se aplica; en modelos no lineales y situaciones singulares hacen falta argumentos adicionales. El teorema de equivalencia de Lax–Richtmyer precisa la relación entre consistencia, estabilidad y convergencia para problemas lineales de valor inicial bien planteados en un marco apropiado. No es una garantía universal para cualquier código que resuelva una EDP. En el capítulo C150 aplicaremos esta disciplina a un caso donde las fronteras alteran de modo demostrable la disponibilidad de un patrón.

En una región acotada, la estabilidad de la fórmula interior no garantiza la del esquema completo. Svärd y Nordström (2014) revisan métodos de suma por partes y tratamientos de frontera que permiten establecer estimaciones discretas de energía. Su análisis aclara por qué la discretización de contorno es parte del método y debe someterse a una prueba propia.

Referencias

• Courant, R., Friedrichs, K. O., y Lewy, H. (1928). Über die partiellen Differenzengleichungen der mathematischen Physik. Mathematische Annalen, 100, 32–74.
• Lax, P. D., y Richtmyer, R. D. (1956). Survey of the stability of linear finite difference equations. Communications on Pure and Applied Mathematics, 9(2), 267–293. https://doi.org/10.1002/cpa.3160090206
• LeVeque, R. J. (2007). Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems. SIAM. https://doi.org/10.1137/1.9780898717839
• Strikwerda, J. C. (2004). Finite Difference Schemes and Partial Differential Equations (2.ª ed.). SIAM.
• Trefethen, L. N. (2000). Spectral Methods in MATLAB. SIAM.
• Cross, M. C., y Hohenberg, P. C. (1993). Pattern formation outside of equilibrium. Reviews of Modern Physics, 65, 851–1112. https://doi.org/10.1103/RevModPhys.65.851
• Svärd, M., y Nordström, J. (2014). Review of summation-by-parts schemes for initial–boundary-value problems. Journal of Computational Physics, 268, 17–38. https://doi.org/10.1016/j.jcp.2014.02.031