Tarea 9: Diagrama de fase IIM3019. Modelación de Sistemas Físicos

Abraham PradoA01213521

Monterrey, Nuevo León.

21 de febrero del 2012

Resumen

En el presente reporte se programa el método Monte Carlo metrópolis en Matlab, el objetivo de dicho programa esmostrar las configuraciones que minimizan la energía del sistema con diagramas de fases, al igual que observar como secomporta el sistema cuando se aplica el método Monte Carlo. También se realizan análisis de la evolución temporal devariables termodinámicas de interés tales como la temperatura y la energía total del sistema.

Palabras Clave: Método Monte Carlo - Diagrama de fase - Optimización estocástica

Índice

1. Teoría 21.1. Método Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2. Métodos para generar diagramas de fase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2.1. Histogramas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2.2. Diagramas de fases para ODE’s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2.3. Scatter Plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2.4. Quiver Plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2. Explicación del Código 32.1. Código 1: Dinámica molecular . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.2. Código 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

3. Análisis de resultados 43.1. Gráficas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

3.1.1. Diagramas de fase del primer código . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43.1.2. Diagramas de fase del segundo código . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53.1.3. Evolución de variables termodinámicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63.1.4. Configuración final . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

4. Conclusiones 9

1

F3019. Modelación de sistemas Físicos 2 ITESM Campus Monterrey

5. Bibliografía 9

1. Teoría

1.1. Método Monte Carlo

El método de Monte Carlo es una técnica numérica para calcular probabilidades y otras cantidades relacionadas,utilizando secuencias de números aleatorios.Además es una herramienta de investigación y planeamiento; básicamentees una técnica de muestreo artificial, empleada para operar numéricamente sistemas complejos que tengan componentesaleatorios. 1

Gracias a la constante evolución de las microcomputadoras, en lo que se refiere a su capacidad de procesamiento de lainformación, el método de Monte Carlo es cada ves más frecuentemente utilizado.

Esta metodología provee como resultado, incorporada a los modelos financieros, aproximaciones para las distribucionesde probabilidades de los parámetros que están siendo estudiados.

Para ello son realizadas diversas simulaciones donde, en cada una de ellas, son generados valores aleatorios para elconjunto de variables de entrada y parámetros del modelo que están sujetos a incertidumbre. Tales valores aleatoriosgenerados siguen distribuciones de probabilidades específicas que deben ser identificadas o estimadas previamente.

1.2. Métodos para generar diagramas de fase

Entre los métodos para generar diagramas de fase se mencionan los siguientes:

1.2.1. Histogramas

Dentro del análisis por Monte Carlo, por lo general se llega a tener n variables independientes y una salida escalar. Unamanera de visualizar los datos es un histograma de dos dimensiones. Desafortunadamente no hay soporte directo en Matlabpara hacer los histogramas, no obstante Laszlo Balkay desarrollo uno y esta disponible en Mathworks File Exchange. Elcomando hist2() recibe como entrada dos variables independientes y la escala de los parámetros. El resultado es una gráficaen dos dimensiones que muestra como estan distribuidos los puntos de acuerdo a esas dos variables.

1.2.2. Diagramas de fases para ODE’s

1.2.3. Scatter Plot

El comando scatter(X,Y,S,C) despliega círculos de color en las coordenadas específicadas por los vectores X y Y(quedeben ser del mismo tamaño). En particular este fue el comando que se utilizo para visualizar la evolución de los diagramasde fase para poder ver si el método Monte-Carlo converge a la mínima configuración.

1.2.4. Quiver Plot

Una manera de visualizar el retrato fase de un sistema dinámico de dos ecuaciones diferenciales autónomas es a travézde los siguientes comandos 2 :

[x1, x2] = meshgrid(-.5:0.05:0.5, -.5:.05:.5);x1dot = -x1 - 2 *x2 .*x1.^2+x2; %Note the use of .* and .^x2dot = -x1-x2;quiver(x1,x2,x1dot, x2dot)

Para visualizar a mayor detalle las ecuaciones resueltas y la gráfica correspondiente favor de checar la referencia [7].1http://www.cyta.com.ar/elearn/proyectoinversion/riesgo/riesgo.htm2Strawbridge , Marie. Using Matlab to draw phase portraits. University of Chicago

2

F3019. Modelación de sistemas Físicos 3 ITESM Campus Monterrey

2. Explicación del Código

2.1. Código 1: Dinámica molecular

El tamaño de paso escogido fue de e − 12 [s] . La selección del tamaño de paso es importante debido a que si no seescoje con las unidades correctas , no se podrá visualizar las energías bien , ni el método de Verlet funcionara. El númerode iteraciones fue de 10.

La función recibe de entrada las constantes del potencial, una matriz de 3 columnas y N renglones de las posiciones dela plata (plata), el tamaño de paso y el número de iteraciones (h).

Para aumentar la eficiencia del código se preasignaron las variables posición (r) , velocidad (v) , aceleración(ac) , energíacinética (T), energía potencial(U), energía total(E) en matrices de ceros. La implementación del método Verlet se realizocon los siguientes comandos:

r(:,:,i+1)=r(:,:,i)+tam_paso.*v(:,:,i)+tam_paso.^2 .*ac(:,:,i);ac(:,:,i+1)=-A01213521_derivada_sch(puntos,a,n,m,eps,c)./masa;v(:,:,i+1)=v(:,:,i)+0.5.*tam_paso.*(ac(:,:,i)+ac(:,:,i+1));

Donde la función A01213521_derivada_sch calcula la fuerza neta de las partículas de plata. Para el cálculo de las energíasse usaron los siguientes comandos:

U(i+1)=A01213521_suttonchen(puntos,eps,c,a,n,m);T(i+1)=masa./2.*sum(sum(v(:,:,i+1).^2));E(i+1)=U(i+1)+T(i+1);

Donde la función A01213521_suttonchen calcula la energía potencial Sutton-Chen con los parámetros anteriormenteexplicados.

Para implementar los diagramas de fases en cada instante se usaron los comandos:

figure(1)subplot(3,1,1)

scatter(V(:,1,i+1), U(:,1, i+1));...pause(0.005);

subplot(3,1,2)scatter(V(:,1,i+1), AC(:,1, i+1));...pause(0.005);

subplot(3,1,3)scatter(U(:,1,i+1), AC(:,1, i+1));...pause(0.005);

Donde V y AC son las velocidades y aceleraciones de cada partícula. Los diagramas de fase implementados permitenvisualizar la dinámica molecular de las partículas desde otra perspectiva, por lo que es muy útil su implementación en elcódigo para hacer análisis con más detalles.

2.2. Código 2:Aplicación de Monte Carlo

1. Con el fin de que el programa converja más rápido a la solución , las posiciones de las partículas se hicieron iguales aplata.mat. Posteriormente se fueron alterando de manera aleatoria sus posiciones. Las velocidades se inicializaron enel orden de 102 , debido a datos obtenidos del programa A0XXXXXXX_dinamica_sch anteriormente.

2. Las perturbaciones se definieron de la siguiente manera:

delta_v(:,:,i) = 1.5*10^(5)*(1- 2*rand([N,3]));delta_r(:,:,i) = 3e-19.*rand(N,3);

3

F3019. Modelación de sistemas Físicos 4 ITESM Campus Monterrey

3. Para llegar al cálculo de la energías , se uso el programa A0XXXXXXX_suttonchen_matriz que calcula la energía porpartícula de los puntos del sistema. Las energías cinéticas se calculan con:

V(:,1,i+1)= sqrt(sum(v(:,:,i+1).^2,2));

4. Si el resultado del muestreo es positivo, se aplica:

E(i+1) > E(i) r(:,:,i+1) = r(:,:,i); v(:,:,i+1) = v(:,:,i) ;

5. Al aplicar el método Monte Carlo , con el fin de obtener resultados negativos aceptados con frecuencia aceptable, semodifico ligeramente la probabilidad de transición como:

0.9 + 0.1*rand >exp(- df1(i)/T_monte )

Dependiendo de los valores de df1 (diferencia de energías), se fue modificando la temperatura de Monte Carlo.

6. Con el programa A0XXXXXXX_derivada_sch se calculo la fuerza que se ejerce en cada partícula por iteración, condicha fuerza se calculo la aceleración para poder observar su diagrama de fase.

7. Para observar los diagramas de fase se uso el comando scatter en las tres variables de interés(aceleración , velocidad yenergía potencial por partícula).

8. Adicionalmente se gráfico la temperatura y la energía total con el fin de poder analizar a detalle la termodinámica delproblema.

3. Análisis de resultados

3.1. Gráficas

3.1.1. Diagramas de fase del primer código

El código regresa en cada instante los diagramas de fase de la energía potencial, la velocidad y la aceleración , en uninstante dado se observan los siguientes diagramas:

Figura 1: Diagrama de fase de f(x)

Se observa que los puntos de mayor conglomeración son en U = −3,5 × 10−19 , AC = 4 × 1015. No obstante la velocidadno se estabiliza en ningún punto en específico.

4

F3019. Modelación de sistemas Físicos 5 ITESM Campus Monterrey

3.1.2. Diagramas de fase del segundo código

Para las siguientes perturbaciones y temperatura se obtienen los diagramas de fase siguientes:

T_monte= 3*10^(-15) ;delta_v(:,:,i) = 1.5*10^(5)*(1- 2*rand([N,3]));delta_r(:,:,i) = 3e-19.*rand(N,3);

Figura 2: Diagrama de fase de f(x)

Para las siguientes perturbaciones y temperatura se obtienen los diagramas de fase siguientes:

T_monte= 3*10^(-15) ;delta_v(:,:,i) = 1.5*10^(5)*(1- 2*rand([N,3]));delta_r(:,:,i) = 3e-19.*rand(N,3);0.6 + 0.1*rand >exp(- df1(i)/T_monte )

Figura 3: Diagrama de fase de f(x)

5

F3019. Modelación de sistemas Físicos 6 ITESM Campus Monterrey

Para las siguientes perturbaciones y temperatura se obtienen los diagramas de fase siguientes:

T_monte= 3*10^(-12) ;0.9 + 0.1*rand >exp(- df1(i)/T_monte )delta_v(:,:,i) = 1.5*10^(5)*(1- 2*rand([N,3]));delta_r(:,:,i) = 3e-19.*rand(N,3);

Figura 4: Diagrama de fase de f(x)

3.1.3. Evolución de variables termodinámicas

En la siguiente figura se muestra con evoluciono la energía potencial en el caso 1 en cada iteración con el fin de observarla naturaleza del potencial y obtener datos como puntos mínimos, máximos, etc:

Figura 5: Gráfica de energía potencial y cinética

6

F3019. Modelación de sistemas Físicos 7 ITESM Campus Monterrey

Se observa que el mínimo de la energía potencial se repite de manera periódica, el valor de U = −4,3 × 10−19J es elmínimo al que connverge. En la energía cinética no se observa que una estabilización , esto se puede ver de igual maneraen los diagrams de fase, las partículas se mueven entre sí dentro de un volumen nanométrico.

Para el caso 2 se modifico la perturbación y la temperatura de Monte Carlo como se explico anteriormente y se obtiene:

Figura 6: Gráfica de energía potencial y cinética

En el caso 3 se graficaron la temperatura y la energía cinética tambien, la temperatura crece de manera monótona , noobstante la energía cinética y la energía fluctuan en un rango proporcional a la perturbación dado al sistema:

Figura 7: Gráfica de variables termodinámicas

7

F3019. Modelación de sistemas Físicos 8 ITESM Campus Monterrey

La energía potencial se comporta igual que en el casos anteriores, no obstante la energía cinética es mucho mayor quela energía potencial y por lo tanto la energía total se observa igual que la cinética.

La simulación anterior se repitió para 3000 iteraciones, no obstante solo se observo una estabilización de la temperatura:

Figura 8: Gráfica de variables termodinámicas

3.1.4. Configuración final

Para las siguientes perturbaciones y temperatura se obtienen los diagramas de fase siguientes:

T_monte= 3*10^(-15) ;delta_v(:,:,i) = 1.5*10^(5)*(1- 2*rand([N,3]));delta_r(:,:,i) = 3e-19.*rand(N,3);

Figura 9: Configuración partículas plata para el caso I

8

F3019. Modelación de sistemas Físicos 9 ITESM Campus Monterrey

Cabe resaltar que el aunque el programa corre de manera indefinida, cuando las partículas alcanzan un punto deestabilidad se tardan en perder dicha configuración. Sin embargo hay partículas que salen del volumen confinado , por loque incrementan el tiempo de corrida del programa.

4. Conclusiones

En conclusión, se pudo simular las partículas y hallar configuraciones estables usando el método de Monte CarloMetrópolis, se gráfico la evolución de la temperatura para diferentes temperaturas y perturbaciones. Se observo la evoluciónde temperatura igualmente. Con los diagramas de fase se observo donde se obtenian las configuraciones más estables deenergía total. En simulaciones futuras se puede modificar el número de partículas para minimizar de manera global, deigual manera se puede agregar una condición que límite el rango de interacción entre las partículas, haciendo la búsquedade mínimos más compacta y rápida. Finalmente la búsqueda por Monte Carlo asegura siempre la existencia de al menos unmínimo , sin necesidad de conocer el potencial a fondo. Más análisis se requieren para determinar la naturaleza periódicadel potencial Sutton-Chen en las corridas.

5. Bibliografía

Referencias

[1] Guevara, E. : Optimización con Monte Carlo [Presentación de PowerPoint]. Recuperada de la base de datos del curso deModelación de Sistemas Físicos en Blackboard.

[2] R2011b MathWorks Documentation.

[3] Liu, Jun S. Monte Carlo strategies in scientific computing. New York : Springer, c2001

[4] Rubinstein, Reuven Y. Simulation and the Monte Carlo method , New York : Wiley, c1981

[5] James A. Glazier, Monte Carlo Methods and Statistical Physics, Mathematical Biology Lecture 4,6, the BiocomplexityInstitute

[6] Stéphane MOINS .Implementation of a simulated annealing algorithm for Matlab . Training performed in Electronicssystems . Linköping Institute of Technology

[7] Strawbridge , Marie. Using Matlab to draw phase portraits. University of Chicago

9