Integración de Verlet

Algoritmo de integración numérica

La integración de Verlet ( pronunciación francesa: [vɛʁˈlɛ] ) es un método numérico utilizado para integrar las ecuaciones de movimiento de Newton . ^[1] Se utiliza con frecuencia para calcular trayectorias de partículas en simulaciones de dinámica molecular y gráficos por computadora . El algoritmo fue utilizado por primera vez en 1791 por Jean Baptiste Delambre y ha sido redescubierto muchas veces desde entonces, la más reciente por Loup Verlet en la década de 1960 para su uso en dinámica molecular. También fue utilizado por PH Cowell y ACC Crommelin en 1909 para calcular la órbita del cometa Halley , y por Carl Størmer en 1907 para estudiar las trayectorias de partículas eléctricas en un campo magnético (de ahí que también se le llame método de Størmer ). ^[2] El integrador de Verlet proporciona una buena estabilidad numérica , así como otras propiedades que son importantes en los sistemas físicos, como la reversibilidad temporal y la preservación de la forma simpléctica en el espacio de fases , sin ningún costo computacional adicional significativo respecto del método de Euler simple .

Básico Størmer-Verlet

Para una ecuación diferencial de segundo orden del tipo con condiciones iniciales y , se puede obtener una solución numérica aproximada en los tiempos con tamaño de paso mediante el siguiente método: ${\displaystyle {\ddot {\mathbf {x} }}(t)=\mathbf {A} {\bigl (}\mathbf {x} (t){\bigr )}}$ ${\displaystyle \mathbf {x} (t_{0})=\mathbf {x} _{0}}$ ${\displaystyle {\dot {\mathbf {x} }}(t_{0})=\mathbf {v} _{0}}$ ${\displaystyle \mathbf {x} _{n}\approx \mathbf {x} (t_{n})}$ ${\displaystyle t_{n}=t_{0}+n\,\Delta t}$ ${\displaystyle \Delta t>0}$

1. colocar , ${\textstyle \mathbf {x} _{1}=\mathbf {x} _{0}+\mathbf {v} _{0}\,\Delta t+{\tfrac {1}{2}}\mathbf {A} (\mathbf {x} _{0})\,\Delta t^{2}}$
2. para n = 1, 2, ... iterar $${\displaystyle \mathbf {x} _{n+1}=2\mathbf {x} _{n}-\mathbf {x} _{n-1}+\mathbf {A} (\mathbf {x} _{n})\,\Delta t^{2}.}$$

Ecuaciones de movimiento

La ecuación de movimiento de Newton para sistemas físicos conservativos es

${\displaystyle {\boldsymbol {M}}{\ddot {\mathbf {x} }}(t)=F{\bigl (}\mathbf {x} (t){\bigr )}=-\nabla V{\bigl (}\mathbf {x} (t){\bigr )},}$

o individualmente

${\displaystyle m_{k}{\ddot {\mathbf {x} }}_{k}(t)=F_{k}{\bigl (}\mathbf {x} (t){\bigr )}=-\nabla _{\mathbf {x} _{k}}V\left(\mathbf {x} (t)\right),}$

dónde

• ${\displaystyle t}$es el momento,
• ${\displaystyle \mathbf {x} (t)={\bigl (}\mathbf {x} _{1}(t),\ldots ,\mathbf {x} _{N}(t){\bigr )}}$es el conjunto del vector de posición de los objetos, ${\displaystyle N}$
• ${\displaystyle V}$es la función potencial escalar,
• ${\displaystyle F}$es el gradiente negativo del potencial , que da el conjunto de fuerzas sobre las partículas,
• ${\displaystyle {\boldsymbol {M}}}$es la matriz de masa , típicamente diagonal con bloques con masa para cada partícula. ${\displaystyle m_{k}}$

Esta ecuación, para diversas opciones de la función potencial , puede utilizarse para describir la evolución de diversos sistemas físicos, desde el movimiento de moléculas en interacción hasta la órbita de los planetas . ${\displaystyle V}$

Después de una transformación para llevar la masa al lado derecho y olvidar la estructura de múltiples partículas, la ecuación se puede simplificar a

${\displaystyle {\ddot {\mathbf {x} }}(t)=\mathbf {A} {\bigl (}\mathbf {x} (t){\bigr )}}$

con alguna función vectorial adecuada que represente la aceleración dependiente de la posición. Normalmente, también se dan una posición inicial y una velocidad inicial . ${\displaystyle \mathbf {A} (\mathbf {x} )}$ ${\displaystyle \mathbf {x} (0)=\mathbf {x} _{0}}$ ${\displaystyle \mathbf {v} (0)={\dot {\mathbf {x} }}(0)=\mathbf {v} _{0}}$

Integración de Verlet (sin velocidades)

Para discretizar y resolver numéricamente este problema de valor inicial , se elige un paso de tiempo y se considera la secuencia de puntos de muestreo. La tarea consiste en construir una secuencia de puntos que sigan de cerca los puntos de la trayectoria de la solución exacta. ${\displaystyle \Delta t>0}$ ${\displaystyle t_{n}=n\,\Delta t}$ ${\displaystyle \mathbf {x} _{n}}$ ${\displaystyle \mathbf {x} (t_{n})}$

Mientras que el método de Euler utiliza la aproximación de diferencia hacia adelante a la primera derivada en ecuaciones diferenciales de orden uno, la integración de Verlet puede verse como el uso de la aproximación de diferencia central a la segunda derivada:

${\displaystyle {\begin{aligned}{\frac {\Delta ^{2}\mathbf {x} _{n}}{\Delta t^{2}}}&={\frac {{\frac {\mathbf {x} _{n+1}-\mathbf {x} _{n}}{\Delta t}}-{\frac {\mathbf {x} _{n}-\mathbf {x} _{n-1}}{\Delta t}}}{\Delta t}}\\[6pt]&={\frac {\mathbf {x} _{n+1}-2\mathbf {x} _{n}+\mathbf {x} _{n-1}}{\Delta t^{2}}}=\mathbf {a} _{n}=\mathbf {A} (\mathbf {x} _{n}).\end{aligned}}}$

La integración de Verlet en la forma utilizada como el método de Størmer ^[3] utiliza esta ecuación para obtener el siguiente vector de posición a partir de los dos anteriores sin utilizar la velocidad como

${\displaystyle {\begin{aligned}\mathbf {x} _{n+1}&=2\mathbf {x} _{n}-\mathbf {x} _{n-1}+\mathbf {a} _{n}\,\Delta t^{2},\\[6pt]\mathbf {a} _{n}&=\mathbf {A} (\mathbf {x} _{n}).\end{aligned}}}$

Error de discretización

La simetría temporal inherente al método reduce el nivel de errores locales introducidos en la integración por la discretización al eliminar todos los términos de grado impar, en este caso los términos de grado tres. El error local se cuantifica insertando los valores exactos en la iteración y calculando las expansiones de Taylor en el momento del vector de posición en diferentes direcciones temporales: ${\displaystyle \Delta t}$ ${\displaystyle \mathbf {x} (t_{n-1}),\mathbf {x} (t_{n}),\mathbf {x} (t_{n+1})}$ ${\displaystyle t=t_{n}}$ ${\displaystyle \mathbf {x} (t\pm \Delta t)}$

${\displaystyle {\begin{aligned}\mathbf {x} (t+\Delta t)&=\mathbf {x} (t)+\mathbf {v} (t)\Delta t+{\frac {\mathbf {a} (t)\Delta t^{2}}{2}}+{\frac {\mathbf {b} (t)\Delta t^{3}}{6}}+{\mathcal {O}}\left(\Delta t^{4}\right)\\\mathbf {x} (t-\Delta t)&=\mathbf {x} (t)-\mathbf {v} (t)\Delta t+{\frac {\mathbf {a} (t)\Delta t^{2}}{2}}-{\frac {\mathbf {b} (t)\Delta t^{3}}{6}}+{\mathcal {O}}\left(\Delta t^{4}\right),\end{aligned}}}$

donde es la posición, la velocidad, la aceleración y el tirón (tercera derivada de la posición con respecto al tiempo). ${\displaystyle \mathbf {x} }$ ${\displaystyle \mathbf {v} ={\dot {\mathbf {x} }}}$ ${\displaystyle \mathbf {a} ={\ddot {\mathbf {x} }}}$ ${\displaystyle \mathbf {b} ={\dot {\mathbf {a} }}={\overset {\dots }{\mathbf {x} }}}$

Al agregar estas dos expansiones se obtiene

${\displaystyle \mathbf {x} (t+\Delta t)=2\mathbf {x} (t)-\mathbf {x} (t-\Delta t)+\mathbf {a} (t)\Delta t^{2}+{\mathcal {O}}\left(\Delta t^{4}\right).}$

Podemos ver que los términos de primer y tercer orden de la expansión de Taylor se cancelan, lo que hace que el integrador de Verlet sea un orden más preciso que la integración mediante la simple expansión de Taylor únicamente.

Se debe tener cuidado con el hecho de que la aceleración aquí se calcula a partir de la solución exacta, , mientras que en la iteración se calcula en el punto de iteración central, . Al calcular el error global, es decir, la distancia entre la solución exacta y la secuencia de aproximación, esos dos términos no se cancelan exactamente, lo que influye en el orden del error global. ${\displaystyle \mathbf {a} (t)=\mathbf {A} {\bigl (}\mathbf {x} (t){\bigr )}}$ ${\displaystyle \mathbf {a} _{n}=\mathbf {A} (\mathbf {x} _{n})}$

Un ejemplo sencillo

Para comprender mejor la relación entre los errores locales y globales, resulta útil examinar ejemplos sencillos en los que la solución exacta, así como la solución aproximada, se pueden expresar en fórmulas explícitas. El ejemplo estándar para esta tarea es la función exponencial .

Considere la ecuación diferencial lineal con una constante . Sus soluciones base exactas son y . ${\displaystyle {\ddot {x}}(t)=w^{2}x(t)}$ ${\displaystyle w}$ ${\displaystyle e^{wt}}$ ${\displaystyle e^{-wt}}$

El método Størmer aplicado a esta ecuación diferencial conduce a una relación de recurrencia lineal

${\displaystyle x_{n+1}-2x_{n}+x_{n-1}=h^{2}w^{2}x_{n},}$

o

${\displaystyle x_{n+1}-2\left(1+{\tfrac {1}{2}}(wh)^{2}\right)x_{n}+x_{n-1}=0.}$

Se puede resolver hallando las raíces de su polinomio característico . Estas son ${\displaystyle q^{2}-2\left(1+{\tfrac {1}{2}}(wh)^{2}\right)q+1=0}$

${\displaystyle q_{\pm }=1+{\tfrac {1}{2}}(wh)^{2}\pm wh{\sqrt {1+{\tfrac {1}{4}}(wh)^{2}}}.}$

Las soluciones base de la recurrencia lineal son y . Para compararlas con las soluciones exactas, se calculan las expansiones de Taylor: ${\displaystyle x_{n}=q_{+}^{n}}$ ${\displaystyle x_{n}=q_{-}^{n}}$

${\displaystyle {\begin{aligned}q_{+}&=1+{\tfrac {1}{2}}(wh)^{2}+wh\left(1+{\tfrac {1}{8}}(wh)^{2}-{\tfrac {3}{128}}(wh)^{4}+{\mathcal {O}}\left(h^{6}\right)\right)\\&=1+(wh)+{\tfrac {1}{2}}(wh)^{2}+{\tfrac {1}{8}}(wh)^{3}-{\tfrac {3}{128}}(wh)^{5}+{\mathcal {O}}\left(h^{7}\right).\end{aligned}}}$

El cociente de esta serie con el de la exponencial comienza con , por lo que ${\displaystyle e^{wh}}$ ${\displaystyle 1-{\tfrac {1}{24}}(wh)^{3}+{\mathcal {O}}\left(h^{5}\right)}$

${\displaystyle {\begin{aligned}q_{+}&=\left(1-{\tfrac {1}{24}}(wh)^{3}+{\mathcal {O}}\left(h^{5}\right)\right)e^{wh}\\&=e^{-{\frac {1}{24}}(wh)^{3}+{\mathcal {O}}\left(h^{5}\right)}\,e^{wh}.\end{aligned}}}$

De allí se deduce que para la primera solución base el error se puede calcular como

${\displaystyle {\begin{aligned}x_{n}=q_{+}^{n}&=e^{-{\frac {1}{24}}(wh)^{2}\,wt_{n}+{\mathcal {O}}\left(h^{4}\right)}\,e^{wt_{n}}\\&=e^{wt_{n}}\left(1-{\tfrac {1}{24}}(wh)^{2}\,wt_{n}+{\mathcal {O}}(h^{4})\right)\\&=e^{wt_{n}}+{\mathcal {O}}\left(h^{2}t_{n}e^{wt_{n}}\right).\end{aligned}}}$

Es decir, aunque el error de discretización local es de orden 4, debido al segundo orden de la ecuación diferencial el error global es de orden 2, con una constante que crece exponencialmente en el tiempo.

Comenzando la iteración

Obsérvese que al comienzo de la iteración de Verlet en el paso , tiempo , cálculo , ya se necesita el vector de posición en el tiempo . A primera vista, esto podría dar problemas, porque las condiciones iniciales se conocen solo en el tiempo inicial . Sin embargo, a partir de ellas se conoce la aceleración , y se puede obtener una aproximación adecuada para la posición en el primer paso de tiempo utilizando el polinomio de Taylor de grado dos: ${\displaystyle n=1}$ ${\displaystyle t=t_{1}=\Delta t}$ ${\displaystyle \mathbf {x} _{2}}$ ${\displaystyle \mathbf {x} _{1}}$ ${\displaystyle t=t_{1}}$ ${\displaystyle t_{0}=0}$ ${\displaystyle \mathbf {a} _{0}=\mathbf {A} (\mathbf {x} _{0})}$

${\displaystyle \mathbf {x} _{1}=\mathbf {x} _{0}+\mathbf {v} _{0}\Delta t+{\tfrac {1}{2}}\mathbf {a} _{0}\Delta t^{2}\approx \mathbf {x} (\Delta t)+{\mathcal {O}}\left(\Delta t^{3}\right).}$

El error en el primer paso de tiempo es entonces del orden . Esto no se considera un problema porque en una simulación a lo largo de un gran número de pasos de tiempo, el error en el primer paso de tiempo es solo una cantidad insignificantemente pequeña del error total, que en el momento es del orden , tanto para la distancia de los vectores de posición a como para la distancia de las diferencias divididas a . Además, para obtener este error global de segundo orden, el error inicial debe ser al menos de tercer orden. ${\displaystyle {\mathcal {O}}\left(\Delta t^{3}\right)}$ ${\displaystyle t_{n}}$ ${\displaystyle {\mathcal {O}}\left(e^{Lt_{n}}\Delta t^{2}\right)}$ ${\displaystyle \mathbf {x} _{n}}$ ${\displaystyle \mathbf {x} (t_{n})}$ ${\displaystyle {\tfrac {\mathbf {x} _{n+1}-\mathbf {x} _{n}}{\Delta t}}}$ ${\displaystyle {\tfrac {\mathbf {x} (t_{n+1})-\mathbf {x} (t_{n})}{\Delta t}}}$

Diferencias horarias no constantes

Una desventaja del método de Størmer-Verlet es que si el paso de tiempo ( ) cambia, el método no aproxima la solución de la ecuación diferencial. Esto se puede corregir utilizando la fórmula ^[4] ${\displaystyle \Delta t}$

${\displaystyle \mathbf {x} _{i+1}=\mathbf {x} _{i}+\left(\mathbf {x} _{i}-\mathbf {x} _{i-1}\right){\frac {\Delta t_{i}}{\Delta t_{i-1}}}+\mathbf {a} _{i}\Delta t_{i}^{2}.}$

Una derivación más exacta utiliza la serie de Taylor (de segundo orden) para tiempos y para obtener después de la eliminación de ${\displaystyle t_{i}}$ ${\displaystyle t_{i+1}=t_{i}+\Delta t_{i}}$ ${\displaystyle t_{i-1}=t_{i}-\Delta t_{i-1}}$ ${\displaystyle \mathbf {v} _{i}}$

${\displaystyle {\frac {\mathbf {x} _{i+1}-\mathbf {x} _{i}}{\Delta t_{i}}}+{\frac {\mathbf {x} _{i-1}-\mathbf {x} _{i}}{\Delta t_{i-1}}}=\mathbf {a} _{i}\,{\frac {\Delta t_{i}+\Delta t_{i-1}}{2}},}$

de modo que la fórmula de iteración se convierte en

${\displaystyle \mathbf {x} _{i+1}=\mathbf {x} _{i}+(\mathbf {x} _{i}-\mathbf {x} _{i-1}){\frac {\Delta t_{i}}{\Delta t_{i-1}}}+\mathbf {a} _{i}\,{\frac {\Delta t_{i}+\Delta t_{i-1}}{2}}\,\Delta t_{i}.}$

Cálculo de velocidades: método Størmer-Verlet

Las velocidades no se dan explícitamente en la ecuación básica de Størmer, pero a menudo son necesarias para el cálculo de ciertas cantidades físicas como la energía cinética. Esto puede crear desafíos técnicos en las simulaciones de dinámica molecular , porque la energía cinética y las temperaturas instantáneas en el tiempo no se pueden calcular para un sistema hasta que se conocen las posiciones en el tiempo . Esta deficiencia se puede solucionar utilizando el algoritmo de velocidad de Verlet o estimando la velocidad utilizando los términos de posición y el teorema del valor medio : ${\displaystyle t}$ ${\displaystyle t+\Delta t}$

${\displaystyle \mathbf {v} (t)={\frac {\mathbf {x} (t+\Delta t)-\mathbf {x} (t-\Delta t)}{2\Delta t}}+{\mathcal {O}}\left(\Delta t^{2}\right).}$

Nótese que este término de velocidad está un paso por detrás del término de posición, ya que esto es para la velocidad en el tiempo , no , lo que significa que es una aproximación de segundo orden a . Con el mismo argumento, pero dividiendo a la mitad el paso de tiempo, es una aproximación de segundo orden a , con . ${\displaystyle t}$ ${\displaystyle t+\Delta t}$ ${\displaystyle \mathbf {v} _{n}={\tfrac {\mathbf {x} _{n+1}-\mathbf {x} _{n-1}}{2\Delta t}}}$ ${\displaystyle \mathbf {v} (t_{n})}$ ${\displaystyle \mathbf {v} _{n+{\frac {1}{2}}}={\tfrac {\mathbf {x} _{n+1}-\mathbf {x} _{n}}{\Delta t}}}$ ${\displaystyle \mathbf {v} \left(t_{n+{\frac {1}{2}}}\right)}$ ${\displaystyle t_{n+{\frac {1}{2}}}=t_{n}+{\tfrac {1}{2}}\Delta t}$

Se puede acortar el intervalo para aproximar la velocidad en el tiempo a costa de la precisión: ${\displaystyle t+\Delta t}$

${\displaystyle \mathbf {v} (t+\Delta t)={\frac {\mathbf {x} (t+\Delta t)-\mathbf {x} (t)}{\Delta t}}+{\mathcal {O}}(\Delta t).}$

Verlet de velocidad

Un algoritmo relacionado y más comúnmente utilizado es el algoritmo de velocidad de Verlet, ^[5] similar al método de salto de rana , excepto que la velocidad y la posición se calculan con el mismo valor de la variable de tiempo (el método de salto de rana no lo hace, como sugiere el nombre). Este utiliza un enfoque similar, pero incorpora explícitamente la velocidad, lo que resuelve el problema del primer paso de tiempo en el algoritmo básico de Verlet:

${\displaystyle {\begin{aligned}\mathbf {x} (t+\Delta t)&=\mathbf {x} (t)+\mathbf {v} (t)\,\Delta t+{\tfrac {1}{2}}\,\mathbf {a} (t)\Delta t^{2},\\[6pt]\mathbf {v} (t+\Delta t)&=\mathbf {v} (t)+{\frac {\mathbf {a} (t)+\mathbf {a} (t+\Delta t)}{2}}\Delta t.\end{aligned}}}$

Se puede demostrar que el error en el algoritmo Verlet de velocidad es del mismo orden que en el algoritmo Verlet básico. Nótese que el algoritmo de velocidad no necesariamente consume más memoria, porque en el algoritmo Verlet básico hacemos un seguimiento de dos vectores de posición, mientras que en el algoritmo Verlet de velocidad hacemos un seguimiento de un vector de posición y un vector de velocidad. El esquema de implementación estándar de este algoritmo es:

1. Calcular . ${\displaystyle \mathbf {v} \left(t+{\tfrac {1}{2}}\,\Delta t\right)=\mathbf {v} (t)+{\tfrac {1}{2}}\,\mathbf {a} (t)\,\Delta t}$
2. Calcular . ${\displaystyle \mathbf {x} (t+\Delta t)=\mathbf {x} (t)+\mathbf {v} \left(t+{\tfrac {1}{2}}\,\Delta t\right)\,\Delta t}$
3. Derive del potencial de interacción utilizando . ${\displaystyle \mathbf {a} (t+\Delta t)}$ ${\displaystyle \mathbf {x} (t+\Delta t)}$
4. Calcular . ${\displaystyle \mathbf {v} (t+\Delta t)=\mathbf {v} \left(t+{\tfrac {1}{2}}\,\Delta t\right)+{\tfrac {1}{2}}\,\mathbf {a} (t+\Delta t)\Delta t}$

Este algoritmo también funciona con pasos de tiempo variables y es idéntico a la forma 'kick-drift-kick' del método de integración de salto.

Eliminando la velocidad de medio paso, este algoritmo se puede acortar a

1. Calcular . ${\displaystyle \mathbf {x} (t+\Delta t)=\mathbf {x} (t)+\mathbf {v} (t)\,\Delta t+{\tfrac {1}{2}}\,\mathbf {a} (t)\,\Delta t^{2}}$
2. Derive del potencial de interacción utilizando . ${\displaystyle \mathbf {a} (t+\Delta t)}$ ${\displaystyle \mathbf {x} (t+\Delta t)}$
3. Calcular . ${\displaystyle \mathbf {v} (t+\Delta t)=\mathbf {v} (t)+{\tfrac {1}{2}}\,{\bigl (}\mathbf {a} (t)+\mathbf {a} (t+\Delta t){\bigr )}\Delta t}$

Tenga en cuenta, sin embargo, que este algoritmo supone que la aceleración solo depende de la posición y no de la velocidad . ${\displaystyle \mathbf {a} (t+\Delta t)}$ ${\displaystyle \mathbf {x} (t+\Delta t)}$ ${\displaystyle \mathbf {v} (t+\Delta t)}$

Se podría observar que los resultados a largo plazo del método de Verlet de velocidad, y de manera similar del método de salto de rana, son un orden mejor que los del método de Euler semiimplícito . Los algoritmos son casi idénticos hasta un cambio de medio paso de tiempo en la velocidad. Esto se puede demostrar girando el bucle anterior para comenzar en el paso 3 y luego notando que el término de aceleración en el paso 1 podría eliminarse combinando los pasos 2 y 4. La única diferencia es que la velocidad del punto medio en el método de Verlet de velocidad se considera la velocidad final en el método de Euler semiimplícito.

El error global de todos los métodos de Euler es de orden uno, mientras que el error global de este método es, de manera similar al método del punto medio , de orden dos. Además, si la aceleración resulta de hecho de las fuerzas en un sistema mecánico conservativo o hamiltoniano , la energía de la aproximación oscila esencialmente alrededor de la energía constante del sistema resuelto exactamente, con un error global acotado nuevamente de orden uno para Euler semiexplícito y de orden dos para salto de Verlet. Lo mismo ocurre con todas las demás cantidades conservadas del sistema, como el momento lineal o angular, que siempre se conservan o casi se conservan en un integrador simpléctico . ^[6]

El método de Verlet de velocidad es un caso especial del método Newmark-beta con y . ${\displaystyle \beta =0}$ ${\displaystyle \gamma ={\tfrac {1}{2}}}$

Representación algorítmica

Dado que el algoritmo de Verlet de velocidad es un algoritmo útil en general en aplicaciones 3D, una solución escrita en C++ podría verse como la siguiente. Este tipo de integración de posición aumentará significativamente la precisión en simulaciones y juegos 3D en comparación con el método de Euler habitual.

Estructura Cuerpo { Vec3d pos { 0.0 , 0.0 , 0.0 };       Vec3d vel { 2.0 , 0.0 , 0.0 }; // 2 m/s a lo largo del eje x        Vec3d acc { 0.0 , 0.0 , 0.0 }; // sin aceleración al principio        doble masa = 1,0 ; // 1 kg     /** * Actualiza pos y vel mediante la integración "Velocity Verlet" * @param dt DeltaTime / paso de tiempo [por ejemplo: 0,01] */ actualización nula ( doble dt )   { Vec3d nueva_pos = pos + vel * dt + acc * ( dt * dt * 0.5 );        Vec3d new_acc = aplicar_fuerzas ();    Vec3d nueva_vel = vel + ( acc + nueva_acc ) * ( dt * 0.5 );      pos = nueva_pos ;   vel = nuevo_vel ;   acc = nuevo_acc ;   } /** * Para aplicar velocidad a sus objetos, calcule el vector de fuerza requerido. * y aplicar las fuerzas acumuladas aquí. */ Vec3d aplicar_fuerzas () constante   { Vec3d new_acc = Vec3d { 0.0 , 0.0 , -9.81 }; // 9.81 m/s² hacia abajo en el eje z        //Aplica cualquier otra fuerza aquí... // NOTA: Evite depender de `vel` porque Velocity Verlet asume que la aceleración depende de la posición. devolver nuevo_acc ;  }};

Términos de error

El error de truncamiento global del método de Verlet es , tanto para la posición como para la velocidad. Esto contrasta con el hecho de que el error local en la posición es solo como se describió anteriormente. La diferencia se debe a la acumulación del error de truncamiento local a lo largo de todas las iteraciones. ${\displaystyle {\mathcal {O}}\left(\Delta t^{2}\right)}$ ${\displaystyle {\mathcal {O}}\left(\Delta t^{4}\right)}$

El error global se puede derivar observando lo siguiente:

${\displaystyle \operatorname {error} {\bigl (}x(t_{0}+\Delta t){\bigr )}={\mathcal {O}}\left(\Delta t^{4}\right)}$

y

${\displaystyle x(t_{0}+2\Delta t)=2x(t_{0}+\Delta t)-x(t_{0})+\Delta t^{2}{\ddot {x}}(t_{0}+\Delta t)+{\mathcal {O}}\left(\Delta t^{4}\right).}$

Por lo tanto

${\displaystyle \operatorname {error} {\bigl (}x(t_{0}+2\Delta t){\bigr )}=2\cdot \operatorname {error} {\bigl (}x(t_{0}+\Delta t){\bigr )}+{\mathcal {O}}\left(\Delta t^{4}\right)=3\,{\mathcal {O}}\left(\Delta t^{4}\right).}$

Similarmente:

${\displaystyle {\begin{aligned}\operatorname {error} {\bigl (}x(t_{0}+3\Delta t){\bigl )}&=6\,{\mathcal {O}}\left(\Delta t^{4}\right),\\[6px]\operatorname {error} {\bigl (}x(t_{0}+4\Delta t){\bigl )}&=10\,{\mathcal {O}}\left(\Delta t^{4}\right),\\[6px]\operatorname {error} {\bigl (}x(t_{0}+5\Delta t){\bigl )}&=15\,{\mathcal {O}}\left(\Delta t^{4}\right),\end{aligned}}}$

lo cual se puede generalizar a (se puede demostrar por inducción, pero se da aquí sin prueba):

${\displaystyle \operatorname {error} {\bigl (}x(t_{0}+n\Delta t){\bigr )}={\frac {n(n+1)}{2}}\,{\mathcal {O}}\left(\Delta t^{4}\right).}$

Si consideramos el error global de posición entre y , donde , está claro que ^{[ cita requerida ]} ${\displaystyle x(t)}$ ${\displaystyle x(t+T)}$ ${\displaystyle T=n\Delta t}$

${\displaystyle \operatorname {error} {\bigl (}x(t_{0}+T){\bigr )}=\left({\frac {T^{2}}{2\Delta t^{2}}}+{\frac {T}{2\Delta t}}\right){\mathcal {O}}\left(\Delta t^{4}\right),}$

y por lo tanto, el error global (acumulativo) durante un intervalo de tiempo constante viene dado por

${\displaystyle \operatorname {error} {\bigr (}x(t_{0}+T){\bigl )}={\mathcal {O}}\left(\Delta t^{2}\right).}$

Dado que la velocidad se determina de forma no acumulativa a partir de las posiciones en el integrador de Verlet, el error global en la velocidad también es . ${\displaystyle {\mathcal {O}}\left(\Delta t^{2}\right)}$

En las simulaciones de dinámica molecular, el error global suele ser mucho más importante que el error local y, por lo tanto, el integrador de Verlet se conoce como un integrador de segundo orden.

Restricciones

Los sistemas de múltiples partículas con restricciones son más sencillos de resolver con la integración de Verlet que con los métodos de Euler. Las restricciones entre puntos pueden ser, por ejemplo, potenciales que los limitan a una distancia específica o fuerzas de atracción. Pueden modelarse como resortes que conectan las partículas. Utilizando resortes de rigidez infinita, el modelo puede entonces resolverse con un algoritmo de Verlet.

En una dimensión, la relación entre las posiciones sin restricciones y las posiciones reales de los puntos en el tiempo , dada una distancia de restricción deseada de , se puede encontrar con el algoritmo ${\displaystyle {\tilde {x}}_{i}^{(t)}}$ ${\displaystyle x_{i}^{(t)}}$ ${\displaystyle i}$ ${\displaystyle t}$ ${\displaystyle r}$

${\displaystyle {\begin{aligned}d_{1}&=x_{2}^{(t)}-x_{1}^{(t)},\\[6px]d_{2}&=\|d_{1}\|,\\[6px]d_{3}&={\frac {d_{2}-r}{d_{2}}},\\[6px]x_{1}^{(t+\Delta t)}&={\tilde {x}}_{1}^{(t+\Delta t)}+{\tfrac {1}{2}}d_{1}d_{3},\\[6px]x_{2}^{(t+\Delta t)}&={\tilde {x}}_{2}^{(t+\Delta t)}-{\tfrac {1}{2}}d_{1}d_{3}.\end{aligned}}}$

La integración de Verlet es útil porque relaciona directamente la fuerza con la posición, en lugar de resolver el problema utilizando velocidades.

Sin embargo, surgen problemas cuando actúan múltiples fuerzas restrictivas sobre cada partícula. Una forma de resolver esto es recorrer cada punto de una simulación, de modo que en cada punto ya se utilice la relajación de la restricción del último para acelerar la propagación de la información. En una simulación, esto se puede implementar utilizando pequeños pasos de tiempo para la simulación, utilizando un número fijo de pasos de resolución de restricciones por paso de tiempo o resolviendo las restricciones hasta que se cumplan con una desviación específica.

Al aproximar las restricciones localmente a primer orden, esto es lo mismo que el método de Gauss-Seidel . Para matrices pequeñas se sabe que la descomposición LU es más rápida. Los sistemas grandes se pueden dividir en grupos (por ejemplo, cada ragdoll  = grupo). Dentro de los grupos se utiliza el método LU, entre grupos se utiliza el método de Gauss-Seidel . El código de la matriz se puede reutilizar: la dependencia de las fuerzas en las posiciones se puede aproximar localmente a primer orden, y la integración de Verlet se puede hacer más implícita.

Existen programas sofisticados, como SuperLU ^[7], que permiten resolver problemas complejos utilizando matrices dispersas. Se pueden utilizar técnicas específicas, como el uso de (grupos de) matrices, para abordar problemas específicos, como el de la fuerza que se propaga a través de una lámina de tela sin formar una onda sonora . ^[8]

Otra forma de resolver restricciones holonómicas es utilizar algoritmos de restricción .

Reacciones de colisión

Una forma de reaccionar ante las colisiones es utilizar un sistema basado en penalizaciones, que básicamente aplica una fuerza determinada a un punto en el momento del contacto. El problema con esto es que es muy difícil elegir la fuerza aplicada. Si se utiliza una fuerza demasiado fuerte, los objetos se volverán inestables, y si es demasiado débil, se penetrarán entre sí. Otra forma es utilizar reacciones de colisión por proyección, que toman el punto ofensivo e intentan moverlo la distancia más corta posible para alejarlo del otro objeto.

En el último caso, la integración de Verlet manejaría automáticamente la velocidad impartida por la colisión; sin embargo, tenga en cuenta que no se garantiza que esto se haga de una manera que sea coherente con la física de colisiones (es decir, no se garantiza que los cambios en el momento sean realistas). En lugar de cambiar implícitamente el término de velocidad, sería necesario controlar explícitamente las velocidades finales de los objetos que colisionan (cambiando la posición registrada del paso de tiempo anterior).

Los dos métodos más simples para decidir una nueva velocidad son las colisiones perfectamente elásticas e inelásticas . Una estrategia un poco más complicada que ofrece más control implicaría utilizar el coeficiente de restitución .

Véase también

• Condición de Courant-Friedrichs-Lewy
• Deriva energética
• Integrador simpléctico
• Integración de Leapfrog
• Algoritmo de Beeman

Literatura

1. Verlet, Loup (1967). "Experimentos" informáticos sobre fluidos clásicos. I. Propiedades termodinámicas de las moléculas de Lennard−Jones". Physical Review . 159 (1): 98–103. Bibcode :1967PhRv..159...98V. doi : 10.1103/PhysRev.159.98 .
2. Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). "Sección 17.4. Ecuaciones conservativas de segundo orden". Recetas numéricas: el arte de la computación científica (3.ª ed.). Nueva York: Cambridge University Press. ISBN 978-0-521-88068-8.
3. Página web archivada el 3 de agosto de 2004 en Wayback Machine con una descripción del método Størmer.
4. Dummer, Jonathan. "Un método de integración de Verlet simple con corrección temporal".
5. Swope, William C.; HC Andersen; PH Berens; KR Wilson (1 de enero de 1982). "Un método de simulación por ordenador para el cálculo de constantes de equilibrio para la formación de agrupaciones físicas de moléculas: aplicación a agrupaciones pequeñas de agua". The Journal of Chemical Physics . 76 (1): 648 (Apéndice). Bibcode :1982JChPh..76..637S. doi :10.1063/1.442716.
6. Peluquero, Ernst; Lubich, cristiano; Wanner, Gerhard (2003). "Integración numérica geométrica ilustrada por el método Størmer/Verlet". Acta Numérica . 12 : 399–450. Código Bib : 2003AcNum..12..399H. CiteSeerX 10.1.1.7.7106 . doi :10.1017/S0962492902000144. S2CID  122016794. 
7. Guía del usuario de SuperLU.
8. Baraff, D.; Witkin, A. (1998). "Grandes pasos en la simulación de telas" (PDF) . Actas de Computer Graphics . Serie de conferencias anuales: 43–54.

Enlaces externos

• Demostración de integración de Verlet y código como subprograma Java
• Física avanzada de personajes por Thomas Jakobsen
• Simulaciones de la teoría de la dinámica molecular – parte inferior de la página
• Integración de Verlet implementada en JavaScript moderno (parte inferior de la página)