Métodos de Runge–Kutta: consistencia, estabilidad y convergencia

Un método consistente reproduce la ecuación cuando el paso tiende a cero. La estabilidad controla la propagación de ciertos errores. Un método convergente aproxima la solución en un intervalo fijo al refinar la malla. Las tres palabras responden preguntas distintas y deben acompañarse de la propiedad que se ha comprobado.

Identidad integral y métodos de un paso

El problema

$$y'=f(t,y),\qquad y(t_n)=y_n$$

satisface formalmente

$$y(t_n+h)=y(t_n)+\int_{t_n}^{t_n+h}f(s,y(s))\,ds.$$

Los métodos de Runge–Kutta aproximan esta integral mediante evaluaciones del campo dentro del paso. Euler usa la pendiente izquierda. Heun predice con Euler y promedia las pendientes de ambos extremos:

$$k_1=f(t_n,y_n),\qquad k_2=f(t_n+h,y_n+hk_1),$$ $$y_{n+1}=y_n+\frac h2(k_1+k_2).$$

El punto medio explícito también tiene orden dos, pero evalúa

$$k_1=f(t_n,y_n),\qquad k_2=f(t_n+h/2,y_n+hk_1/2), \qquad y_{n+1}=y_n+hk_2.$$

Ambos comparten orden y pueden tener constantes de error distintas. Decir «Runge–Kutta de orden dos» no identifica un algoritmo único.

RK4 clásico usa cuatro etapas:

$$k_1=f(t_n,y_n),$$ $$k_2=f(t_n+h/2,y_n+hk_1/2),\qquad k_3=f(t_n+h/2,y_n+hk_2/2),$$ $$k_4=f(t_n+h,y_n+hk_3),$$ $$y_{n+1}=y_n+\frac h6(k_1+2k_2+2k_3+k_4).$$

Los pesos se eligen para igualar la expansión de Taylor hasta orden cuatro sin calcular derivadas de $f$ explícitamente. Fuchs presenta estos esquemas como herramientas para modelos dinámicos y vincula las expansiones de Taylor con EDO no lineales. 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

Para la ecuación de prueba $y'=\lambda y$, Euler produce

$$y_{n+1}=(1+h\lambda)y_n.$$

Si $\operatorname{Re}\lambda<0$, la solución exacta decae. Euler reproduce ese decaimiento sólo cuando $|1+h\lambda|<1$. Por eso un sistema continuo estable puede generar una sucesión numérica creciente si el paso está fuera de la región de estabilidad.

Cada método de un paso aplicado a la ecuación de prueba tiene la forma $y_{n+1}=R(z)y_n$, con $z=h\lambda$. Para Euler,

$$R_E(z)=1+z;$$

para Heun,

$$R_H(z)=1+z+\frac{z^2}{2};$$

y para RK4,

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

La región de estabilidad absoluta es $\{z:|R(z)|<1\}$. Aumentar el orden agranda ciertas porciones de la región, pero no vuelve a RK4 estable para todo el semiplano izquierdo. Un autovalor muy negativo puede imponer un paso diminuto aunque su modo desaparezca rápidamente: ésa es una de las señales de rigidez.

Consistencia local y acumulación global

El defecto local de un método se obtiene insertando la solución exacta en una actualización. Para Euler, Taylor da

$$y(t_{n+1})=y(t_n)+hy'(t_n)+\frac{h^2}{2}y''(\xi_n),$$

así que el error local es $O(h^2)$. En un intervalo fijo hay $N\approx(T-t_0)/h$ pasos; bajo estabilidad y regularidad, la acumulación produce error global $O(h)$. En general, un método de orden global $p$ tiene error local $O(h^{p+1})$.

La palabra «estabilidad» aquí no describe el equilibrio del sistema físico. Describe cómo el esquema propaga perturbaciones o errores dentro de una clase de problemas. Un método puede ser estable para cierto $h$ y todavía demasiado impreciso; también puede ser de orden alto y operar fuera de su región estable.

Hirsch, Smale y Devaney tratan los métodos numéricos junto con existencia y dependencia continua. La conexión importa: convergencia supone que el problema inicial define una solución suficientemente regular y que el método aproxima esa solución, no sólo una curva suave. 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

Una prueba de refinamiento conserva ecuación, condición inicial e intervalo, y repite con $h$, $h/2$, $h/4$. Con solución exacta, calcula $E(h)=\|x_h(T)-x(T)\|$. Cuando falta una solución exacta, usa una referencia más fina o diferencias sucesivas y declara esa condición. El orden observado es

$$p_{\mathrm{obs}}\approx \frac{\log(E(h)/E(h/2))}{\log 2}.$$

Con tres mallas y sin solución exacta puede estimarse

$$p_{\mathrm{obs}}\approx \log_2\left( \frac{\|y_h-y_{h/2}\|}{\|y_{h/2}-y_{h/4}\|} \right),$$

si las soluciones se comparan en tiempos comunes y ya se encuentran en el régimen asintótico. Una única razón no demuestra el orden. Se buscan varios refinamientos con pendiente estable; al final pueden dominar redondeo, tolerancias internas o error de la referencia.

El error en una trayectoria caótica crece por discretización y sensibilidad dinámica. Una prueba de orden usa tiempos suficientemente cortos para observar el régimen numérico. Para tiempos largos, fundamenta la comparación en propiedades estadísticas, secciones o distribuciones.

Tres métodos en una ecuación lineal

Para $y'=-2y$, $y(0)=1$ y $h=0.1$, Euler multiplica por $R_E(-0.2)=0.8$, Heun por

$$R_H(-0.2)=1-0.2+0.02=0.82,$$

y RK4 por

$$R_4(-0.2)=1-0.2+0.02-0.001333\ldots+0.0000667\ldots \approx0.818733.$$

El factor exacto por paso es $e^{-0.2}\approx0.818731$. RK4 casi coincide; Heun queda cerca y Euler sobre-amortigua. A $t=1$, cada factor se eleva a la décima potencia, por lo que un pequeño error de paso se acumula.

Ahora tome $h=1.1$. Euler usa $R_E(-2.2)=-1.2$ y crece alternando signos. La solución exacta sigue decayendo. El fallo se predice antes de ejecutar al ubicar $z=-2.2$ fuera del disco de Euler. Este ejemplo separa precisión, signo y estabilidad en una ecuación cuya solución no plantea ambigüedad.

Implementación explícita de las etapas

Antes de llamar a una biblioteca, implementa Euler, Heun y RK4 con la misma interfaz step(f, t, y, h). Prueba que no muten el estado de entrada, que acepten vectores y que evalúen $f$ el número esperado de veces. Para una integración completa, construye la malla de forma que el último paso termine exactamente en $T$ o declara el intervalo realmente alcanzado.

La documentación oficial de integración de SciPy reúne solve_ivp, cuadratura y herramientas relacionadas. Sus métodos adaptativos no son equivalentes a repetir RK4 fijo: usan pares embebidos, interpolación densa y control de error. El experimento de paso fijo sirve como referencia para auditar esas rutinas.

Alligood, Sauer y Yorke incluyen integradores y control adaptativo dentro de su apéndice computacional. La fuente justifica el uso en dinámica y los ejercicios de refinamiento; no se le atribuye una cobertura completa de álgebra numérica o cuadratura. 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

Error de fase y error de amplitud

En un oscilador, medir sólo $\|y_h(T)-y(T)\|$ mezcla amplitud y fase. Una frecuencia numérica ligeramente distinta produce un error grande después de muchas vueltas aunque el radio permanezca casi correcto. Por eso se registran por separado los tiempos de cruce, la envolvente y una cantidad de energía.

Para el oscilador armónico $x'=v$, $v'=-x$, Euler actualiza

$$\begin{pmatrix}x_{n+1}\\v_{n+1}\end{pmatrix} = \begin{pmatrix}1&h\\-h&1\end{pmatrix} \begin{pmatrix}x_n\\v_n\end{pmatrix}.$$

Los multiplicadores tienen módulo $\sqrt{1+h^2}>1$, de modo que la energía crece en cada paso. La figura puede parecer una órbita durante pocas vueltas, pero el análisis de la matriz predice una espiral hacia fuera para cualquier $h>0$. Euler simpléctico modifica el orden de las actualizaciones y mantiene la energía acotada con oscilación; no es el mismo método aunque use una evaluación por componente.

Esta comparación introduce una idea importante: dos algoritmos del mismo orden pueden preservar estructuras distintas. En problemas hamiltonianos, un método simpléctico puede ser preferible para tiempos largos aunque un RK de orden alto tenga menor error puntual al inicio. En problemas disipativos se vigila que el método no añada o elimine amortiguación en una escala comparable con la física.

Costo por evaluaciones de la función

Euler usa una evaluación de $f$ por paso, Heun dos y RK4 cuatro. Comparar todos con el mismo $h$ responde qué gana el orden a una malla fija; comparar con el mismo número de evaluaciones responde qué gana por costo de campo. Ambas preguntas son válidas y deben distinguirse.

Si evaluar $f$ implica resolver un sistema, consultar una tabla o ejecutar un modelo grande, el costo de administración del integrador es secundario. Si $f$ es barato y el estado es enorme, almacenamiento e interpolación pueden dominar. Una tabla útil incluye evaluaciones, tiempo, memoria, error final, error máximo y propiedades específicas como fase o energía.

El orden observado sólo se interpreta donde el error se comporta como $Ch^p$. Con pasos gruesos puede dominar la inestabilidad o un término de orden superior; con pasos diminutos aparecen redondeo y error de referencia. Seleccionar a posteriori dos puntos que produzcan la pendiente deseada no valida el método: la región debe sostenerse en varios refinamientos contiguos y bajo una norma declarada.

Prueba empírica del orden de convergencia

Consistencia conecta el método con la ecuación; estabilidad limita la propagación del error; convergencia se observa al refinar sobre un intervalo fijo. Una tabla con error, costo y función de estabilidad permite distinguir esas tres afirmaciones.