Integración 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 del tiempo y la preservación de la forma simpléctica en el espacio de fases , sin ningún costo computacional adicional significativo sobre el método simple de Euler .

Störmer básico–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 conservadores 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 la hora,
• ${\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 varias opciones de la función potencial , se puede utilizar para describir la evolución de diversos sistemas físicos, desde el movimiento de moléculas que interactúan 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 en 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 en diferencias directas a la primera derivada en ecuaciones diferenciales de orden uno, se puede considerar que la integración de Verlet utiliza la aproximación en diferencias centrales 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 método de Störmer ^[3] usa esta ecuación para obtener el siguiente vector de posición de los dos anteriores sin usar 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, aquí 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 está 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} }}}$

Agregar estas dos expansiones da

${\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 una simple expansión de Taylor únicamente.

Se debe tener precaución 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 central de la iteración . 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, pueden expresarse 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 de base exacta 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 encontrando las raíces de su polinomio característico . Estos 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 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 empieza por , así 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 ahí se deduce que para la solución de primera 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

Tenga en cuenta que al inicio de la iteración de Verlet en paso , tiempo y cálculo , ya se necesita el vector de posición en tiempo . A primera vista, esto podría plantear problemas, porque las condiciones iniciales sólo se conocen en el momento inicial . Sin embargo, a partir de estos 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 del tiempo es entonces de orden . Esto no se considera un problema porque en una simulación de un gran número de pasos de tiempo, el error en el primer paso de tiempo es sólo una cantidad insignificante del error total, que en el tiempo es del orden , tanto para la distancia del vectores de posición para la distancia de las diferencias divididas para . 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 Störmer-Verlet es que si el paso de tiempo ( ) cambia, el método no aproxima la solución a la ecuación diferencial. Esto se puede corregir usando 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) en 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}},}$

para que la fórmula de iteración se convierta 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 de 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 conozcan las posiciones en el momento . 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).}$

Tenga en cuenta que este término de velocidad está un paso por detrás del término de posición, ya que 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 reduciendo a la mitad el paso de tiempo, se obtiene 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).}$

Velocidad Verlet

Un algoritmo relacionado, y más comúnmente utilizado, es el algoritmo de velocidad Verlet , ^[5] similar al método de salto , excepto que la velocidad y la posición se calculan con el mismo valor de la variable de tiempo (el salto no lo hace, como sugiere el nombre). . Esto utiliza un enfoque similar, pero incorpora explícitamente la velocidad, resolviendo 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 Verlet de velocidad es del mismo orden que en el Verlet básico. Tenga en cuenta que el algoritmo de velocidad no consume necesariamente más memoria porque, en Verlet básico, realizamos un seguimiento de dos vectores de posición, mientras que en Verlet de velocidad, realizamos 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 usando . ${\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" de integración del método de salto .

Eliminando la velocidad de medio paso, este algoritmo puede acortarse 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 usando . ${\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 sólo 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 señalar que los resultados a largo plazo de la velocidad de Verlet y, de manera similar, del salto son un orden mejores que el método semiimplícito de Euler . Los algoritmos son casi idénticos hasta un cambio de medio paso de tiempo en la velocidad. Esto se puede probar girando el bucle anterior para comenzar en el paso 3 y luego notar 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 la velocidad Verlet 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, similar al método del punto medio , de orden dos. Además, si la aceleración realmente resulta de las fuerzas en un sistema mecánico conservador o hamiltoniano , la energía de la aproximación oscila esencialmente alrededor de la energía constante del sistema exactamente resuelto, con un error global nuevamente limitado de orden uno para Euler y semiexplícitos. Pida dos para Verlet-leapfrog. 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 velocidad Verlet 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 la velocidad Verlet es un algoritmo generalmente útil en aplicaciones 3D, una solución general escrita en C++ podría verse como a continuación. Se utiliza una fuerza de arrastre simplificada para demostrar el cambio en la aceleración; sin embargo, solo es necesaria si la aceleración no es constante.

estructura del 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 cuenta { 0.0 , 0.0 , 0.0 }; // sin aceleración al principio        doble masa = 1,0 ; // 1 kg     doble arrastre = 0,1 ; // rho*C*Area – arrastre simplificado para este ejemplo     /** * Actualice pos y vel usando la integración "Velocity Verlet" * @param dt DeltaTime / paso de tiempo [por ejemplo: 0,01] */ actualización nula ( doble dt )   { Vec3d new_pos = pos + vel * dt + acc * ( dt * dt * 0.5 );        Vec3d new_acc = aplicar_fuerzas (); // sólo es necesario si la aceleración no es constante     Vec3d nueva_vel = vel + ( acc + nueva_acc ) * ( dt * 0.5 );      pos = nueva_pos ;   vel = nueva_vel ;   cuenta = nueva_cuenta ;   } Vec3d apply_forces () const   { Vec3d grav_acc = Vec3d { 0,0 , 0,0 , -9,81 }; // 9,81 m/s² hacia abajo en el eje z        Vec3d drag_force = 0.5 * arrastrar * ( vel * vel ); // D = 0,5 * (rho * C * Área * vel^2)           Vec3d drag_acc = fuerza_arrastre / masa ; // a = F/m       devolver grav_acc - drag_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 sólo el descrito anteriormente. La diferencia se debe a la acumulación del error de truncamiento local en 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}}}$

que 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 en la posición entre y , donde , está claro que ^{[ cita necesaria ]} ${\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 (acumulado) 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).}$

Debido a 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 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 resolverse con un algoritmo de Verlet.

En una dimensión, la relación entre las posiciones no restringidas 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 usando velocidades.

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

Cuando se aproximan las restricciones localmente al 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 muñeco de trapo  = grupo). Dentro de los conglomerados se utiliza el método LU, entre conglomerados se utiliza el método de Gauss-Seidel . El código de la matriz se puede reutilizar: la dependencia de las fuerzas de las posiciones se puede aproximar localmente al primer orden y la integración de Verlet se puede hacer más implícita.

Existe software sofisticado, como SuperLU ^[7] , para resolver problemas complejos utilizando matrices dispersas. Se pueden utilizar técnicas específicas, como el uso de (grupos de) matrices, para abordar el problema específico, como el de la fuerza que se propaga a través de una hoja 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 impartida. Utilice una fuerza demasiado fuerte y los objetos se volverán inestables, demasiado débiles y se penetrarán entre sí. Otra forma es utilizar reacciones de colisión de proyección, que toma el punto infractor e intenta moverlo la distancia más corta posible para sacarlo del otro objeto.

La integración de Verlet manejaría automáticamente la velocidad impartida por la colisión en el último caso; sin embargo, tenga en cuenta que no se garantiza que esto sea coherente con la física de colisiones (es decir, no se garantiza que los cambios en el impulso 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 chocan (cambiando la posición registrada desde el paso de tiempo anterior).

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

Ver también

• Condición de Courant-Friedrichs-Lewy
• Deriva de energía
• Integrador simpléctico
• Integración avanzada
• 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". Revisión física . 159 (1): 98-103. Código Bib : 1967PhRv..159...98V. doi : 10.1103/PhysRev.159.98 .
2. Prensa, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). "Sección 17.4. Ecuaciones conservadoras de segundo orden". Recetas numéricas: el arte de la informática 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. Más tonto, Jonathan. "Un método simple de integración Verlet con corrección de tiempo".
5. Golpe, William C.; HC Andersen; PH Berens; KR Wilson (1 de enero de 1982). "Un método de simulación por computadora para el cálculo de constantes de equilibrio para la formación de grupos físicos de moléculas: aplicación a pequeños grupos de agua". La Revista de Física Química . 76 (1): 648 (Apéndice). Código bibliográfico : 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 gráficos por computadora . Serie de conferencias anuales: 43–54.

enlaces externos

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