Integración numérica y control del error

Una computadora no dibuja directamente la solución continua: construye una sucesión de aproximaciones. En dinámica no lineal, el método y el paso pueden introducir amortiguación, crecimiento o desplazamientos de fase que no pertenecen al modelo. Por eso forman parte del experimento y deben quedar junto a la ecuación, no escondidos como una preferencia del software.

Para $\dot{\mathbf{x}}=\mathbf{f}(t,\mathbf{x})$, Euler avanza mediante

$$\mathbf{x}_{n+1}=\mathbf{x}_n+h\mathbf{f}(t_n,\mathbf{x}_n).$$

Su error local de truncamiento es $O(h^2)$ y el acumulado en un intervalo fijo es $O(h)$. Heun usa una predicción

$$\mathbf x_{n+1}^{P}=\mathbf x_n+h\mathbf f(t_n,\mathbf x_n)$$

y promedia la pendiente inicial con la evaluada al final:

$$\mathbf x_{n+1}=\mathbf x_n+\frac h2 \left[\mathbf f(t_n,\mathbf x_n)+ \mathbf f(t_{n+1},\mathbf x_{n+1}^{P})\right].$$

Así alcanza orden global dos. RK4 combina cuatro evaluaciones:

$$\mathbf{x}_{n+1}=\mathbf{x}_n+\frac{h}{6}(\mathbf{k}_1+2\mathbf{k}_2+2\mathbf{k}_3+\mathbf{k}_4),$$

con error global $O(h^4)$ bajo hipótesis regulares. Un orden alto se acompaña de un paso compatible con la estabilidad del método.

Las constantes ocultas en $O(h^p)$ importan. Un método de orden mayor no garantiza menor error con cualquier paso, y cuatro evaluaciones de $f$ cuestan más que una. La comparación debe fijar si se mide error por paso, por evaluación o por tiempo de cómputo.

La ecuación $y'=-2y$ permite aislar el error del algoritmo porque su solución exacta es conocida. Con el mismo paso, Euler, Heun y RK4 producen aproximaciones distintas; el panel de error revela diferencias que la superposición de curvas puede ocultar. La comparación es justa porque usa la misma malla y la misma condición inicial.

Refinar $h$ convierte esa comparación puntual en una prueba de orden. En ejes logarítmicos, una pendiente uno, dos o cuatro corresponde al comportamiento global esperado de Euler, Heun o RK4. La región útil es aquella donde las curvas siguen una ley aproximadamente recta; fuera de ella pueden dominar inestabilidad, redondeo o una malla todavía demasiado gruesa.

Error local, error global y estabilidad

El error local pregunta cuánto se desvía un solo paso si empieza desde el valor exacto. El error global incorpora la propagación de todos los pasos hasta un tiempo fijo. La estabilidad absoluta pregunta si el método reproduce la contracción de la ecuación de prueba

$$y'=\lambda y, \qquad \operatorname{Re}\lambda<0.$$

Euler produce $y_{n+1}=R(z)y_n$ con $R(z)=1+z$ y $z=h\lambda$. La aproximación se contrae sólo si $|1+z|<1$. Para $\lambda$ real negativa esto exige $0<h<2/|\lambda|$. La solución exacta es estable para cualquier $h$ porque no usa una malla; la restricción pertenece al algoritmo.

RK4 tiene función de estabilidad

$$R(z)=1+z+\frac{z^2}{2}+\frac{z^3}{6}+\frac{z^4}{24}.$$

Su región es mayor, pero tampoco cubre todo el semiplano izquierdo. Los problemas rígidos contienen escalas muy distintas y pueden obligar a pasos diminutos en métodos explícitos aunque el comportamiento de interés sea lento. Alligood, Sauer y Yorke resumen solución de EDO, error y adaptación dentro del instrumental computacional para dinámica. Alligood, 1996 Kathleen T. Alligood, Tim D. Sauer y James A. Yorke (1996) Chaos: An Introduction to Dynamical Systems Springer Ubicación consultada: apéndice B, pp. 567–576

Método de Euler en el decaimiento exponencial

Para $\dot x=-2x$, $x(0)=1$, la solución exacta es $e^{-2t}$. Euler produce $x_n=(1-2h)^n$. Si $h=0.1$, a $t=1$ hay diez pasos y

$$x_{10}=0.8^{10}\approx0.1074, \qquad e^{-2}\approx0.1353.$$

Con $h=0.05$, $x_{20}=0.9^{20}\approx0.1216$, más cerca del valor exacto. Si $h>1$, el factor $1-2h$ tiene magnitud mayor que uno y el método diverge aunque la solución real se contraiga. El fallo es numérico.

Para $h=0.75$, el factor es $-0.5$: Euler converge, pero alterna de signo aunque la solución exacta permanece positiva. Para $h=1$, el factor es $-1$ y aparece una oscilación numérica de amplitud constante. Para $h=1.1$, $|1-2h|=1.2$ y la amplitud crece. Esta secuencia muestra por qué «no explota» es un criterio demasiado débil: un paso puede ser estable y todavía deformar signo y fase.

El error en $t=1$ puede estimarse expandiendo

$$(1-2h)^{1/h}=\exp\left(\frac{\log(1-2h)}{h}\right) =e^{-2}\left(1-2h+O(h^2)\right).$$

La diferencia principal es proporcional a $h$, en acuerdo con el orden global uno.

Orden observado mediante refinamiento

La prueba de convergencia incorpora dos curvas en una ventana corta y la norma $\|\mathbf{x}_{h}(t)-\mathbf{x}_{h/2}(t)\|$. En un régimen sensible, las trayectorias puntuales terminarán separándose; aun así pueden coincidir en estadísticas, geometría y escalas temporales. Distingue la fidelidad de corto plazo del acuerdo cualitativo de largo plazo.

Si no hay solución exacta, se calculan aproximaciones con $h$, $h/2$ y $h/4$ sobre instantes comunes. Para un método de orden $p$ en su régimen asintótico,

$$Q_h=\frac{\|x_h-x_{h/2}\|}{\|x_{h/2}-x_{h/4}\|} \approx2^p, \qquad p_{\mathrm{obs}}=\log_2Q_h.$$

La norma, el intervalo y las variables escaladas deben declararse. Si $p_{\mathrm{obs}}$ no se estabiliza, puede faltar refinamiento, dominar redondeo o existir una discontinuidad que degrade el orden. Hirsch, Smale y Devaney colocan los métodos numéricos junto a existencia y dependencia de datos: la convergencia del algoritmo presupone que el problema continuo está bien planteado. Hirsch, 2013 Morris W. Hirsch, Stephen Smale y Robert L. Devaney (2013) Differential Equations, Dynamical Systems, and an Introduction to Chaos 3.ª ed. · Academic Press Ubicación consultada: cap. 7, pp. 139–158

Convergencia de corto plazo en Lorenz

Selecciona Lorenz, $\sigma=10$, $\rho=28$, $\beta=8/3$, $(1,1,1)$ y tiempo $40$. Ejecuta RK4 con $dt=0.02$, $0.01$ y $0.005$. Si el selector principal ofrece Euler o Heun, repite $dt=0.01$ con cada uno; si no, usa el comparador numérico del Explorador Sprott para el mismo propósito.

Compara las series durante $0\le t\le5$ y luego las proyecciones de toda la ventana. Anota el instante aproximado en que dos trazas dejan de coincidir visualmente. La selección del método usa refinamiento, acotación y costo por paso, además de la suavidad visual.

No uses la separación tardía como prueba de que una corrida es errónea. En un régimen sensible, dos errores iniciales pequeños crecen y las trayectorias dejan de coincidir aunque ambos cálculos resuelvan adecuadamente el flujo durante un horizonte finito. Primero mide error en una ventana corta; después compara propiedades de largo plazo como acotación, distribución espacial, espectro o promedios, siempre con varias mallas.

El transitorio tampoco debe elegirse porque «la curva ya se ve bonita». Define una cantidad, por ejemplo energía, distancia a una sección o media móvil, y comprueba que su estadística se estabilice al aumentar el descarte. Si el sistema posee transitorios largos o caos transitorio, la ventana necesaria puede depender mucho de la condición inicial.

Paso adaptativo y muestreo de salida

Un integrador adaptativo estima error local comparando aproximaciones de distinto orden y ajusta $h$ para mantenerlo cerca de una tolerancia. La malla interna deja de ser uniforme. Para Fourier, retornos o comparación entre señales, se necesita una salida densa evaluada en tiempos uniformes; seleccionar los puntos internos del integrador sesga el muestreo hacia regiones donde el algoritmo tomó pasos cortos.

Las tolerancias relativa y absoluta cumplen funciones distintas. Una regla típica controla cada componente mediante

$$\mathrm{escala}_i=\mathrm{atol}_i+ \mathrm{rtol}\max(|x_i^{(n)}|,|x_i^{(n+1)}|).$$

Si las variables tienen magnitudes muy diferentes, una sola tolerancia absoluta puede ignorar las pequeñas o sobrerresolver las grandes. Escalar el modelo o usar tolerancias por componente forma parte del análisis. Fuchs ofrece un resumen breve de estos instrumentos numéricos junto a sus ejemplos dinámicos. Fuchs, 2013 Armin Fuchs (2013) Nonlinear Dynamics in Complex Systems: Theory and Applications for the Life-, Neuro- and Natural Sciences Springer Ubicación consultada: apéndice B, pp. 205–212

Pruebas de integración mediante refinamiento

El objetivo de un integrador no es producir la figura esperada, sino aproximar el problema continuo con un error auditable. Orden, estabilidad, tolerancias y muestreo convierten una trayectoria calculada en evidencia numérica utilizable.