

---
## 01. Integradores Simplécticos 

Eduard Larrañaga (ealarranaga@unal.edu.co)

---


### Resumen

Se presentan los integradores LeapFrog y Verlet para resolver sistemas de ODEs.

`A. Garcia. Numerical Methods for Physics. (1999). Chapter 3 `

---

## Integrador LeapFrog 

Considere el siguiente sistema de ODEs,

\begin{equation}
\frac{d^2 x}{dt^2} = a(x),
\end{equation}

que puede ser escrito como el sistema diferencial de primer orden

\begin{align}
\begin{cases}
\frac{d x}{dt} = &v \\
\frac{dv}{dt} = & a(x).
\end{cases}
\end{align}


**Note la dependencia de la función $a = a(x(t))$ el el lado derecho de las ecuaciones.**

Discretización de la derivada temporal (utilizando diferencias centradas) en la segunda ecuación da como resultado

\begin{equation}
\frac{v(t+\Delta t) - v(t-\Delta t)}{2\Delta t} + \mathcal{O}(\Delta t^2)= a(x).
\end{equation}

En la primera ecuación diferencial también se utilizará una discretización centrada de la derivada, pero en este caso se centrará en el punto $t+\Delta t$ (esto quiere decir que allí se utilizaran los puntos $t+2\Delta t$ y $t$ en la diferencia finita). Con ello se obtiene

\begin{equation}
\frac{x(t+2\Delta t) - x(t)}{2\Delta t} + \mathcal{O}(\Delta t^2)= v(t+\Delta t).
\end{equation}

Estas dos ecuaciones de diferencias se escriben en forma iterativa asi,

\begin{align}
\begin{cases}
\frac{x_{n+2} - x_n}{2\Delta t} + \mathcal{O}(\Delta t^2) &=  v_{n+1} \\
\frac{v_{n+1} - v_{n-1}}{2\Delta t} + \mathcal{O}(\Delta t^2) &= a_n.
\end{cases}
\end{align}

Re-acomodando estas ecuaciones, se obtiene el método iteratio conocido como **Leap-Frog**,

\begin{align}
\begin{cases}
v_{n+1} &=  v_{n-1} + 2 a_n \Delta t  + \mathcal{O}(\Delta t^3)\\
x_{n+2} &=  x_n + 2 v_{n+1}\Delta t   + \mathcal{O}(\Delta t^3)
\end{cases}
\end{align}

Este es un método de orden $\mathcal{O}(\Delta t^3)$ y su nombre se debe a que la solución se construye con pasos de tamaño $2\Delta t$, con la función de posición $x$ avanzando en valores pares, $x_0, x_2, x_4, ...$, mientras que la función velocidad $v$ avanza en valores impares $v_1, v_3, v_5, ...$. 

En algunas referencias, las ecuaciones del método LeapFrog se escriben en la siguiente forma (equivalente al método de *Euler-Cromer*):

\begin{align}
\begin{cases}
v_{n+\frac{1}{2}} &=  v_{n-\frac{1}{2}} +  a_n \Delta t  + \mathcal{O}(\Delta t^3)\\
x_{n+1} &=  x_n +  v_{n+\frac{1}{2}}\Delta t   + \mathcal{O}(\Delta t^3)
\end{cases}
\end{align}



El origen del nombre de este método está en la dependencia de la función $a=a(x)$, la cual implica que su valor debe ser evaluado en un instante de tiempo centrado entre los valores antiguo y nuevo de la función $v$.

<center>
</a>
</center>

El método LeapFrog es computacionalmente mas eficiente que el método Predictor-Corrector  y requiere menos almacenamiento de información. Al ser uno de los algoritmo **simplecticos** (i.e. que preservan el área en el espacio de fase §), la conservación de la energía es respetada, aún con pasos de tiempo grandes.De esta forma, este tipo de algoritmos puede lograr una disminución apreciable en el tiempo de computo. Sin embargo, en algunos casos en los que se requieren valores de posición y velocidad mas precisos, se suelen preferir algoritmos como el Predictor-Corrector o los métodos RK. 


§ Se puede demostrar que el algoritmo Leapfrog preserva las areas encerradas por las trayectorias en el espacio de fase (p vs. x) a lo largo de la evolución del sistema. Esto implica que la energía es conservada durante la evolución.

### Iniciación del Método LeapFrog

El método Leap-Frog no se *auto-inicia* porque, además de los valores iniciales ,$(x_0, v_0)$, es necesario conocer el valore $v_{-1}$. Una de las posibilidades para inciar el algoritmo es realizar un paso de Euler hacia atrás para obtener la información requerida,
\begin{align}
v_{-1} = &v_0 -  a_0\Delta t,
\end{align}

o medio paso hacia atras si se escriben las ecuaciones en la segunda forma presentada,
\begin{align}
v_{-\frac{1}{2}} = &v_0 - \frac{1}{2} a_0 \Delta t.
\end{align}

Se puede mostrar que el método LeapFrog es mucho más rápido y más estable que el método de Euler.

---
## Método Velocidad-Verlet

Una modificación del método LeapFrog consiste en agregar una ecuación extra para poder calcular la velocidad en los mismos instantes que la posición. Este se conoce como método **velocidad-Verlet** y las ecuaciones que lo definen se obtienen a partir de las de LeapFrog, modificando la ecuación de velocidades y agregando una ecuación que incrementa en medio paso la velocidad utilizando Euler,

\begin{equation}
\begin{cases}
v_{n+\frac{1}{2}} = & v_n + \frac{1}{2} a_n \Delta t\\
x_{n+1} = & x_n + v_{n+\frac{1}{2}} \Delta t \\
v_{n+1} = & v_{n+\frac{1}{2}} + \frac{1}{2} a_{n+1} \Delta t
\end{cases}
\end{equation}

---
## Método Posición-Verlet

De forma similar,tmabién se puede modificar el método LeapFrog para obtener ls posiciones en los mismos instantes que la velocidad. Este se conoce como método **posición-Verlet** y las ecuaciones que lo definen se obtienen a partir de las de LeapFrog, modificando la ecuación de velocidades y agregando una ecuación que incrementa en medio paso la velocidad utilizando Euler,

\begin{equation}
\begin{cases}
x_{n+\frac{1}{2}} = & x_n + \frac{1}{2}v_{n} \Delta t\\
v_{n+1} = & v_n + a_{n+\frac{1}{2}} \Delta t\\
x_{n+1} = & x_{n+\frac{1}{2}} + \frac{1}{2}v_{n+1} \Delta t \\
\end{cases}
\end{equation}

---
## Integrador Verlet 

Considere de nuevo el sistema de ODEs


\begin{equation}
\frac{d^2 x}{dt^2} = a(x)
\end{equation}

o reescrito, introduciendo el concepto de velocidad, en la forma

\begin{align}
\begin{cases}
\frac{d x}{dt} &= v \\
\frac{d^2 x}{dt^2} &=  a(x).
\end{cases}
\end{align}

Utilizaremos diferencias centradas para la primera y segunda derivada y con ello obtenemos las formulas iterativas

\begin{align}
\begin{cases}
\frac{x_{n+1} - x_{n-1}}{2\Delta t}  + \mathcal{O} (\Delta t^2) &= v_n \\
\frac{x_{n+1} - 2x_n + x_{n-1}}{\Delta t^2} + \mathcal{O} (\Delta t^2)  &= a_n.
\end{cases}
\end{align}

Re-acomodando estas ecuaciones da como resultado

\begin{align}
\begin{cases}
v_n = &\frac{x_{n+1} - x_{n-1}}{2\Delta t}  + \mathcal{O} (\Delta t^2) \\
x_{n+1} = &2x_n - x_{n-1} +  a_n\Delta t^2 + \mathcal{O} (\Delta t^4) .  
\end{cases}
\end{align}

Este esquema se conoce como el método de **Verlet**. Nótese que, de nuevo, el método no es auto-iniciable debido a que se necesita los valores de $x_n$ y $x_{n-1}$ para obtener $x_{n+1}$. Además, es importante resaltar que este método puede obtener la evolución de la posición de forma independiente de los valores de la velocidad.

---

## Ventajas de los Métodos Simplécticos

Los métodos simplécticos presentados poseen las siguientes ventajas sobre otros algoritmos:

### 1. Invarianza bajo inversión temporal

Al igual que la mecánica Newtoniana, los algoritmos simplécticos son invariantes bajo inversión temporal. Para comprender esta característica, considere que se incia el algoritmo en un tiempo $t_0$ con una condición $(x_0,v_0) = (X,V)$ y se avanza un paso para encontrar los valores $(x_1,v_1)$ mediante alguno de los métodos (e.g. velocidad-Verlet). Al reemplazar en las ecuaciones se tiene  

\begin{equation}
\begin{cases}
v_{\frac{1}{2}} = & V + \frac{1}{2} a_0 \Delta t\\
x_{1} = & X + v_{\frac{1}{2}} \Delta t \\
v_{1} = & v_{\frac{1}{2}} + \frac{1}{2} a_{1} \Delta t
\end{cases}
\end{equation}

y por ello

\begin{equation}
\begin{cases}
x_{1} = & x_0 + v_0 \Delta t + \frac{1}{2} a_0 \Delta t^2 \\
v_{1} = & v_0 + \frac{1}{2} (a_0 + a_{1}) \Delta t.
\end{cases}
\end{equation}

Si ahora se inicia el algoritmo con la condición inicial $(x_0^r, v_0^r) = (x_1, -v_1)$ obtenidos arriba, las ecuaciones del método resultan en

\begin{equation}
\begin{cases}
v^r_{\frac{1}{2}} = & -v_1 + \frac{1}{2} a^r_0 \Delta t\\
x^r_{1} = & x_1 + v^r_{\frac{1}{2}} \Delta t \\
v^r_{1} = & v^r_{\frac{1}{2}} + \frac{1}{2} a^r_{1} \Delta t
\end{cases}
\end{equation}

con lo que se obtiene

\begin{equation}
\begin{cases}
x^r_{1} = & x_1 -v_1 \Delta t + \frac{1}{2} a^r_0 \Delta t^2 \\
v^r_{1} = & -v_1 + \frac{1}{2} a^r_0 \Delta t + \frac{1}{2} a^r_{1} \Delta t
\end{cases}
\end{equation}

Al reemplazar los valores obtenidos en la primera parte,

\begin{equation}
\begin{cases}
x^r_{1} = & x_0 + v_0 \Delta t + \frac{1}{2} a_0 \Delta t^2 - \left( v_0 + \frac{1}{2} (a_0 + a_{1}) \Delta t \right) \Delta t + \frac{1}{2} a^r_0 \Delta t^2 \\
v^r_{1} = & -v_0 - \frac{1}{2} (a_0 + a_{1}) \Delta t   + \frac{1}{2} a^r_0 \Delta t + \frac{1}{2} a^r_{1} \Delta t
\end{cases}
\end{equation}

\begin{equation}
\begin{cases}
x^r_{1} = & x_0   - \frac{1}{2}  a_{1} \Delta t^2 + \frac{1}{2} a^r_0 \Delta t^2 \\
v^r_{1} = & -v_0 - \frac{1}{2} (a_0 + a_{1} - a^r_0 - a^r_{1}) \Delta t  . 
\end{cases}
\end{equation}

Finalmente, identificando que $a^r_0 \equiv a_1$ y $a^r_1 \equiv a_0$ se tiene el resultado buscado

\begin{equation}
\begin{cases}
x^r_{1} = & x_0    \\
v^r_{1} = & -v_0. 
\end{cases}
\end{equation}






## Conservación del Área en el Espacio de Fase (Conservación de la Energía)

El adjetivo *Simpléctico* hace referencia a la conservación del área en el espacio de fase del sistema físico. 

Considere los cuatro vértices de un cuadrado en el espacio de fase en un instante de tiempo $t$: $(x,p)$, $(x+dx,p)$, $(x,p+dp)$ y $(x+dx,p+dp)$, los cuales encierran un área $dA=dxdp$. Para un tiempo posterior, $t'$, los cuatro puntos cambiarán de ubicación, pero si la evolución es conservativa (de acuerdo con el teorema de Liouville), la nueva área encerrada será igual a la inicial, $dA' = dA$.

La evolución temporal se puede representar como una transformación de coordenadas

\begin{equation}
\begin{cases}
x' =& x'(x,p)\\
p' =& p'(x,p)
\end{cases}
\end{equation}

Bajo la transformación de coordenadas, el área encerrada en el espacio de fase cambia de acuerdo con

\begin{equation}
dA' = \det{J} dA
\end{equation}

donde $\det{J}$ es el *Jacobiano* de la transformación y

\begin{equation}
J = \begin{pmatrix} 
\frac{\partial x'}{\partial x} & \frac{\partial x'}{\partial p} \\
\frac{\partial p'}{\partial x} & \frac{\partial p'}{\partial p}
\end{pmatrix}.
\end{equation}

Si la transformación (i.e. la evolución) es simpléctica, se tiene $\det{J} = 1$.

Para comprobar que los algoritmos presentados son simplécticos, calcularemos el Jacobiano. Por ejemplo, tomaremos las ecuaciones del método velocidad-Verlet. Al aplicar la primera de ellas a los puntos $(x_0, v_0)$ y $(x_0+\delta x_0, v_0 + \delta v_0)$, se obtienen las ecuaciones


\begin{equation}
v_{\frac{1}{2}} =  v_0 + \frac{1}{2} a_0 \Delta t \hspace{2cm}; \hspace{2cm}
v_{\frac{1}{2}} + \delta v_{\frac{1}{2}} =  v_0 +\delta v_0 + \frac{1}{2} a(x_0 + \delta x_0) \Delta t
\end{equation}

Restando estas relaciones se tiene

\begin{align}
\delta v_{\frac{1}{2}} = & \delta v_0 + \frac{1}{2} a(x_0 + \delta x_0) \Delta t - \frac{1}{2} a_0 \Delta t\\
\delta v_{\frac{1}{2}} = & \delta v_0 + \frac{1}{2} \left(a(x_0 + \delta x_0)  - a_0\right) \Delta t\\
\delta v_{\frac{1}{2}} = & \delta v_0 + \frac{1}{2} \frac{\left(a(x_0 + \delta x_0)  - a_0\right)}{\delta x_0} \delta x_0 \Delta t.
\end{align}

Tomando el límite cuando $\delta x_0 \rightarrow 0$, se puede escribir

\begin{align}
\delta v_{\frac{1}{2}} = & \delta v_0 + \frac{1}{2} \left. \frac{da}{dx} \right|_{x_0} \delta x_0 \Delta t.
\end{align}

De esta forma, se puede escribir la transformación en forma matricial,

\begin{equation}
\begin{pmatrix}
\delta x_0 \\
\delta v_{\frac{1}{2}}
\end{pmatrix} = B 
\begin{pmatrix}
\delta x_0 \\
\delta v_0
\end{pmatrix}
\end{equation}

donde 

\begin{equation}
B = 
\begin{pmatrix}
1 & 0 \\
\frac{1}{2}a'(x_0) \Delta t & 1
\end{pmatrix}.
\end{equation}

En la misma forma, pero utilizando la segunda y tercera ecuación del método de velocidad-Verlet, se obtienen las relaciones

\begin{equation}
\begin{pmatrix}
\delta x_1 \\
\delta v_{\frac{1}{2}}
\end{pmatrix} = C 
\begin{pmatrix}
\delta x_0 \\
\delta v_{\frac{1}{2}}
\end{pmatrix}
\hspace{2cm} ; \hspace{2cm}
\begin{pmatrix}
\delta x_1 \\
\delta v_1
\end{pmatrix} = D 
\begin{pmatrix}
\delta x_1 \\
\delta v_{\frac{1}{2}}
\end{pmatrix}
\end{equation}

donde 

\begin{equation}
C = 
\begin{pmatrix}
1 & \Delta t \\
0 & 1
\end{pmatrix}
\hspace{2cm} y \hspace{2cm}
D = 
\begin{pmatrix}
1 & 0 \\
\frac{1}{2}a'(x_1) \Delta t & 1
\end{pmatrix}
.
\end{equation}

De estas realciones se obtiene la transformación 

\begin{equation}
\begin{pmatrix}
\delta x_1 \\
\delta v_1
\end{pmatrix} = J 
\begin{pmatrix}
\delta x_0 \\
\delta v_0
\end{pmatrix}
\end{equation}

con $J = DCB$. El determinante de esta matriz (Jacobiano) es

\begin{equation}
\det J = \det D \det C \det B.
\end{equation}

Ya que cada una de las matrices construidas satisface $\det B = \det C = \det D = 1$, se tiene que $\det J = 1$ y por lo tanto el método de itnegración es *simpléctico*. El mismo proceso puede ser aplicado a los metodos LeapFrog, posición-Verlet y Verlet.


```python

```

73a0d4a08f901a12b8c083047e3e84eb35aa08a2
ipynb
Jupyter Notebook
07. ODEs II. Integradores Simplecticos/01. Integradores Simplecticos.ipynb
ashcat2005/AstrofisicaComputacional2022
67463ec4041eb08c0f326792fed0dcf9e970e9b7
MIT
2022-03-08T06:18:56.000Z
2022-03-10T04:55:53.000Z
07. ODEs II. Integradores Simplecticos/01. Integradores Simplecticos.ipynb
ashcat2005/AstrofisicaComputacional2022
67463ec4041eb08c0f326792fed0dcf9e970e9b7
MIT
07. ODEs II. Integradores Simplecticos/01. Integradores Simplecticos.ipynb
ashcat2005/AstrofisicaComputacional2022
67463ec4041eb08c0f326792fed0dcf9e970e9b7
MIT
2022-03-09T17:47:43.000Z
2022-03-21T02:29:36.000Z
Qwen/Qwen-72B
1. YES
2. YES