Algoritmo De Verlet Para Movimiento Parabólico Y Semiparabólico

Algoritmo De Verlet Para Movimiento Parab´olico Y Semiparab´olico Laura Victoria Trujillo Taborda 2 de Octubre de 2016 Resumen En este informe se implementa el algoritmo de Verlet con el fin de dar soluci´on num´erica a problemas planteados en mec´anica cl´asica tales como el movimiento parab´olico y semiparab´ olico de una part´ıcula mediante la plataforma virtual Jupyter.

1.

Introducci´ on

1.1.

Algoritmo B´ asico De Verlet

Uno de los algoritmos usados con mayor frecuencia en el c´alculo de trayectorias de part´ıculas, especialmente en din´amica molecular, es el algoritmo de Verlet [1]. Este algoritmo se deduce a partir de la expansi´on en serie de Taylor de la funci´on de posici´on ~r(t) en ~r(t + ∆t) y ~r(t − ∆t), de la siguiente forma: 1 ~r(t + ∆t) = ~r(t) + ~v (t)∆t + ~a(t)∆t2 + 2

1~ b(t)∆t3 + O(∆t4 ) 6

(1)

1 1 ~r(t − ∆t) = ~r(t) − ~v (t)∆t + ~a(t)∆t2 − ~b(t)∆t3 + O(∆t4 ) (2) 2 6 Siendo ~v (t) la velocidad, ~a(t) la aceleraci´on y ~b(t) la tercera derivada de la posici´on ~r(t) con respecto al tiempo. Sumando las expresiones (1) y (2) se obtiene, ~r(t + ∆t) = 2~r(t) − ~r(t − ∆t) + ~a(t)∆t2 + O(∆t4 )

(3)

El error estimado que contiene la nueva posici´on ~r es del orden de ∆t4 , donde ∆t se define como el paso de tiempo en la simulaci´on. N´otese que para evaluar la nueva posici´on s´olo se necesita tener conocimiento de la posici´on anterior ~r(t − ∆t) y la posici´on actual ~r(t) [2]. Adem´as, una buena aproximaci´on para la velocidad ~v (t) (a segundo orden de precisi´on ∆t2 ) es obtenida de la resta de las ecuaciones (1) y (2),

1

~v (t) =

~r(t + ∆t) − ~r(t − ∆t) + o(∆t2 ) 2∆t

(4)

Las expresiones (3) y (4) conforman el algoritmo b´asico de Verlet o Leap-Frog [4].

1.2.

Aplicaci´ on Algoritmo De Verlet Para Simulaci´ on Movimiento Parab´ olico

Una aplicaci´on interesante para el algoritmo de Verlet es el c´alculo de trayectorias de una part´ıcula que interacciona con un campo gravitacional constante, como es el caso de un movimiento parab´olico o movimiento de un proyectil. Para tal movimiento se admitir´a, en concordancia con [3],que no existen resistencias del medio, por lo que la u ´nica fuerza actuante sobre la part´ıcula es la gravedad, que suponemos ˆ Las ecuaciones son: definida por el campo gravitatorio simplificado (−mg k). m¨ x = 0; m¨ y = 0; m¨ z = −mg Si tomamos unos ejes en los que el plano vertical Oxz contenga a la velocidad inicial ~v0 , es f´acil comprobar que el movimiento se desarrollar´a dentro del mismo plano vertical: y˙ = y¨ = 0 ⇒ y = 0

Figura 1: Trayectoria parab´olica de una part´ıcula en el vac´ıo (Goicolea,2010, p.2.18). Denominando ϕ0 el a´ngulo de la velocidad en el lanzamiento con la horizontal (Figura 1) las ecuaciones en x y z se integran de manera trivial: x˙ = v0 cos(ϕ0 ); z˙ = −gt + v0 sin(ϕ0 ) 1 x = v0 cos(ϕ0 )t; z = − gt2 + v0 sin(ϕ0 )t 2

2

(5)

Las expresiones obtenidas en (5) son las llamadas ecuaciones horarias de la trayectoria, es decir, las ecuaciones que permiten obtener la posici´on en funci´on del tiempo. La trayectoria descrita por la part´ıcula queda definida por las ecuaciones horarias de forma param´etrica, mediante el param´etro t. Despejando t de la ecuaci´on de x y sustituyendo en la expresi´on de z, se obtiene la ecuaci´on impl´ıcita de la trayectoria: z=−

gx2 1 + x tan(ϕ0 ) 2 (v0 cos(ϕ0 ))2

(6)

ecuaci´on que representa una par´abola de eje vertical y por cuya raz´on se conoce a este movimiento como parab´olico. Con lo mencionado anteriormente, el esquema del algoritmo de Verlet ser´ıa el siguiente: ~x(t + ∆t) = 2x~0 − ~x(t − ∆t) + ~a(t)∆t2

(7)

con condiciones iniciales x~0 = (0, 0, 0) y v~0 = (v0x , 0, v0z )

1.3.

Aplicaci´ on Algoritmo De Verlet Para Simulaci´ on Movimiento Semiparab´ olico

El movimiento semiparab´olico de una part´ıcula es un caso especial del movimiento parab´olico. Es decir, las ecuaciones que describen el movimiento son las mismas, sin embargo las condiciones iniciales son diferentes. Para este caso, se considerar´a una rampa encargada de acelerar la part´ıcula desde una altura h (Figura 2) donde, p (8) v0 = 2 2ghr cuya expresi´on es obtenida a partir del principio de conservaci´on de energ´ıa.

Figura 2: Montaje experimental movimiento semiparab´olico.

3

Con ello, el esquema del algoritmo de Verlet ser´ıa, ~x(t + ∆t) = 2x~0 − ~x(t − ∆t) + ~a(t)∆t2

(9)

con condiciones iniciales θ = 0◦ , x~0 = (0, 0, h) y v~0 = (v0 , 0, 0)

2.

Metodolog´ıa

En la realizaci´on de este informe se consideraron dos experimentos tratados en mec´anica cl´asica conocidos como movimiento parab´olico y semiparab´olico de una part´ıcula, donde se infiri´o una aceleraci´on gravitatoria constante cerca de la superficie de Marte de magnitud 3,71 sm2 [7]. Para llevar a cabo ambos experimentos, se implement´o el algoritmo de Verlet bajo condiciones iniciales dadas para cada uno de ´estos; las cuales ser´an detalladas en las siguientes secciones.

2.1.

Experimento No.1: Movimiento Parab´ olico

Para este experimento se consider´o una part´ıcula lanzada inicialmente a 1000 ms con un a´ngulo de tiro de 51◦ sobre la superficie marciana. Es decir, que las condiciones iniciales est´an dadas como, x~0 = (0, 0, 0)m v~0 = (1000 cos 51◦ , 0, 1000 sin 51◦ )

m s

(10)

Para su implementaci´on se calcul´o en primera instancia la posici´on anterior ~x(t − ∆t) transcurrido cierto ∆t, para este caso se eligi´o ∆t = 0,01s, mediante la expresi´on (2). Obtenida tal posici´on y dadas las condiciones iniciales (10), se determina la posici´on posterior ~x(t + ∆t) a partir de (7). Adem´as, se tendr´a en cuenta la posibilidad de que exista una barrera o pared que limite la trayectoria horizontal de la part´ıcula para as´ı observar la altura h en la que ´esta part´ıcula colisiona con tal barrera. Para ello, se escogieron seis distancias en las cuales se ubic´o dicha pared que fueron respectivamente: 50000m, 250000m, 8000m, 180000m, 500m y 300000m. Dicha implementaci´on fue realizada en el lenguaje de programaci´on Python 3.4, utilizando librerias como matplotlib, numpy y scipy.

2.2.

Experimento No.2: Movimiento semiparab´ olico

Para este experimento se consider´o una part´ıcula acelerada mediante una rampa cuya altura es 0,66m, donde a partir de (8), es lanzada aproximadamente a 2,21 ms desde una altura de 7,10m sobre la superficie marciana. Las condiciones iniciales para tal experimento son, x~0 = (0, 0, 7,10)m v~0 = (2,21, 0, 0) 4

m s

(11)

Para su implementaci´on se calcul´o en primera instancia la posici´on anterior ~x(t − ∆t) transcurrido cierto ∆t, para este caso se eligi´o ∆t = 0,001s, mediante la expresi´on (2). Obtenida tal posici´on y dadas las condiciones iniciales (11), se determina la posici´on posterior ~x(t + ∆t) a partir de (9). Adem´as, se tuvo en cuenta la posibilidad de que existiera una barrera o pared que limitara la trayectoria horizontal de la part´ıcula para as´ı observar la altura h en la que ´esta part´ıcula colisiona con tal barrera. Para ello, se escogieron seis distancias en las cuales se ubic´o dicha pared que fueron respectivamente: 4,0m, 1,0m, 3,0m, 0,7m, 2,5m y 0,4m. Dicha implementaci´on fue realizada en el lenguaje de programaci´on Python 3.4, utilizando librerias como matplotlib, numpy y scipy.

3.

Resultados

A continuaci´on se muestran los resultados obtenidos para el experimento no. 1 movimiento parab´olico y experimento no. 2 movimiento semiparab´olico.

Experimento No. 1

Figura 3: Trayectoria part´ıcula sin ruido. En la Figura 3 se obtiene la ecuaci´on de movimiento de la part´ıcula como z = −1,86t2 + 777,14t, para la cual los valores de gravedad y velocidad se obtienen a partir de (5) como,

5

1 m − g = −1,86 ⇒ g = −3,72 2 2 s m s lo cual coincide con los valores te´oricos con un porcentaje de error menor al 1 %. Esto era de esperarse ya que en el c´alculo de la trayectoria no se agreg´o ruido alguno que generar´a un error mayor en el resultado esperado. A su vez, se obtiene el punto de impacto en el eje horizontal a 263651,2m cuyo error absoluto es 0,4 con el obtenido te´oricamente 263651,6m y un tiempo de ca´ıda de 418,9s con un error del 0,003. V0y = 777,14

Figura 4: Colisi´on barrera-part´ıcula con un ruido de 0.1 Para cada una de las distancias de barrera propuestas se obtiene el punto de impacto de la part´ıcula con ´esta, como se muestra en Tabla 1. Como se aprecia en Tabla 1, los puntos de impacto concuerdan con los obtenidos en la Figura 4. As mismo, se observa que la precisi´on del impacto es del orden de d´ecimas y su 6

Cuadro 1: Punto de impacto para determinada barrera x (m) Impacto Z (m) 50000 50034.92 250000 15986.15 8000 9579.320 180000 70525.92 500 615.8844 300000 0.520017

distribuci´on no tiene la exactitud deseada debido al ruido agregado, pues entre menor ruido o perturbaci´on se a˜ nada al sistema, mayor exactitud y estabilidad se tendr´a.

Figura 5: Colisi´on barrera-part´ıcula con un ruido de 0.01 Para la figura 5 se ve claramente como la disminuci´on de ruido mejora considerablemente la precisi´on del impacto(ver Tabla 1). En la 5.1 se ve la existencia de cierta bimodalidad, la 7

cual se presume es debida al valor de ∆t, ya que si tal ∆t se hiciese mas pequeo se observar´ıa una monomodalidad como en las dem´as figuras. Sin embargo, para la figura 5.6 se ve una dispersi´on por fuera del rango esperado (entre 0 y 2) debido a la perturbaci´on dada.

Figura 6: Colisi´on barrera-part´ıcula con un ruido de 0.001 Al igual que en las anteriores gr´aficas, se verifica la hipot´esis de a menor ruido, mayor precisi´on. Por ejemplo, en la figura 6.1 vemos que el punto de impacto de la part´ıcula oscila entre 50034,75m y 50035,10m valores que est´an, en comparaci´on con la figura 4 (50000m y

8

50070m), m´as cercanos con el valor de impacto te´orico.

Figura 7a: Histograma 2D altura Z en funci´on de la distancia Y para 10.000 repeticiones del experimento. 9

Para un ruido de 0.001 se obtiene para determinada barrera el punto de impacto m´as frecuente de la part´ıcula. Para una distancia de barrera de 50000m (Figura 7a.1)se obtiene que el valor m´as frecuente est´a entre 50034,9m y 50035m en la posici´on 0, 0m, resultado que es satisfactorio con el punto de impacto esperado 50034,92m. Para una distancia de barrera de 250000m (Figura 7a.2)se obtiene que el valor m´as frecuente es aproximadamente 15988,2m en la posici´on 0, 0m de ancho de barrera, resultado que es satisfactorio con el punto de impacto esperado 15986,15m. Para una distancia de barrera de 8000m (Figura 7a.3)se obtiene que el valor m´as frecuente es aproximadamente 9577,86m en posici´on 0, 0m,resultado que es satisfactorio con Tabla 1. Para una distancia de barrera de 180000m (Figura 7a.3)se obtiene que el valor m´as frecuente es aproximadamente 70526,5m en posici´on 0, 0m,resultado que es satisfactorio con Tabla 1. Para una distancia de barrera de 500m (Figura 7a.3)se obtiene que el valor m´as frecuente es aproximadamente 612,78m en posici´on 0, 0m,resultado que es satisfactorio con Tabla 1. Para una distancia de barrera de 300000m (Figura 7a.3)se obtiene que el valor m´as frecuente es aproximadamente 5,1m en posici´on 0, 0m,resultado que es satisfactorio con Tabla 1. Para un ruido de 0.001, se obtuvieron los puntos de impacto m´as frecuente para la altura Z donde para la figura 7b.1 tiende a estar a una altura de 50034,9m en concordancia con la figura 7a.Con esto, se afirma que el error absoluto obtenido para una barrera ubicada a 50000m del eje X es de 0,02. Mientras que los puntos de impacto m´as frecuente para la altura Z en la figura 7b.2 tiende a estar a una altura de 15988,3m en concordancia con la figura 7a.Con esto, se afirma que el error absoluto obtenido para una barrera ubicada a 250000m del eje X es de 2,2. Los puntos de impacto m´as frecuente para la altura Z en la figura 7b.3 tiende a estar a una altura de 9577, 87m en concordancia con la figura 7a.Con esto, se afirma que el error absoluto obtenido para una barrera ubicada a 8000m del eje X es de 1,5. Para la figura 7b.4, los puntos de impacto m´as frecuente para la altura Z en la tiende a estar a una altura de 70526, 5m en concordancia con la figura 7a.Con esto, se afirma que el error absoluto obtenido para una barrera ubicada a 180000m del eje X es de 0,6. Mientras que los puntos de impacto m´as frecuente para la altura Z en la figura 7b.5 tiende a estar a una altura de 612,78m en concordancia con la figura 7a. Con esto, se afirma que el error absoluto obtenido para una barrera ubicada a 500m del eje X es de 3,1. Los puntos de impacto m´as frecuente para la altura Z en la figura 7b.6 tiende a estar a una altura de 5,1m en concordancia con la figura 7a.Con esto, se afirma que el error absoluto obtenido para una barrera ubicada a 300000m del eje X es de 4,6.

10

Figura 7b: Histograma altura Z.

11

Figura 7c: Histograma distancia Y (ancho barrera). En la figura 7c se ve claramente para cada una de las distancia de barrera que existe mayor probabilidad de que la part´ıcula impacto en el punto 0,0 de la barrera. Es decir, hay mayor tendencia a que la part´ıcula choque exactamente en la mitad del ancho de ´esta.

12

Figura 7d: Histograma Tiempo Ca´ıda. Como se observa en la figura 7d, se obtiene para cada una de las gr´aficas una distribuci´on uniforme lo que indica que los datos obtenidos para el tiempo de ca´ıda es consistente, es decir, la frecuencia es la misma para cada dato. En pocas palabras, esta distribuci´on uniforme nos dice que los datos obtenidos no se pueden dividir en intervalos suficientes.

13

Experimento No. 2 Movimiento Semiparab´ olico

Figura 8: Trayectoria part´ıcula sin ruido. En la Figura 8 se obtiene la ecuaci´on de movimiento de la part´ıcula como z = −1,86t2 + 7,10, para la cual los valores de gravedad y velocidad se obtienen a partir de (5) como, 1 m − g = −1,86 ⇒ g = −3,72 2 2 s V0y = 0

m s

m s lo cual coincide con los valores te´oricos con un porcentaje de error menor al 1 %. Esto era de esperarse ya que en el c´alculo de la trayectoria no se agreg´o ruido alguno que generar´a un error mayor en el resultado esperado. A su vez, se obtiene el punto de impacto en el eje horizontal a 4,3227m cuyo error absoluto es 0,0057 con el obtenido te´oricamente 4,3284m y un tiempo de ca´ıda de 1,9563s con un error del 0,0003. h0 = 7,10

14

Figura 9: Colisi´on barrera-part´ıcula con un ruido de 0.1 En la siguiente tabla se registran los puntos de impacto para cada una de las barreras elegidas. Cuadro 2: Punto de impacto para determinada barrera x (m) Impacto Z (m) 4.0 1.029547 1.0 6.721016 3.0 3.684112 0.7 6.914767 2.5 4.727156 0.4 7.039898

Como se ve en la figura 9 los datos obtenidos oscilan entre el valor esperado en Tabla 2.

15

Figura 10: Colisi´on barrera-part´ıcula con un ruido de 0.01 Como se mencion´o anteriormente, se aprecia que para figura 10 la precisi´on del impacto es mayor debido a que el ruido disminuye de 0,1 a 0,01. Para estas gr´aficas se repiti´o el experimento 100 veces, motivo por el cual no se observa una distribuci´on uniforme deseada. En contraste con la figura 12a que se necesito como m´ınimo de repetir este experimento unas 10000 para obtener el dato mas frecuente de impacto de manera clara.

16

Figura 11: Colisi´on barrera-part´ıcula con un ruido de 0.001 Al igual que en las gr´aficas anteriores,se advierte que el punto de impacto es cada vez m´as exacto al esperado en Tabla 2. Por ejemplo, en la figura 11.1 vemos que el punto de impacto de la part´ıcula oscila entre 1,055m y 1,030m lo que lo hace m´as exacto en contraste con la figura 9(2,0m y 0,0m). Tambi´en se observa cierta bimodalidad en la figura 11.1 y figura 11.2, para lo cual ser´ıa necesario aumentar el ∆t para garantizar una distribuci´on uniforme unimodal.

17

Figura 12a: Histograma altura Z en funci´on de distancia Y (ancho barrera).

18

Figura 12b: Histograma altura Z. Para un ruido de 0.001, se obtuvieron los puntos de impacto m´as frecuente para la altura Z donde para la figura 12b.1 tiende a estar a una altura de 1,047m en concordancia con la figura 12a.Con esto, se afirma que el error absoluto obtenido para una barrera ubicada a 4,0m del eje X es de 0,017. Mientras que los puntos de impacto m´as frecuente para la altura Z en la figura 12b.2 tiende a estar a una altura de 6,723m en concordancia con la figura 12a.Con esto, se afirma que el error absoluto obtenido para una barrera ubicada a 1,0m del eje X es de 0,002. Los puntos de impacto m´as frecuente para la altura Z en la figura 12b.3 tiende a estar a una altura de 3,695m en concordancia con la figura 12a.Con esto, se afirma que el error 19

absoluto obtenido para una barrera ubicada a 3,0m del eje X es de 0,01. Para la figura 12b.4, los puntos de impacto m´as frecuente para la altura Z en la tiende a estar a una altura de 6,9145m en concordancia con la figura 12a.Con esto, se afirma que el error absoluto obtenido para una barrera ubicada a 0,7m del eje X es de 0,002.

Figura 12c: Histograma distancia Y (ancho barrera). Mientras que los puntos de impacto m´as frecuente para la altura Z en la figura 12b.5 tiende a estar a una altura de 4,735m en concordancia con la figura 12a.Con esto, se afirma que el error absoluto obtenido para una barrera ubicada a 2,5m del eje X es de 0,0012. Los puntos de impacto m´as frecuente para la altura Z en la figura 12b.6 tiende a estar a una 20

altura de 7,048m en concordancia con la figura 12a.Con esto, se afirma que el error absoluto obtenido para una barrera ubicada a 0,4m del eje X es de 0,009.

Figura 12d: Histograma tiempo ca´ıda. Cabe a˜ nadir que los gr´aficos de puntos bimodales como los observados en figura 12a.1,2,3,4 son reflejo del ∆t escogido, pues si el paso de tiempo se hace menor, tal distribuci´on tendr´a una uniformidad mayor y con ello pasar´an a una monomodalidad, obteniendo entonces una tendencia mayor en cierta regi´on del espacio al momento del impacto. A pesar de que se observa claramente este comportamiento en tales gr´aficas de histogramas bidimensionales, fue necesario realizar 10000 veces el experimento para su correcta visualizaci´on. 21

Referencias [1] L. Verlet, Phys Rev 159,98,Phys Rev 165,201 (1967).

[2] Algoritmo Verlet Repositorio.pedagogica.edu.co. Retrieved 25 September 2016, from http://repositorio.pedagogica.edu.co/xmlui/bitstream/handle/123456789/1353/ANEXO20A. [3] J. Goicolea Curso De Mec´anica. (Colecci´on Escuelas,2.10). [4] D. Frenkel B. Smith, Understanding Molecular Simulation From Algorithms To Applications. (Academic Press,1996). [5] Resnick,R. and Halliday,D.(1960) PHYSICS.-For Students of Science and Engineering (Part I). (New York, John Wiley & SONS inc), pp.83,89. [6] Rojas,J. and Martinez,R.(2014) Mec´anica 3d: Python y el algoritmo de Verlet. Revista Mexicana de F´ısica http://www.scielo.org.mx/pdf/rmfe/v60n1/v60n1a5.pdf

[7] Williams, D.(2016) Mars Fact Sheet. Nasa. http://nssdc.gsfc.nasa.gov/planetary/factsheet/ma 3 Oct. 2016].

22