Resumen

Este proyecto evalúa el desempeño de tres integradores Runge-Kutta clásicos —Euler explícito, Heun (RK2) y Runge-Kutta de cuarto orden (RK4)— al resolver el seguimiento lagrangiano $\dot{\vec{x}}(t) = \vec{u}(\vec{x})$ de partículas en dos flujos 2D incompresibles: (A) un flujo radial-rotacional con punto silla, lineal y con solución cerrada $\vec{x}(t)=e^{At}\vec{x}_0$; y (B) una celda de convección Rayleigh-Bénard no lineal. Cuatro experimentos aíslan trayectoria individual, orden de convergencia empírico, conservación de área (Liouville) y advección de una nube de 1500 partículas. Los resultados confirman las pendientes teóricas $1$, $2$ y $4$ con desviaciones inferiores al $7\%$, y la predicción analítica de pérdida de área para Euler ($82.81\%$ del área inicial) coincide con la medición a seis cifras significativas.

Palabras clave: integración numérica de EDO · Runge-Kutta · seguimiento lagrangiano · flujo incompresible · teorema de Liouville · orden de convergencia.

1Introducción

En mecánica de fluidos coexisten dos descripciones complementarias del movimiento. La descripción euleriana fija un sistema de referencia espacial y describe los campos $\vec{u}(\vec{x},t)$, $p(\vec{x},t)$ en cada punto del dominio. La descripción lagrangiana, en cambio, sigue a una partícula material individual:

$$ \frac{d\vec{x}}{dt} = \vec{u}\bigl(\vec{x}(t)\bigr), \qquad \vec{x}(0) = \vec{x}_0 $$(1)

La descripción lagrangiana es la herramienta natural para estudiar fenómenos donde importa la historia individual: dispersión de contaminantes, tiempos de residencia en reactores, transporte de sedimentos, modelado de aerosoles, agregación en flujos turbulentos. En todos estos casos el problema se reduce a integrar (1) para múltiples condiciones iniciales — lo que hace del estudio comparativo de integradores una pregunta directamente relevante.

2Marco teórico · fluidos

El movimiento de un fluido viscoso de densidad constante obedece Navier-Stokes incompresibles:

$$ \rho(\partial_t \vec{u} + \vec{u}\cdot\nabla\vec{u}) = -\nabla p + \mu \nabla^2 \vec{u} + \vec{f}, \qquad \nabla \cdot \vec{u} = 0 $$(2)

2.1 · Teorema de Liouville

Sea $\vec{u}$ incompresible y $\Phi_t$ el flujo asociado. Entonces para cualquier región $\Omega_0$:

$$ \text{area}(\Phi_t(\Omega_0)) = \text{area}(\Omega_0) \quad \forall t \geq 0 $$(3)

Demostración: por la fórmula de transporte de Reynolds, $\frac{d}{dt}\text{area}(\Phi_t(\Omega_0)) = \int \nabla \cdot \vec{u}\, dA = 0$. Este teorema actúa como test físico en el Experimento 3.

3Marco teórico · integración numérica

3.1 · Euler explícito (orden 1)

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

Trunca Taylor en el primer término no trivial. Error de truncamiento $O(h^2)$, orden global $p=1$.

3.2 · Heun (RK2, orden 2)

$$\begin{aligned} \vec{k}_1 &= \vec{f}(\vec{x}_n) \\ \vec{k}_2 &= \vec{f}(\vec{x}_n + h\, \vec{k}_1) \\ \vec{x}_{n+1} &= \vec{x}_n + \tfrac{h}{2}(\vec{k}_1 + \vec{k}_2) \end{aligned}$$ (5)

Predictor-corrector con promedio. Error local $O(h^3)$, orden global $p=2$.

3.3 · Runge-Kutta clásico (RK4, orden 4)

$$\begin{aligned} \vec{k}_1 &= \vec{f}(\vec{x}_n) \\ \vec{k}_2 &= \vec{f}(\vec{x}_n + \tfrac{h}{2}\vec{k}_1) \\ \vec{k}_3 &= \vec{f}(\vec{x}_n + \tfrac{h}{2}\vec{k}_2) \\ \vec{k}_4 &= \vec{f}(\vec{x}_n + h\, \vec{k}_3) \\ \vec{x}_{n+1} &= \vec{x}_n + \tfrac{h}{6}(\vec{k}_1 + 2\vec{k}_2 + 2\vec{k}_3 + \vec{k}_4) \end{aligned}$$ (6)

Los pesos $\tfrac{1}{6}, \tfrac{2}{6}, \tfrac{2}{6}, \tfrac{1}{6}$ son los de la regla de Simpson. Orden global $p=4$.

4Casos de estudio

4.1 · Caso A — punto silla hiperbólico

Campo lineal $\vec{u}_A(\vec{x}) = A\vec{x}$ con $A$ matriz simétrica de traza cero:

$$ A = \begin{pmatrix} a & c \\ c & -a \end{pmatrix}, \quad c = \sqrt{\alpha/\rho} $$(7)

Solución analítica cerrada: $\vec{x}(t) = e^{At}\vec{x}_0$. Banco de pruebas ideal. Se fija $a=0.5$, $c=1.0$.

4.2 · Caso B — celdas de Rayleigh-Bénard

Campo no lineal y periódico, sin solución cerrada. Mosaico de celdas convectivas con vorticidad alternante. La única referencia práctica es la propia solución de RK4 a paso fino.

5Experimento 1 · Trayectoria individual

Comparación punto a punto contra $\vec{x}(t) = e^{At}\vec{x}_0$. Condición inicial $\vec{x}_0 = (1.5, -0.8)$, $h=0.25$, $t_f = 2$.

Trayectoria comparada — Euler, Heun, RK4 vs analítica
Fig 1. Caso A — trayectoria desde $\vec{x}_0 = (1.5, -0.8)$, $h=0.25$, $t_f=2.0$. Euler se aleja visiblemente de la curva analítica; Heun y RK4 son indistinguibles a esta escala.
Error global vs tiempo
Fig 2. Error global $\|\vec{x}_{\text{num}}(t)-\vec{x}_{\text{exact}}(t)\|_2$ con $h=0.25$. RK4 mantiene error 3–4 órdenes de magnitud por debajo de Euler.

6Experimento 2 · Orden de convergencia

Ajuste log-log del error final $\|\vec{x}_{\text{num}}(t_f) - \vec{x}_{\text{exact}}(t_f)\|_2$ vs paso $h$. Las pendientes empíricas deben coincidir con $p = 1, 2, 4$.

Orden de convergencia empírico
Fig 3. Pendientes medidas: Euler $0.93$ (teórico $1$), Heun $1.94$ (teórico $2$), RK4 $3.93$ (teórico $4$). Desviaciones inferiores al $7\%$ confirman la teoría.

7Experimento 3 · Conservación de área (Liouville)

Se sigue un cuadrado material y se mide $A(t)/A(0)$. La teoría predice conservación exacta para flujos incompresibles.

Para Euler en el Caso A se puede derivar analíticamente la pérdida sistemática:

$$ \frac{A(t_N)}{A(0)} = (1 - h^2(a^2 + c^2))^N $$(8)

Con $a=0.5$, $c=1.0$, $h=0.1$, $N=15$ predice $A/A_0 = 0.8281$ — confirmado experimentalmente a seis cifras.

Conservación de área
Fig 4. Cociente $A(t)/A(0)$ con $h=0.1$. Heun y RK4 preservan el área (Liouville exacta); Euler exhibe deriva sistemática predicha analíticamente.
Deformación del cuadrado material
Fig 5. Deformación del cuadrado material. Los tres integradores producen estiramientos cualitativamente similares, pero solo Heun y RK4 preservan el área cuantitativamente.

8Experimento 4 · Advección de una nube

Se advectan 1500 partículas pasivas en el campo no lineal del Caso B. Sin solución cerrada, la referencia práctica es RK4 a paso fino.

Advección lagrangiana de 1500 partículas
Fig 6. Caso B — advección lagrangiana de 1500 partículas. La nube evoluciona desde un disco hasta una estructura ramificada; los tres integradores producen mapas cualitativamente similares para $t=5$.

9Conclusiones

Los cuatro experimentos verifican empíricamente las predicciones teóricas del Capítulo 25 de Chapra-Canale:

  • Las pendientes de convergencia coinciden con los enteros $1$, $2$, $4$ a menos del $7\%$.
  • La pérdida de área analítica para Euler en el Caso A se reproduce con seis cifras significativas.
  • Heun y RK4 preservan el invariante físico de Liouville hasta el error de máquina; Euler exhibe deriva sistemática.
  • En flujo no lineal (Caso B), Heun y RK4 producen trayectorias mutuamente indistinguibles para $h$ moderado, mientras Euler diverge cualitativamente para $t$ grande.

El estudio confirma que para problemas lagrangianos en flujos incompresibles, RK4 es la opción por defecto: ofrece error de máquina sobre el invariante físico con un costo computacional de cuatro evaluaciones por paso.