Euler, Runge-Kutta y simulación dinámica

Integración temporal de ecuaciones de movimiento mediante esquemas de Euler, Verlet y Runge-Kutta RK4.

1. Planteamiento de Sistemas Dinámicos y Ecuaciones Diferenciales Ordinarias

En la física teórica y computacional, la evolución temporal del estado de un sistema clásico, cuántico o termodinámico fuera del equilibrio se formula comúnmente mediante un sistema de ecuaciones diferenciales ordinarias (EDO) de primer orden de la forma:

$\frac{d\vec{y}}{dt} = \vec{f}(t, \vec{y}), \quad \vec{y}(t_0) = \vec{y}_0$

donde $\vec{y}(t) \in \mathbb{R}^N$ representa el vector de variables de estado (como posiciones $\vec{r}$, velocidades $\vec{v}$, cargas o corrientes) y $\vec{f}: \mathbb{R} \times \mathbb{R}^N \to \mathbb{R}^N$ es el campo vectorial que define las ecuaciones de movimiento. Cualquier ecuación diferencial de orden superior, como la segunda ley de Newton $m \frac{d^2 x}{dt^2} = F(t, x, \dot{x})$, se transforma de manera canónica en un sistema de primer orden mediante la introducción de variables auxiliares:

$\vec{y} = \begin{pmatrix} x \\ v \end{pmatrix}, \quad \frac{d\vec{y}}{dt} = \begin{pmatrix} v \\ \frac{1}{m} F(t, x, v) \end{pmatrix}$

El objetivo de la integración numérica temporal es construir una sucesión discreta de estados $\{\vec{y}_n\}$ que aproxime la solución exacta $\vec{y}(t_n)$ en tiempos discretizados $t_n = t_0 + n h$, donde $h = \Delta t$ es el paso de integración.

2. El Método de Euler Explícito y sus Límites en Dinámica

El método de Euler hacia adelante (Forward Euler) surge de truncar la serie de Taylor de $\vec{y}(t+h)$ en el término de primer orden:

$\vec{y}(t+h) = \vec{y}(t) + h \frac{d\vec{y}}{dt} + \frac{h^2}{2} \frac{d^2\vec{y}}{dt^2} + \mathcal{O}(h^3)$

Eliminando los términos de segundo orden y superiores, se obtiene el esquema iterativo:

$\vec{y}_{n+1} = \vec{y}_n + h \vec{f}(t_n, \vec{y}_n)$

Propiedades y Limitaciones Críticas:

1. Error Local de Truncamiento (LTE): En cada paso individual, el error introducido es de orden $\mathcal{O}(h^2)$.
2. Error Global de Truncamiento (GTE): Tras integrar un intervalo finito $T = N h$, la acumulación de errores produce un error global proporcional a $\mathcal{O}(h)$. Por tanto, Euler es un método de primer orden.
3. No Preservación del Espacio de Fases (Falta de Simpléctica): Para un oscilador armónico simple $\dot{x} = v$, $\dot{v} = -\omega^2 x$, la matriz de evolución por paso es: $\begin{pmatrix} x_{n+1} \\ v_{n+1} \end{pmatrix} = \begin{pmatrix} 1 & h \\ -\omega^2 h & 1 \end{pmatrix} \begin{pmatrix} x_n \\ v_n \end{pmatrix}$ El determinante de la transformación es $\det(\mathbf{M}) = 1 + \omega^2 h^2 > 1$. Dado que el determinante excede la unidad, el volumen en el espacio de fases se expande artificialmente en cada paso, inyectando energía no física al sistema sin límite.

3. Integración Simpléctica: Euler Semicomplícito y Verlet

Para sistemas hamiltonianos conservativos donde $\mathcal{H} = T(\vec{p}) + V(\vec{q})$, los métodos simplécticos garantizan la conservación del volumen del espacio de fases (Teorema de Liouville) y preservan un hamiltoniano sombra modificado a lo largo de tiempos de simulación muy extensos.

Euler Simpléctico (Euler-Cromer)

Actualiza primero la velocidad y utiliza el nuevo valor de velocidad inmediatamente para actualizar la posición:

$v_{n+1} = v_n + h \, a(t_n, x_n)$ $x_{n+1} = x_n + h \, v_{n+1}$

La matriz de transferencia tiene determinante exactamente igual a 1 ($\det(\mathbf{M}) = 1$), lo que estabiliza las órbitas y acota las oscilaciones de energía sin deriva secular divergente.

Algoritmo de Verlet de Velocidades (Velocity-Verlet)

Ampliamente utilizado en dinámica molecular:

$x_{n+1} = x_n + v_n h + \frac{1}{2} a_n h^2$ $v_{n+1} = v_n + \frac{1}{2} (a_n + a_{n+1}) h$

Es simétrico respecto a inversión temporal (time-reversible) y tiene un error local de posición de orden $\mathcal{O}(h^4)$, con error global $\mathcal{O}(h^2)$ para trayectorias conservativas.

4. Métodos de Runge-Kutta de Alto Orden

Los métodos de Runge-Kutta logran la precisión de una serie de Taylor de orden superior evaluando la función $\vec{f}(t, \vec{y})$ en puntos intermedios dentro del intervalo $[t_n, t_{n+1}]$, sin requerir el cálculo analítico de derivadas de orden superior.

Método de Runge-Kutta de Segundo Orden (RK2 / Método del Punto Medio)

$k_1 = \vec{f}(t_n, \vec{y}_n)$ $k_2 = \vec{f}\left(t_n + \frac{h}{2}, \vec{y}_n + \frac{h}{2} k_1\right)$ $\vec{y}_{n+1} = \vec{y}_n + h k_2$ Posee error local $\mathcal{O}(h^3)$ y error global $\mathcal{O}(h^2)$.

Método Clásico de Runge-Kutta de Cuarto Orden (RK4)

Es el estándar fundamental para ecuaciones diferenciales no rígidas en física:

$k_1 = \vec{f}(t_n, \vec{y}_n)$ $k_2 = \vec{f}\left(t_n + \frac{h}{2}, \vec{y}_n + \frac{h}{2} k_1\right)$ $k_3 = \vec{f}\left(t_n + \frac{h}{2}, \vec{y}_n + \frac{h}{2} k_2\right)$ $k_4 = \vec{f}\left(t_n + h, \vec{y}_n + h k_3\right)$ $\vec{y}_{n+1} = \vec{y}_n + \frac{h}{6}\left(k_1 + 2k_2 + 2k_3 + k_4\right)$

• Error local de truncamiento: $\mathcal{O}(h^5)$.
• Error global de truncamiento: $\mathcal{O}(h^4)$. Si se reduce el paso a la mitad ($h \to h/2$), el error global disminuye en un factor de $2^4 = 16$.

5. Control Adaptativo del Paso de Integración (RK45)

En sistemas con dinámicas multiescala (por ejemplo, órbitas elípticas excéntricas con aceleraciones violentas en el periastro y lentitud en el apoastro), un paso fijo $h$ resulta ineficiente o inestable.

Los métodos integrados de Runge-Kutta-Fehlberg (RKF45) o Dormand-Prince (DOPRI5) calculan simultáneamente dos estimaciones: una de orden 4 ($\vec{y}_{n+1}$) y una de orden 5 ($\hat{\vec{y}}_{n+1}$), reutilizando las mismas evaluaciones intermedias $k_i$. La estimación del error local es:

$\epsilon = \|\hat{\vec{y}}_{n+1} - \vec{y}_{n+1}\|$

Si $\epsilon > \text{TOL}$, el paso se rechaza y se repite con un paso reducido. Si $\epsilon \le \text{TOL}$, el paso se acepta y el siguiente paso se optimiza según:

$h_{\text{nuevo}} = S \cdot h \left(\frac{\text{TOL}}{\epsilon}\right)^{1/5}$

donde $S \approx 0.84 - 0.9$ es un factor de seguridad numérico.

6. Rigidez (Stiffness) y Estabilidad Numérica

Una EDO o sistema se califica como rígido (stiff) si existen escalas de tiempo internas con tasas de decaimiento enormemente dispares (autovalores $\lambda_i$ del Jacobiano $J = \partial \vec{f}/\partial \vec{y}$ con partes reales negativas y ratios $\max |\text{Re}(\lambda_i)| / \min |\text{Re}(\lambda_i)| \gg 1$).

Para métodos explícitos como RK4, la estabilidad exige que el paso $h$ respete la restricción impuesta por el autovalor más rápido:

$h < \frac{C}{|\lambda_{\max}|}$

incluso si la dinámica de interés a largo plazo está dominada por los modos lentos. En estos casos, se hace indispensable recurrir a métodos implícitos (como Radau o BDF) o métodos simplécticos para hamiltonianos.

Procedimiento de Implementación: De EDOs Físicas a Integración Numérica RK4

Guía Metodológica Paso a Paso

Para resolver computacionalmente la evolución temporal de cualquier sistema físico modelado por EDOs, se debe seguir un protocolo estructurado que asegure precisión, estabilidad y reproducibilidad.

Paso 1: Vectorización y Reducción a Primer Orden

1. Identificar las variables dinámicas de orden superior.
2. Definir el vector de estado $\vec{y}(t) = [y_1, y_2, \dots, y_N]^T$.
3. Expresar todas las derivadas primeras como funciones de $\vec{y}$ y $t$: $\frac{d\vec{y}}{dt} = \vec{f}(t, \vec{y})$

Paso 2: Elección del Algoritmo Adecuado

• Si el problema es no rígido y requiere alta precisión en trayectorias: emplear RK4 con paso fijo o adaptativo.
• Si el problema es hamiltoniano conservativo a largo plazo (mecánica celeste, dinámica molecular): emplear integradores simplécticos (Euler-Cromer o Verlet de velocidades).
• Si el sistema involucra tasas de decaimiento sumamente dispares (reacciones químicas, circuitos con constantes de tiempo extremas): emplear esquemas implícitos (BDF, Radau).

Paso 3: Implementación del Esquema RK4

Para cada paso de tiempo de longitud $h$:

1. Evaluar la pendiente inicial en el punto actual: $k_1 = \vec{f}(t_n, \vec{y}_n)$
2. Desplazar el estado un semipaso temporal y espacial usando $k_1$, y evaluar la segunda pendiente: $\vec{y}^{(1)} = \vec{y}_n + \frac{h}{2} k_1, \quad k_2 = \vec{f}\left(t_n + \frac{h}{2}, \vec{y}^{(1)}\right)$
3. Desplazar el estado un semipaso usando $k_2$, y evaluar la tercera pendiente: $\vec{y}^{(2)} = \vec{y}_n + \frac{h}{2} k_2, \quad k_3 = \vec{f}\left(t_n + \frac{h}{2}, \vec{y}^{(2)}\right)$
4. Proyectar el estado al final del intervalo usando $k_3$, y evaluar la cuarta pendiente: $\vec{y}^{(3)} = \vec{y}_n + h k_3, \quad k_4 = \vec{f}(t_n + h, \vec{y}^{(3)})$
5. Calcular el promedio ponderado ponderando los puntos medios: $\vec{y}_{n+1} = \vec{y}_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$
6. Actualizar el tiempo $t_{n+1} = t_n + h$.

Paso 4: Monitoreo de Invariantes Físicos y Diagnóstico de Error

En cada iteración temporal se debe verificar la conservación de magnitudes físicas invariantes del sistema:

• Energía Mecánica: $E(t) = T(t) + V(t)$. Monitorear la deriva relativa: $\Delta E(t) = \frac{|E(t) - E(0)|}{|E(0)|}$
• Momento Angular o Impulso Lineal: Confirmar que no existan torques o fuerzas espurias introducidas por el truncamiento.
• Convergencia del Paso: Si se reduce $h$ a $h/2$, el error residual debe dividirse aproximadamente por 16 en RK4.

Ejemplos Resueltos: Osciladores y Métodos Numéricos

Ejemplo 1: Oscilador Armónico Simple - Comparación Analítica y Numérica

Considérese una masa $m = 1\text{ kg}$ unida a un resorte de constante $k = 1\text{ N/m}$, de modo que la frecuencia natural es $\omega = \sqrt{k/m} = 1\text{ rad/s}$. El sistema de ecuaciones es: $\frac{dx}{dt} = v, \quad \frac{dv}{dt} = -x$ Condiciones iniciales: $x(0) = 1.0\text{ m}$, $v(0) = 0.0\text{ m/s}$. La solución analítica exacta es: $x(t) = \cos(t), \quad v(t) = -\sin(t), \quad E(t) = \frac{1}{2}(v^2 + x^2) = 0.5\text{ J}$

Efectuemos un paso de integración numérico con paso $h = 0.1\text{ s}$ comparando Euler Explícito vs RK4:

A. Método de Euler Explícito:

$x_1 = x_0 + h v_0 = 1.0 + (0.1)(0.0) = 1.000000\text{ m}$ $v_1 = v_0 + h (-x_0) = 0.0 + (0.1)(-1.0) = -0.100000\text{ m/s}$ Energía en $t = 0.1\text{ s}$: $E_{\text{Euler}} = \frac{1}{2}[(-0.1)^2 + (1.0)^2] = \frac{1}{2}(0.01 + 1.0) = 0.505000\text{ J}$ ¡La energía se ha incrementado artificialmente en un 1.0% en un solo paso!

B. Método de Runge-Kutta 4 (RK4):

Calculamos las etapas paso a paso:

1. $k_1^{(x)} = v_0 = 0.0$, $\quad k_1^{(v)} = -x_0 = -1.0$

2. Puntos intermedios con $h/2 = 0.05$: $x^{(1)} = x_0 + 0.05 k_1^{(x)} = 1.0 + 0 = 1.0$ $v^{(1)} = v_0 + 0.05 k_1^{(v)} = 0 - 0.05 = -0.05$ $k_2^{(x)} = v^{(1)} = -0.05, \quad k_2^{(v)} = -x^{(1)} = -1.0$

3. Segunda evaluación intermedia: $x^{(2)} = x_0 + 0.05 k_2^{(x)} = 1.0 + 0.05(-0.05) = 0.9975$ $v^{(2)} = v_0 + 0.05 k_2^{(v)} = 0.0 + 0.05(-1.0) = -0.05$ $k_3^{(x)} = v^{(2)} = -0.05, \quad k_3^{(v)} = -x^{(2)} = -0.9975$

4. Evaluación en el extremo $h = 0.1$: $x^{(3)} = x_0 + 0.1 k_3^{(x)} = 1.0 + 0.1(-0.05) = 0.9950$ $v^{(3)} = v_0 + 0.1 k_3^{(v)} = 0.0 + 0.1(-0.9975) = -0.09975$ $k_4^{(x)} = v^{(3)} = -0.09975, \quad k_4^{(v)} = -x^{(3)} = -0.9950$

5. Combinación ponderada final: $x_1 = 1.0 + \frac{0.1}{6}[0 + 2(-0.05) + 2(-0.05) + (-0.09975)] = 1.0 - \frac{0.1}{6}(0.29975) = 0.99500417\text{ m}$ $v_1 = 0.0 + \frac{0.1}{6}[-1.0 + 2(-1.0) + 2(-0.9975) + (-0.9950)] = -\frac{0.1}{6}(5.990) = -0.09983333\text{ m/s}$

Solución exacta en $t = 0.1\text{ s}$: $x(0.1) = \cos(0.1) = 0.995004165\text{ m}, \quad v(0.1) = -\sin(0.1) = -0.099833417\text{ m/s}$ Error absoluto en $x$ con RK4: $|0.99500417 - 0.995004165| \approx 5 \times 10^{-9}\text{ m}$. La precisión es de 8 cifras significativas.

Ejemplo 2: Péndulo Rígido No Lineal de Gran Amplitud

Consideremos un péndulo de longitud $L = 1\text{ m}$, gravedad $g = 9.8\text{ m/s}^2$: $\frac{d\theta}{dt} = \omega, \quad \frac{d\omega}{dt} = -\frac{g}{L} \sin\theta = -9.8\sin\theta$ Condiciones iniciales de gran amplitud: $\theta(0) = \frac{\pi}{2}\text{ rad}$ ($90^\circ$), $\omega(0) = 0\text{ rad/s}$. Paso de integración: $h = 0.05\text{ s}$.

Cálculo de etapas RK4:

• $k_1^{(\theta)} = \omega_0 = 0$

• $k_1^{(\omega)} = -9.8 \sin(\pi/2) = -9.8000$

• $\theta^{(1)} = \pi/2 + 0.025(0) = 1.570796$

• $\omega^{(1)} = 0 + 0.025(-9.8) = -0.245000$

• $k_2^{(\theta)} = -0.245000$

• $k_2^{(\omega)} = -9.8 \sin(1.570796) = -9.8000$

• $\theta^{(2)} = \pi/2 + 0.025(-0.245000) = 1.564671$

• $\omega^{(2)} = 0 + 0.025(-9.8000) = -0.245000$

• $k_3^{(\theta)} = -0.245000$

• $k_3^{(\omega)} = -9.8 \sin(1.564671) = -9.799015$

• $\theta^{(3)} = \pi/2 + 0.05(-0.245000) = 1.558546$

• $\omega^{(3)} = 0 + 0.05(-9.799015) = -0.489951$

• $k_4^{(\theta)} = -0.489951$

• $k_4^{(\omega)} = -9.8 \sin(1.558546) = -9.796062$

Actualización final: $\theta_1 = \frac{\pi}{2} + \frac{0.05}{6}[0 + 2(-0.245) + 2(-0.245) + (-0.489951)] = 1.570796 - 0.0122496 = 1.558547\text{ rad}$ $\omega_1 = 0 + \frac{0.05}{6}[-9.8 + 2(-9.8) + 2(-9.799015) + (-9.796062)] = -0.489951\text{ rad/s}$ La aceleración angular captura exactamente la dependencia no lineal del seno sin requerir aproximación de ángulos pequeños.

Errores Comunes en Integración Temporal y Simulación Numérica

1. Emplear Euler Explícito para Integrar Sistemas Conservativos a Largo Plazo

• Error: Utilizar el método de Euler hacia adelante para simular órbitas planetarias o moléculas durante miles de períodos.
• Consecuencia: Como el determinante del Jacobiano del método excede 1 ($\det(\mathbf{M}) = 1 + \omega^2 h^2$), la energía total del sistema diverge exponencialmente y las órbitas sufren una espiral de eyección no física.
• Corrección: Usar integradores simplécticos (como Euler-Cromer o Verlet de velocidades) para dinámica hamiltoniana, o RK4 con paso fino para intervalos cortos.

2. Confundir el Error Local de Truncamiento con el Error Global

• Error: Creer que un método con error local $\mathcal{O}(h^p)$ tendrá un error de $\mathcal{O}(h^p)$ al final de toda la simulación.
• Consecuencia: En un intervalo total $T$, el número de pasos es $N = T/h$. La acumulación de $N$ errores locales de orden $\mathcal{O}(h^{p+1})$ produce un error global de orden $\mathcal{O}(h^p)$. En RK4, el error local es $\mathcal{O}(h^5)$, pero el error acumulado global es $\mathcal{O}(h^4)$.
• Corrección: Diseñar estudios de convergencia midiendo la reducción del error acumulado global tras duplicar la resolución temporal.

3. Ignorar la Rigidez (Stiffness) del Sistema

• Error: Intentar integrar un sistema con reacciones químicas o disipación muy fuerte empleando RK4 con un paso $h$ adaptado únicamente a la escala macroscópica visible.
• Consecuencia: El autovalor rápido $\lambda_{\max} \ll 0$ fuerza $h < 2.8/|\lambda_{\max}|$ para mantener estabilidad. Si se supera este límite, la simulación diverge con valores numéricos infinitos o NaN de forma explosiva.
• Corrección: Identificar el espectro de autovalores del Jacobiano y adoptar esquemas implícitos (como BDF o Radau) para sistemas stiff.

4. No Actualizar el Tiempo $t$ en las Etapas Intermedias No Autónomas

• Error: En sistemas dependientes explícitamente del tiempo $\vec{f}(t, \vec{y})$, evaluar $k_2$ o $k_3$ manteniendo $t_n$ en lugar de $t_n + h/2$, o $k_4$ en $t_n$ en lugar de $t_n + h$.
• Consecuencia: El orden de convergencia de RK4 colapsa de cuarto orden a segundo o primer orden.
• Corrección: Asegurar rigurosamente que $k_2 = \vec{f}(t_n + h/2, \vec{y} + \frac{h}{2} k_1)$ y $k_4 = \vec{f}(t_n + h, \vec{y} + h k_3)$.

5. Falta de Verificación de Invariantes y Simetrías Físicas

• Error: Dar por válida una simulación dinámica simplemente porque el código no arroja errores ni genera NaNs, sin medir cantidades conservadas.
• Consecuencia: Derivas silenciosas en la energía mecánica, momento lineal, momento angular o volumen de fase pasan desapercibidas, falseando conclusiones científicas.
• Corrección: Programar aserciones y gráficos continuos de $\Delta E(t)/E(0)$ y de los invariantes de Noether del sistema.

Aplicaciones Físicas de Métodos de Integración Numérica

1. Mecánica Celeste y Problema de los Tres Cuerpos Gravitacional

En astrodinámica y mecánica planetaria, las trayectorias de $N$ cuerpos celestes bajo gravitación newtoniana mutua no admiten solución analítica general: $m_i \frac{d^2 \vec{r}_i}{dt^2} = \sum_{j \ne i} \frac{G m_i m_j}{|\vec{r}_j - \vec{r}_i|^3} (\vec{r}_j - \vec{r}_i)$ La integración con esquemas RK4 adaptativos o integradores simplécticos de orden superior permite calcular trayectorias de sondas interplanetarias (como las maniobras de asistencia gravitacional de las misiones Voyager o Cassini) y predecir resonancias orbitales caóticas en cinturones de asteroides.

2. Dinámica Molecular y Termodinámica Estadística

En la simulación microscópica de materiales, proteínas y fluidos densos, millones de átomos interactúan a través de potenciales interatómicos (Lennard-Jones, Coulomb): $m_i \ddot{\vec{r}}_i = -\vec{\nabla}_i \sum_{j < k} V(|\vec{r}_j - \vec{r}_k|)$ El algoritmo de Verlet de velocidades es el estándar indiscutible en este campo: su naturaleza simpléctica y reversibilidad temporal permiten conservar la energía en ensambles microcanónicos (NVE) durante millones de pasos temporales de femtosegundos sin calentamiento numérico espurio.

3. Dinámica Caótica y Convección Atmosférica (Atractor de Lorenz)

El modelo simplificado de convección térmica atmosférica de Edward Lorenz consta de tres variables acopladas no lineales: $\frac{dx}{dt} = \sigma (y - x), \quad \frac{dy}{dt} = x(\rho - z) - y, \quad \frac{dz}{dt} = x y - \beta z$ La integración precisa con RK4 revela la estructura fractal del atractor extraño y permite cuantificar los exponentes de Lyapunov, demostrando la sensibilidad extrema a las condiciones iniciales (efecto mariposa).

4. Circuitos Electrónicos No Lineales y Oscilador de Van der Pol

En ingeniería eléctrica y física de semiconductores, los osciladores auto-mantenidos (como el oscilador de Van der Pol) modelan la transición a ciclos límite: $\ddot{x} - \mu (1 - x^2) \dot{x} + x = 0$ Para parámetros $\mu \gg 1$, el sistema entra en régimen de oscilaciones de relajación altamente rígidas (stiff), donde métodos implícitos y RK adaptativo son cruciales para resolver las rápidas descargas del ciclo límite.

5. Dinámica de Haces de Partículas en Aceleradores

En física de aceleradores de alta energía (como el Gran Colisionador de Hadrones, LHC), billones de protones circulan por anillos magnéticos durante horas a velocidades ultra-relativistas. La integración de la fuerza de Lorentz relativista: $\frac{d\vec{p}}{dt} = q(\vec{E} + \vec{v} \times \vec{B})$ mediante integración simpléctica garantiza que las dimensiones transversales del haz (emitancia) no sufran dispersión artificial por acumulación de errores numéricos en miles de millones de revoluciones.

Problemas Prácticos: Simulación Dinámica e Integración Temporal

Problema 1: Decaimiento Radiactivo y Ley de Escalamiento del Error en RK4

Una muestra de radioisótopos decae según la ecuación diferencial: $\frac{dN}{dt} = -\lambda N, \quad \lambda = 0.5\text{ s}^{-1}$ con población inicial $N(0) = 1000.0\text{ átomos}$.

1. Deduzca analíticamente la expresión exacta para el factor de amplificación de un paso $N_{n+1} = G(z) N_n$ con $z = -\lambda h$ empleando el método RK4.
2. Para un paso $h = 0.2\text{ s}$, calcule $N_1$ tras el primer paso y determine el error relativo respecto a la solución analítica exacta $N(0.2) = 1000 e^{-0.1}$.

Solución Paso a Paso:

1. Las etapas RK4 son:
   • $k_1 = -\lambda N_n$
   • $k_2 = -\lambda (N_n + \frac{h}{2} k_1) = -\lambda N_n (1 - \frac{\lambda h}{2})$
   • $k_3 = -\lambda (N_n + \frac{h}{2} k_2) = -\lambda N_n [1 - \frac{\lambda h}{2} + \frac{(\lambda h)^2}{4}]$
   • $k_4 = -\lambda (N_n + h k_3) = -\lambda N_n [1 - \lambda h + \frac{(\lambda h)^2}{2} - \frac{(\lambda h)^3}{4}]$ Sustituyendo en la combinación de RK4 con $z = -\lambda h$: $G(z) = 1 + z + \frac{z^2}{2!} + \frac{z^3}{3!} + \frac{z^4}{4!}$ ¡Coincide exactamente con el polinomio de Taylor de grado 4 de la función exponencial $e^z$!
2. Con $\lambda = 0.5$ y $h = 0.2$, se tiene $z = -0.1$: $G(-0.1) = 1 - 0.1 + \frac{0.01}{2} - \frac{0.001}{6} + \frac{0.0001}{24}$ $G(-0.1) = 1 - 0.1 + 0.005 - 0.000166667 + 0.000004167 = 0.90483750$ $N_1 = 1000 \times 0.90483750 = 904.83750$ Solución exacta: $N(0.2) = 1000 e^{-0.1} = 1000 \times 0.904837418 = 904.83742$ Error absoluto: $\Delta N = |904.83750 - 904.83742| = 8 \times 10^{-5}$ Error relativo menor al $0.00001\%$.

Problema 2: Oscilador de Duffing y Dinámica No Lineal

La ecuación de un oscilador no lineal biestable sin amortiguamiento viene dada por: $\ddot{x} - x + x^3 = 0$ Transforme la ecuación en un sistema de dos variables de primer orden y evalúe la primera derivada de las variables intermedias $k_1$ y $k_2$ para la condición inicial $x_0 = 0.5\text{ m}$, $v_0 = 1.0\text{ m/s}$ con paso $h = 0.05\text{ s}$.

Solución Paso a Paso:

1. Vector de estado: $\vec{y} = \begin{pmatrix} x \\ v \end{pmatrix}$, con derivadas: $\frac{dx}{dt} = v, \quad \frac{dv}{dt} = x - x^3$
2. Etapa 1: $k_1^{(x)} = v_0 = 1.0$ $k_1^{(v)} = x_0 - x_0^3 = 0.5 - (0.5)^3 = 0.5 - 0.125 = 0.375$
3. Estados en el semipaso ($h/2 = 0.025$): $x^{(1)} = x_0 + 0.025 k_1^{(x)} = 0.5 + 0.025(1.0) = 0.525$ $v^{(1)} = v_0 + 0.025 k_1^{(v)} = 1.0 + 0.025(0.375) = 1.009375$
4. Etapa 2: $k_2^{(x)} = v^{(1)} = 1.009375$ $k_2^{(v)} = x^{(1)} - (x^{(1)})^3 = 0.525 - (0.525)^3 = 0.525 - 0.144703 = 0.380297$

Problema 3: Péndulo Invertido Estabilizado por Excitación Rápida (Péndulo de Kapitza)

El péndulo de Kapitza posee una aceleración en la base $y_b(t) = a \cos(\Omega t)$ con $\Omega \gg \omega_0$. La ecuación es: $\ddot{\theta} = -\left(\frac{g}{L} - \frac{a \Omega^2}{L}\cos(\Omega t)\right)\sin\theta$ Si $L = 0.2\text{ m}$, $g = 9.8\text{ m/s}^2$, $a = 0.01\text{ m}$ y $\Omega = 100\pi\text{ rad/s}$ ($50\text{ Hz}$):

1. Determine la frecuencia temporal más rápida presente en el sistema.
2. Establezca el límite superior del paso de integración $h$ para garantizar que el integrador RK4 capture al menos 20 pasos por ciclo de oscilación rápida.

Solución Paso a Paso:

1. El término de forzamiento tiene una frecuencia de oscilación rápida: $f_{\text{rápida}} = 50\text{ Hz} \implies T_{\text{rápida}} = \frac{1}{50} = 0.02\text{ s}$
2. Para muestrear adecuadamente la oscilación armónica de la aceleración de la base y evitar aliasing numérico o inestabilidad inducida: $h \le \frac{T_{\text{rápida}}}{20} = \frac{0.02\text{ s}}{20} = 0.001\text{ s} = 1\text{ ms}$ Si se usase un paso de $h = 0.02\text{ s}$, el integrador registraría siempre $\cos(\Omega t) = 1$, perdiendo por completo el promedio temporal que crea el pozo de potencial efectivo estabilizador.