Integración de Verlet

El algoritmo de Verlet es un procedimiento para la integración numérica de ecuaciones diferenciales ordinarias de segundo orden con valores iniciales conocidos (problema de Cauchy). Es particularmente apropiado en las situaciones en que la expresión de la segunda derivada solo es función de las variables, dependiente o independiente, sin participar la primera derivada. Este es el caso de numerosos problemas de la dinámica newtoniana, por lo que se emplea frecuentemente en astronomía y mecánica molecular.

Dado un sistema de coordenadas, la posición de una partícula a lo largo del tiempo puede ser expresada por una función de las coordenadas del vector y del tiempo. En el instante t

Ec.(1)  : Posición = ${\displaystyle {\vec {x}}(t)}$; Velocidad = ${\displaystyle {\frac {\mathrm {d} {\vec {x}}(t)}{\mathrm {d} t}}}$; Aceleración = ${\displaystyle {\frac {\mathrm {d} ^{2}{\vec {x}}(t)}{\mathrm {d} t^{2}}}}$.

En 1967, el matemático francés Loup Verlet presentó la primera versión de su modelo,[1]​ denominada Integración de Verlet, caracterizada por su simplicidad sin pérdida de exactitud y estabilidad. Posteriormente, en 1985,[2]​ se propuso una ligera corrección a la Integración de Verlet, conocida como algoritmo de Verlet con velocidad, que mejora la precisión y estabilidad de las soluciones.

Algoritmo de Verlet. Integración

Desarrollando en serie de Taylor el valor de la función para dos incrementos de la variable independiente: dt y - dt se obtiene:

Ec.( 2 )  : ${\displaystyle {\vec {x}}(t+\Delta t)={\vec {x}}(t)+{\vec {v}}(t)\Delta t+{\frac {{\vec {a}}(t)\Delta t^{2}}{2}}+{\frac {{\vec {b}}(t)\Delta t^{3}}{6}}+O(\Delta t^{4})\,}$

Ec.( 3 )  : ${\displaystyle {\vec {x}}(t-\Delta t)={\vec {x}}(t)-{\vec {v}}(t)\Delta t+{\frac {{\vec {a}}(t)\Delta t^{2}}{2}}-{\frac {{\vec {b}}(t)\Delta t^{3}}{6}}+O(\Delta t^{4}).\,}$

Donde ${\displaystyle {\vec {x}}}$ es la posición, ${\displaystyle {\vec {v}}}$ la velocidad, ${\displaystyle {\vec {a}}}$ la aceleración y ${\displaystyle {\vec {b}}}$ la sobreaceleración .

Sumando ambas ecuaciones, limitando la serie hasta el término de cuarto grado, y reagrupando términos se obtiene la ley de recursión para las posiciones:

Ec.( 4 )  : ${\displaystyle {\vec {x}}(t+\Delta t)=2{\vec {x}}(t)-{\vec {x}}(t-\Delta t)+{\vec {a}}(t)\Delta t^{2}+O(\Delta t^{4}).\,}$

Al despreciar el último término el error cometido es de cuarto grado, inferior a la mayoría de los algoritmos convencionales.

Como puede comprobarse, la Ec.(4) proporciona la trayectoria sin intervenir la velocidad, lo que simplifica los cálculos. Caso de necesitar conocer las velocidades, estas se obtienen directamente por simple diferencia de las Ec.(2) y Ec.(3):

Ec. ( 5 )  : ${\displaystyle {\vec {v}}(t)=[{\vec {x}}(t+\Delta t)-{\vec {x}}(t-\Delta t)]/(2\Delta t)}$

Inicialización

Las Ec.(4) y Ec.(5) son característica de un método de resolución implícito (multipaso). Al comenzar la serie, se halla el tercer término (${\displaystyle {\vec {x}}(t+\Delta t)}$), en función del segundo (${\displaystyle {\vec {x}}(t)}$) y el primero ( ${\displaystyle {\vec {x}}(t-\Delta t)}$). El primer término, es un dato del problema, pero el segundo es desconocido y la Ec.(5) no permite su cálculo. Para obviar esta limitación, el segundo valor de la serie (${\displaystyle {\vec {x}}(t)}$) se calcula a partir de los valores iniciales con ayuda de algún método explícito que permita esta operación, por ejemplo el de Euler,

Ec. ( 6 )  : ${\displaystyle {\vec {x}}(t+\Delta t)={\vec {x}}(t)+\Delta t{\vec {v}}(t)}$

Y ya se tienen todos los términos necesarios para proseguir el cálculo recursivo.

Esta discontinuidad puede ser origen de cierta inestabilidad, aunque, si el valor de ${\displaystyle \Delta t}$ es suficientemente pequeño, suele corregirse automáticamente.

Aplicación: Problema de los dos cuerpos

En el denominado problema de los dos cuerpos se trata de hallar las trayectorias de dos cuerpos sometidos exclusivamente a su fuerza gravitatoria mutua. Como es conocido el centro de masas del conjunto puede ser tomado como origen de coordenadas y las órbitas se mueven en un plano. Esto permite simplificar el cálculo. Adoptando un sistema de coordenadas y posiciones iniciales como se muestra en el gráfico, el problema se reduce a un sistemas de cuatro variables dependientes, las coordenadas de posición de los dos cuerpos, una variable independiente, el tiempo, todo ello ligado por cuatro ecuaciones que se deben resolver secuencialmente en cada paso de la iteración.

La figura muestra los cálculos realizados en una Hoja de Cálculo y la salida de un programa gráfico escrito en javascript

El dibujo se ha obtenido trazando varios periodos con el propósito de mostrar la falta de coincidencia entre órbitas consecutivas, motivado por haber seleccionado un ${\displaystyle \Delta t}$ excesivamente grande. Resulta sencillo mejorar la precisión reduciendo el intervalo o, como se muestra a continuación, aplicando el siguiente algoritmo.

Algoritmo de Verlet con velocidad

Como mejora del original algoritmo, se ha propuesto una fórmula aplicable desde el principio donde interviene explicitamente la velocidad, que es calculada por un procedimiento algo más refinado.

En la expresión del desarrollo en serie,

Ec.( 7 )  : ${\displaystyle {\vec {x}}(t+\Delta t)={\vec {x}}(t)+{\vec {v}}(t)\Delta t+{\frac {{\vec {a}}(t)\cdot \Delta t^{2}}{2}}+O(\Delta t^{3})\ +...,}$

Limitando la serie al término de segundo grado,[3]​ el error en la posición es de orden tercero, se encuentra la forma iterativa para determinar la posición en función de velocidad (d${\displaystyle {\vec {x}}(t)}$ / dt) y aceleración (d^2${\displaystyle {\vec {x}}(t)}$ / dt^2), ambos datos del problema). Para calcular las velocidades, operando sobre las expresiones obtenidas en las Ec.(5) y Ec.(6)[4]​ se llega a la siguiente expresión:

Ec.(8)  : ${\displaystyle {\vec {v}}(t+\Delta t)={\vec {v}}(t)+\Delta t\cdot [d^{2}{\vec {x}}(t+\Delta t)/dt^{2}+d^{2}{\vec {x}}(t)/dt^{2}]/2}$

Las Ec.(7) y Ec.(8) presentan gran analogía con las utilizadas en el procedimiento de Integración del salto de rana; en Verlet, posición y velocidad se obtienen en el mismo punto de la trayectoria, mientras que en salto de rana se calculan intercalados.

Inicialización

En la Ec.(8), para calcular ${\displaystyle {\vec {v}}(t+\Delta t)}$ se necesita el valor previo de ${\displaystyle {\vec {a}}(t+\Delta t)}$. Pero a diferencia de la Integración de Verlet, no es necesario recurrir a ninguna suposición externa al propio procedimiento

Llamando t0, v0 y a0 a los valores de la respectivas variables en el instante inicial y t1 = t0 + ${\displaystyle \Delta t}$, v1 y a1 para el siguiente término:

A partir de la Ec.( 1 )  : ${\displaystyle {\vec {a}}_{1}(t_{0})}$

A partir de la Ec.( 7 )  : ${\displaystyle {\vec {x}}_{1}(t_{0}+\Delta t)={\vec {x}}_{0}(t_{0})+{\vec {v}}(t_{0})\cdot \Delta t+{\frac {{\vec {a}}_{0}(t_{0})\cdot \Delta t^{2}}{2}}}$

Con el valor de ${\displaystyle {\vec {x}}_{1}(t_{1})}$  : ${\displaystyle {\vec {a}}_{1}(t_{0}+\Delta t)}$

Continua la secuencia  : ${\displaystyle {\vec {x}}(t+\Delta t)={\vec {x}}(t)+{\vec {v}}(t+\Delta t)+{\frac {{\vec {a}}(t)\cdot \Delta t^{2}}{2}}}$

Aplicación: Problema de los dos cuerpos

Para una primera comparación de la precisión de ambos procedimientos en la siguiente figura se ha resuelto el mismo caso, con dos valores de ${\displaystyle \Delta t}$, 0,005 y 0,001.

Como puede comprobarse, para ${\displaystyle \Delta t=0,005}$, el trazo de la figura es más neto que con el algoritmo de Integración de Verlet, lo que implica mayor precisión. No obstante, no puede considerarse completamente satisfactorio. En cambio, haciendo ${\displaystyle \Delta t=0,001}$ el perfil de las órbitas parece nítido por completo, lo que significa que en las sucesivas revoluciones se ha utilizado idéntico camino(hasta donde alcanza la precisión observable), siendo preferible esta vía para mejorar los resultados.

Referencias

1. Loup Verlet "Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules", Physical Review 159 pp. 98-103 (1967)
2. William C. Swope, Hans C. Andersen, Peter H. Berens and Kent R. Wilson "A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters", Journal of Chemical Physics 76 pp. 637-649 (1982)
3. Obsérvese la analogía con la fórmula que expresa el espacio recorrido en el caso del movimiento uniformemente acelerado
4. Mecanica 3d: python y el algoritmo de Verlet. J. F. Rojas, R. Martínez y M. A. Morales

• Datos: Q2984997