Trabajo Práctico 1

Asignatura: Introducción al Modelado Continuo
Tema: Sistemas Dinámicos, Espacio de Fases y Problema de Kepler

---

1. Modelado en el Espacio de Estados (Primer Orden)

Partimos de la ecuación diferencial ordinaria (EDO) de segundo orden no lineal que describe el inverso del radio $u(\theta) = 1/r(\theta)$:

$$\ddot{u}(\theta) + u(\theta) - \frac{1}{\alpha} - \delta u^2(\theta) = 0$$

Para transformar esta ecuación escalar de segundo orden en un sistema de primer orden, definimos el vector de estado $\mathbf{y}(\theta) \in \mathbb{R}^2$ como:

$$\mathbf{y}(\theta) = \begin{pmatrix} y_1(\theta) \\ y_2(\theta) \end{pmatrix} \equiv \begin{pmatrix} u(\theta) \\ \dot{u}(\theta) \end{pmatrix}$$

Derivando el vector de estado respecto a la variable independiente $\theta$, obtenemos el sistema dinámico $\dot{\mathbf{y}} = \mathbf{f}(\theta, \mathbf{y})$:

$$\dot{\mathbf{y}} = \begin{pmatrix} \dot{y}_1 \\ \dot{y}_2 \end{pmatrix} = \begin{pmatrix} y_2 \\ \ddot{u} \end{pmatrix}$$

Sustituyendo la aceleración geométrica $\ddot{u}$ despejada de la ecuación original, el sistema resulta:

$$\dot{\mathbf{y}} = \begin{pmatrix} y_2 \\ -y_1 + \delta y_1^2 + \frac{1}{\alpha} \end{pmatrix}$$

---

2. Resolución del Problema de Kepler Clásico ($\delta = 0$)

Al anular el término relativista ($\delta = 0$), el sistema se torna lineal y no homogéneo. La dinámica queda gobernada por la siguiente ecuación matricial:

$$\dot{\mathbf{y}} = \mathbf{A}\mathbf{y} + \mathbf{b}$$

Donde la matriz de coeficientes $\mathbf{A}$ (Jacobiana del sistema) y el vector de carga $\mathbf{b}$ son:

$$\mathbf{A} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}, \quad \mathbf{b} = \begin{pmatrix} 0 \\ \frac{1}{\alpha} \end{pmatrix}$$

2.1. Solución del Sistema Homogéneo ($\dot{\mathbf{y}}_h = \mathbf{A}\mathbf{y}_h$)

Buscamos los autovalores de $\mathbf{A}$ resolviendo el polinomio característico $\det(\mathbf{A} - \lambda \mathbf{I}) = 0$:

$$\det \begin{pmatrix} -\lambda & 1 \\ -1 & -\lambda \end{pmatrix} = \lambda^2 + 1 = 0 \implies \lambda_{1,2} = \pm i$$

Para el autovalor $\lambda_1 = i$, hallamos su autovector asociado $\mathbf{v}_1$ resolviendo $(\mathbf{A} - i\mathbf{I})\mathbf{v}_1 = \mathbf{0}$:

$$\begin{pmatrix} -i & 1 \\ -1 & -i \end{pmatrix} \begin{pmatrix} v_{11} \\ v_{12} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} \implies v_{12} = i v_{11} \implies \mathbf{v}_1 = \begin{pmatrix} 1 \\ i \end{pmatrix}$$

La solución compleja es $\mathbf{z}(\theta) = \mathbf{v}_1 e^{i\theta}$. Utilizando la identidad de Euler ($e^{i\theta} = \cos\theta + i\sin\theta$):

$$\mathbf{z}(\theta) = \begin{pmatrix} 1 \\ i \end{pmatrix} (\cos\theta + i\sin\theta) = \begin{pmatrix} \cos\theta \\ -\sin\theta \end{pmatrix} + i \begin{pmatrix} \sin\theta \\ \cos\theta \end{pmatrix}$$

Las partes real e imaginaria constituyen la Matriz Fundamental de Soluciones $\mathbf{\Psi}(\theta)$:

$$\mathbf{\Psi}(\theta) = \begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix}$$

La solución general de la homogénea es $\mathbf{y}_h(\theta) = \mathbf{\Psi}(\theta) \mathbf{c}$, con $\mathbf{c} = \begin{pmatrix} c_1 \\ c_2 \end{pmatrix}$.

2.2. Solución Particular y Equilibrio ($\mathbf{y}_p$)

Dado que $\mathbf{b}$ es constante, buscamos un punto de equilibrio $\mathbf{y}_p$ tal que $\dot{\mathbf{y}} = \mathbf{0}$:

$$\mathbf{A}\mathbf{y}_p + \mathbf{b} = \mathbf{0} \implies \mathbf{y}_p = -\mathbf{A}^{-1}\mathbf{b}$$

Calculando la inversa $\mathbf{A}^{-1} = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}$:

$$\mathbf{y}_p = - \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} 0 \\ 1/\alpha \end{pmatrix} = \begin{pmatrix} 1/\alpha \\ 0 \end{pmatrix}$$

2.3. Solución General del Sistema

$$\mathbf{y}(\theta) = \mathbf{\Psi}(\theta) \mathbf{c} + \mathbf{y}_p = \begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} c_1 \\ c_2 \end{pmatrix} + \begin{pmatrix} 1/\alpha \\ 0 \end{pmatrix}$$

---

3. Condiciones Iniciales y Determinación de la Órbita

El problema impone condiciones sobre el radio físico $r(\theta)$ y su derivada en el perihelio ($\theta = 0$):

• $r(0) = \frac{\alpha}{1 + \epsilon} \implies y_1(0) = \frac{1 + \epsilon}{\alpha}$
• $\dot{r}(0) = 0 \implies y_2(0) = 0$ (ya que $y_2 = -r^{-2}\dot{r}$)

Planteamos el sistema algebraico para $\mathbf{c}$ evaluando en $\theta = 0$ (donde $\mathbf{\Psi}(0) = \mathbf{I}$):

$$\mathbf{y}(0) = \mathbf{I} \mathbf{c} + \mathbf{y}_p = \begin{pmatrix} c_1 \\ c_2 \end{pmatrix} + \begin{pmatrix} 1/\alpha \\ 0 \end{pmatrix} = \begin{pmatrix} \frac{1 + \epsilon}{\alpha} \\ 0 \end{pmatrix}$$

De aquí se desprende inmediatamente:

1. $c_1 + \frac{1}{\alpha} = \frac{1}{\alpha} + \frac{\epsilon}{\alpha} \implies c_1 = \frac{\epsilon}{\alpha}$
2. $c_2 + 0 = 0 \implies c_2 = 0$

Vector de solución final (para $\alpha = 1$):

$$\mathbf{y}(\theta) = \begin{pmatrix} \epsilon \cos\theta + 1 \\ -\epsilon \sin\theta \end{pmatrix}$$

---

4. Análisis Geométrico y Periodicidad

La variable de interés físico es $r(\theta) = \frac{1}{y_1(\theta)}$:

$$r(\theta) = \frac{1}{1 + \epsilon \cos\theta}$$

Conclusiones del Modelo:

1. Periodicidad: La solución es una función de $\cos\theta$, la cual posee un período de $2\pi$. Por lo tanto, $\mathbf{y}(\theta) = \mathbf{y}(\theta + 2\pi k)$ para cualquier vuelta $k$. Esto implica que, en el límite clásico ($\delta=0$), la órbita es estrictamente cerrada.
2. Espacio de Fases: En el plano $(y_1, y_2)$, la trayectoria describe una elipse (o circunferencia si $\epsilon=0$) centrada en el punto de equilibrio $(1, 0)$.
3. Clasificación de la Órbita:
   • Si $0 \le \epsilon < 1$: El denominador es siempre positivo. La trayectoria es una elipse cerrada.
   • Si $\epsilon \ge 1$: Existen ángulos donde $r(\theta) \to \infty$. La trayectoria es abierta (parábola o hipérbola).

Este desarrollo confirma que bajo las leyes de Newton, la posición del perihelio es fija, lo cual se observa gráficamente al superponerse las trayectorias para $\theta > 2\pi$.

---

Parte 3: Corrección Relativista y Precesión del Perihelio ($\delta = 0.05$)

Objetivo: Analizar el efecto del término no lineal $\delta u^2$ en la dinámica del sistema, encontrar los nuevos puntos de equilibrio, evaluar su estabilidad local y demostrar teóricamente el fenómeno de precesión del perihelio antes de proceder a la integración numérica.

El Sistema No Lineal: Fijando $\alpha = 1$ y $\delta = 0.05$, nuestro sistema vectorial autónomo es:

$$\begin{cases} \dot{y}_1 = y_2 \\ \dot{y}_2 = -y_1 + 1 + 0.05 y_1^2 \end{cases}$$

3.1. Análisis Cualitativo: Puntos de Equilibrio

Buscamos los puntos de equilibrio $\mathbf{y}_{eq} = (y_{1eq}, y_{2eq})$ donde el campo vectorial se anula ($\dot{\mathbf{y}} = \mathbf{0}$):

1. De la primera ecuación: $\dot{y}_1 = 0 \implies y_{2eq} = 0$. (La velocidad radial es nula).
2. De la segunda ecuación: $\dot{y}_2 = 0 \implies 0.05 y_{1eq}^2 - y_{1eq} + 1 = 0$.

Resolvemos esta ecuación cuadrática para la componente radial $y_{1eq}$:

$$y_{1eq} = \frac{1 \pm \sqrt{1 - 4(0.05)(1)}}{2(0.05)} = \frac{1 \pm \sqrt{0.8}}{0.1}$$

Obtenemos dos puntos de equilibrio físico:

• $\mathbf{y}_{eq1} \approx (1.056, 0)$
• $\mathbf{y}_{eq2} \approx (18.944, 0)$

Interpretación Física:

• El punto $\mathbf{y}_{eq1} \approx (1.056, 0)$ es la órbita casi circular perturbada. Notemos que si $\delta \to 0$, este punto tiende a $(1,0)$, recuperando la órbita circular de Kepler pura.
• El punto $\mathbf{y}_{eq2} \approx (18.944, 0)$ es un nuevo equilibrio puramente relativista. Corresponde a un radio muy pequeño ($r \approx 0.052$). Representa una órbita inestable muy cercana a la masa central (relacionado con el límite de captura de un agujero negro).

3.2. Linealización y Demostración Teórica de la Precesión

Para entender el movimiento alrededor de la órbita estable ($\mathbf{y}_{eq1}$), calculamos la matriz Jacobiana del sistema:

$$\mathbf{J}(y_1, y_2) = \begin{pmatrix} \frac{\partial \dot{y}_1}{\partial y_1} & \frac{\partial \dot{y}_1}{\partial y_2} \\ \frac{\partial \dot{y}_2}{\partial y_1} & \frac{\partial \dot{y}_2}{\partial y_2} \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 + 2\delta y_1 & 0 \end{pmatrix}$$

Evaluamos el Jacobiano en nuestro equilibrio de interés $\mathbf{y}_{eq1} \approx (1.056, 0)$:

$$\mathbf{J}(\mathbf{y}_{eq1}) = \begin{pmatrix} 0 & 1 \\ -1 + 2(0.05)(1.056) & 0 \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 + 0.1056 & 0 \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -0.8944 & 0 \end{pmatrix}$$

Calculamos los autovalores $\lambda$ de esta matriz resolviendo $\det(\mathbf{J} - \lambda\mathbf{I}) = 0$:

$$\lambda^2 + 0.8944 = 0 \implies \lambda_{1,2} \approx \pm i \sqrt{0.8944} \approx \pm 0.9457 i$$

El Origen de la Precesión: El sistema linealizado alrededor del equilibrio presenta autovalores imaginarios puros $\lambda = \pm i \omega$, con una frecuencia angular "espacial" $\omega \approx 0.9457$. La solución aproximada alrededor del equilibrio oscila entonces como $\cos(\omega \theta)$.

En la mecánica newtoniana ($\delta=0$), teníamos $\omega = 1$, lo que daba un período orbital exacto de $\Delta\theta = 2\pi$. Con la corrección relativista, la nueva "frecuencia" es $\omega < 1$. Por lo tanto, el "período" angular necesario para que el planeta vaya de un perihelio al próximo es:

$$\Delta\theta_{perihelio} = \frac{2\pi}{\omega} = \frac{2\pi}{0.9457} \approx 1.057 \times (2\pi) > 2\pi$$

Conclusión Analítica: El planeta necesita girar un ángulo mayor a $2\pi$ para volver a su distancia mínima. Esto significa que el eje de la elipse orbital rota en la misma dirección del movimiento en cada vuelta. A este desfase $\Delta\theta - 2\pi > 0$ se lo denomina avance o precesión del perihelio.

3.3. Resolución Numérica y Comparación

Debido a la no linealidad, la órbita exacta $\mathbf{y}(\theta)$ desde $y_1(0) = 1+\epsilon$ y $y_2(0) = 0$ se obtiene integrando numéricamente el sistema.

Metodología: Se debe aplicar un método de integración explícito de orden superior (como Runge-Kutta de orden 2 o 4) con un tamaño de paso $h \le 0.01$ en el dominio de $\theta \in [0, 6\pi]$ (equivalente a unas tres "vueltas" newtonianas).

Resultados Esperados en la Comparación:

1. Caso $\delta = 0$: Al graficar la trayectoria en el plano polar transformado a cartesiano $(x=r\cos\theta, y=r\sin\theta)$, la curva dibuja una elipse singular. Todas las vueltas se superponen exactamente (órbita invariante).
2. Caso $\delta = 0.05$: La trayectoria describe una forma similar a una elipse pero su eje mayor rota (precesa) en sentido antihorario. Las vueltas no se superponen, formando un patrón en forma de "roseta". Visualmente, el punto de máxima cercanía al Sol (perihelio) se desplaza un ángulo $\approx 0.057 \times 2\pi$ radianes en cada revolución, confirmando el cálculo de estabilidad lineal realizado previamente.