La integración de Verlet ( pronunciación francesa: [vɛʁlɛ] ) es un método numérico utilizado para integrar de Newton ecuaciones de movimiento . [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 Delambre y ha sido redescubierto muchas veces desde entonces, más recientemente por Loup Verlet en la década de 1960 para su uso en dinámica molecular. También fue utilizado por Cowell y Crommelin en 1909 para calcular la órbita del cometa Halley , y porCarl 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 Verlet proporciona una buena estabilidad numérica , así como otras propiedades que son importantes en sistemas físicos como la reversibilidad en el tiempo y la preservación de la forma simpléctica en el espacio de fase , sin costo computacional adicional significativo sobre el método de Euler simple .
Störmer – Verlet básico
Para una ecuación diferencial de segundo orden del tipo con condiciones iniciales y , una solución numérica aproximada en los tiempos con tamaño de paso se puede obtener mediante el siguiente método:
- colocar ,
- para n = 1, 2, ... iterar
Ecuaciones de movimiento
La ecuación de movimiento de Newton para sistemas físicos conservadores es
o individualmente
dónde
- t es el tiempo,
- es el conjunto del vector de posición de N objetos,
- V es la función de potencial escalar,
- F es el gradiente negativo del potencial , que da el conjunto de fuerzas sobre las partículas,
- M es la matriz de masa , típicamente diagonal con bloques con masa para cada partícula.
Esta ecuación, para varias elecciones de la función potencial V , 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 .
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
con alguna función de valor vectorial adecuada A que represente la aceleración dependiente de la posición. Normalmente, una posición inicial y una velocidad inicial también se dan.
Integración Verlet (sin velocidades)
Para discretizar y resolver numéricamente este problema de valor inicial , un paso de tiempo se elige, y la secuencia de puntos de muestreo considerado. La tarea es construir una secuencia de puntos. que siguen de cerca los puntos en la trayectoria de la solución exacta.
Donde el método de Euler usa la aproximación de diferencia directa a la primera derivada en ecuaciones diferenciales de orden uno, se puede considerar que la integración de Verlet usa la aproximación de diferencia central a la segunda derivada:
La integración de Verlet en la forma utilizada como el 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
Error de discretización
La simetría de tiempo 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 en de grado tres. El error local se cuantifica insertando los valores exactosen la iteración y el cálculo de las expansiones de Taylor en el tiempo del vector de posición en diferentes direcciones de tiempo:
dónde es la posición, la velocidad, la aceleración, y el tirón (tercera derivada de la posición con respecto al tiempo).
Agregar estas dos expansiones da
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 por simple expansión de Taylor.
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, que es 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.
Un simple ejemplo
Para comprender mejor la relación de los errores locales y globales, es útil examinar ejemplos simples 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 w . Sus soluciones base exactas son y .
El método de Störmer aplicado a esta ecuación diferencial conduce a una relación de recurrencia lineal
o
Se puede resolver encontrando las raíces de su polinomio característico. . Estos son
Las soluciones de base de la recurrencia lineal son y . Para compararlos con las soluciones exactas, se calculan las expansiones de Taylor:
El cociente de esta serie con el de la exponencial comienza con , entonces
De ahí se deduce que para la primera solución básica, el error se puede calcular como
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 comienzo de la iteración de Verlet en el paso , hora , informática , uno ya necesita el vector de posición en el momento . A primera vista, esto podría dar problemas, porque las condiciones iniciales se conocen solo en el momento inicial.. Sin embargo, de estos la aceleracinse conoce, 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:
El error en el primer paso de tiempo es entonces de orden . Esto no se considera un problema porque en una simulación sobre un gran número de pasos de tiempo, el error en el primer paso de tiempo es solo una cantidad insignificante del error total, que en el momento es de la orden , tanto para la distancia de los vectores de posición a en cuanto a la distancia de las diferencias divididas a . Además, para obtener este error global de segundo orden, el error inicial debe ser de al menos tercer orden.
Diferencias de tiempo 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 mediante la fórmula [4]
Una derivación más exacta usa la serie de Taylor (de segundo orden) en por tiempos y para obtener después de la eliminación de
para que la fórmula de iteración se convierta en
Calcular 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 puede calcular para un sistema hasta que se conozcan las posiciones en el momento . Esta deficiencia puede tratarse utilizando el algoritmo de Verlet de velocidad o estimando la velocidad utilizando los términos de posición y el teorema del valor medio :
Tenga en cuenta que este término de velocidad es un paso por detrás del término de posición, ya que esto es para la velocidad en el tiempo , no , significa que es una aproximación de segundo orden a . Con el mismo argumento, pero reduciendo a la mitad el paso de tiempo, es una aproximación de segundo orden a , con .
Se puede acortar el intervalo para aproximar la velocidad en el tiempo a costa de la precisión:
Verlet de velocidad
Un algoritmo relacionado, y de uso más común, es el algoritmo de velocidad Verlet , [5] similar al método de salto , excepto que la velocidad y la posición se calculan al mismo valor de la variable de tiempo (el salto no lo hace, como sugiere el nombre) . Esto usa un enfoque similar, pero incorpora explícitamente la velocidad, resolviendo el problema del primer paso de tiempo en el algoritmo básico de Verlet:
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, hacemos un seguimiento de dos vectores de posición, mientras que en 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:
- Calcular .
- Calcular .
- Derivar del potencial de interacción usando .
- Calcular .
Eliminando la velocidad de medio paso, este algoritmo puede acortarse a
- Calcular .
- Derivar del potencial de interacción usando .
- Calcular .
Sin embargo, tenga en cuenta que este algoritmo asume que la aceleración solo depende de la posición y no depende de la velocidad .
Se podría notar que los resultados a largo plazo de la velocidad Verlet y, de manera similar, del salto de rana son un orden mejores que el método de Euler semi-implícito . Los algoritmos son casi idénticos hasta un cambio de medio paso de tiempo en la velocidad. Esto se prueba fácilmente girando el ciclo anterior para comenzar en el paso 3 y luego observando 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 semi-implí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 de hecho 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 resuelto exactamente, con un error global ligado nuevamente de orden uno para Euler y semi-explícito. ordene 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 .
Representación algorítmica
Dado que velocity 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 cuerpo{ Vec3d pos { 0.0 , 0.0 , 0.0 }; Vec3d vel { 2.0 , 0.0 , 0.0 }; // 2 m / sa lo largo del eje x Vec3d acc { 0.0 , 0.0 , 0.0 }; // sin aceleración al principio masa doble = 1,0 ; // 1 kilogramo doble arrastre = 0,1 ; // rho * C * Area - arrastre simplificado para este ejemplo / ** * Actualizar 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 = apply_forces (); // solo es necesario si la aceleración no es constante Vec3d nuevo_nivel = vel + ( acc + nuevo_acc ) * ( dt * 0.5 ); pos = new_pos ; vel = nuevo_nivel ; acc = new_acc ; } Vec3d apply_forces () const { Vec3d grav_acc = Vec3d { 0.0 , 0.0 , -9.81 }; // 9,81 m / s ^ 2 hacia abajo en el eje Z Vec3d drag_force = 0.5 * drag * ( vel * abs ( vel )); // D = 0.5 * (rho * C * Área * vel ^ 2) Vec3d drag_acc = drag_force / mass ; // a = F / m return grav_acc - drag_acc ; }};
Términos de error
El error local en la posición del integrador Verlet es como se describió anteriormente, y el error local en la velocidad es .
El error global en la posición, por el contrario, es , y el error global en la velocidad es . Estos se pueden derivar observando lo siguiente:
y
Por lo tanto
Similar:
que se puede generalizar a (se puede demostrar por inducción, pero se da aquí sin prueba):
Si consideramos el error global en la posición entre y , dónde , está claro que
- [ cita requerida ]
y por lo tanto, el error global (acumulativo) sobre un intervalo de tiempo constante viene dado por
Debido a que la velocidad se determina de manera no acumulativa a partir de las posiciones en el integrador Verlet, el error global en la velocidad también es .
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 partículas múltiples con restricciones son más simples 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 restringen a una distancia específica o fuerzas de atracción. Pueden modelarse como resortes que conectan las partículas. Usando resortes de rigidez infinita, el modelo puede resolverse con un algoritmo de Verlet.
En una dimensión, la relación entre las posiciones ilimitadas y las posiciones reales de los puntos i en el tiempo t se puede encontrar con el algoritmo
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 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 use 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.
Al aproximar 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 de LU es más rápida. Los sistemas grandes se pueden dividir en grupos (por ejemplo, cada ragdoll = cluster). Dentro de los clústeres se utiliza el método LU, entre los clústeres se utiliza el método Gauss-Seidel . El código de la matriz se puede reutilizar: la dependencia de las fuerzas en las posiciones se puede aproximar localmente al primer orden, y la integración de Verlet se puede hacer más implícita.
Existe un 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 de sonido . [8]
Otra forma de resolver las 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 fija a un punto al hacer contacto. El problema con esto es que es muy difícil elegir la fuerza impartida. Use una fuerza demasiado fuerte y los objetos se volverán inestables, demasiado débiles y los objetos se penetrarán entre sí. Otra forma es usar reacciones de colisión de proyección, que toma el punto ofensivo 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 compatible 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, se necesitaría controlar explícitamente las velocidades finales de los objetos que chocan (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 .
Ver también
- Condición de Courant-Friedrichs-Lewy
- Deriva de energía
- Integrador simpléctico
- Integración de Leapfrog
- El algoritmo de Beeman
Literatura
- ^ Verlet, Loup (1967). "Computadora" Experimentos "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 Bibliográfico : 1967PhRv..159 ... 98V . doi : 10.1103 / PhysRev.159.98 .
- ^ Presione, 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.
- ^ página web Archivada el 3 de agosto de 2004 en la Wayback Machine con una descripción del método Störmer.
- ^ Dummer, Jonathan. "Un método de integración de Verlet con corrección de tiempo simple" .
- ^ Swope, 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 .
- ^ Hairer, Ernst; Lubich, Christian; Wanner, Gerhard (2003). "Integración numérica geométrica ilustrada por el método Störmer / Verlet". Acta Numerica . 12 : 399–450. Bibcode : 2003AcNum..12..399H . CiteSeerX 10.1.1.7.7106 . doi : 10.1017 / S0962492902000144 .
- ^ Guía del usuario de SuperLU .
- ^ Baraff, D .; Witkin, A. (1998). "Grandes pasos en la simulación de tela" (PDF) . Procedimientos de gráficos por computadora . Serie de conferencias anuales: 43–54.
enlaces externos
- Demostración y código de integración de Verlet como un 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