montecarlo

This commit is contained in:
Vladimir Lemus 2023-10-05 15:04:53 -06:00
parent 844356c3c6
commit 44e182faf6
17 changed files with 512 additions and 10 deletions

21
#deca2.py# Normal file
View File

@ -0,0 +1,21 @@
from random import random
import matplotlib.pyplot as plt
import numpy as np
#Constantes
NTl = 1000 #numero inicial de thalios
tau = 3.053*60 #vida media del thalio en seg.
#los puntos a graficar
temps = []
#Determinando tiempo de decaimiento
for i in range(NTl):
t = -tau*np.log(1-random()/np.log(2))
temps.append(t)
#La grafica
(counts, bins, patches) = plt.hist(temps)
plt.xlabel("Tiempos")
plt.ylabel("Eventos")
plt.show()

1
.#deca2.py Symbolic link
View File

@ -0,0 +1 @@
devuan@bcm2711.5923:1696492413

15
axarinp.py Normal file
View File

@ -0,0 +1,15 @@
from random import random
N=100
agui=0
so=0
for i in range(N):
if random()<0.2:
print("Aguila")
agui+=1
else:
print("Sol")
so+=1
print("Aguilas: ", agui)
print("Soles: ", so)

7
axarinp.py~ Normal file
View File

@ -0,0 +1,7 @@
from random import random,
for i in range [0,10]:
if random()<0.2:
print("Aguila")
else:
print("Sol")

View File

@ -8,7 +8,7 @@ x=1
results=[]
for i in range(N):
x=(a*x+c)\%m
x=(a*x+c)%m
results.append(x)
plot(results,"o")

38
cascada.py Normal file
View File

@ -0,0 +1,38 @@
import numpy as np
import matplotlib.pyplot as plt
ZZ= 7 #nitrógeno
E0 = 20*10**9 #eV
EC = 600*10**6/ZZ #eV
def rand():return(np.random.uniform(0,1))
E=[E0]
X=[0]
N=[1]
T=[]
tmax = int(np.log(E0/EC)/np.log(2))
for t in range(0,tmax):
X.append(t+1)
#s=N[t]
for i in range(len(E)):
x=rand()
e=0.5*x*E[i]
e1=0.5*(1-x)*E[i]
if e>EC:
T.append(e)
if e1>EC:
T.append(e1)
prof=len(T)
N.append(prof)
E=T
T=[]
plt.plot(X,np.log(N))
plt.xlabel('Longitudes de radiación')
plt.ylabel('Número de partículas')
plt.title('Cascada electromagnética')
plt.grid(True)

33
cascada.py~ Normal file
View File

@ -0,0 +1,33 @@
import numpy as np
import matpltlib.pyplot as plt
ZZ= 7 #nitrógeno
E0 = 20*10**9 #eV
EC = 600*10**6/ZZ #eV
def rand():return(np.random.uniform(0,1))
E=[E0]
X=[0]
N=[1]
T=[]
tmax = int(np.log(E0/EC)/np.log(2))
for t in range(0,tmax):
X.append(t+1)
s=N[t]
for i in range(len(E)):
x=rand()
e=0.5*x*E[i]
e1=0.5*(1-x)*E[i]
if e>EC:
T.append(e)
if e1>EC:
T.append(e1)
prof=len(T)
N.append(prof)
E=T
T=[]

23
deca2.py Normal file
View File

@ -0,0 +1,23 @@
from random import random
import matplotlib.pyplot as plt
import numpy as np
#Constantes
NTl = 1000 #numero inicial de thalios
tau = 3.053*60 #vida media del thalio en seg.
#los puntos a graficar
temps = []
#Determinando tiempo de decaimiento
for i in range(NTl):
t = -tau*np.log(1-random()/np.log(2))
temps.append(t)
#Armamos el histograma
#La grafica
(counts, bins, patches) = plt.hist(temps)
plt.xlabel("Tiempos")
plt.ylabel("Eventos")
plt.show()

36
deca2.py~ Normal file
View File

@ -0,0 +1,36 @@
from random import random
from numpy import arange
from pylab import plot,xlabel,ylabel,show
#Constantes
NTl = 1000 #numero inicial de thalios
NPb = 0 #numero inicial de plomos
tau = 3.053*60 #vida media del thalio en seg.
th = 1.0 #tamano de paso en el tiempo seg.
p = 1-2**(-th/tau) #probabilidad de que decaiga un Tl
tmax = 1000 #tiempo total seg
#los puntos a graficar
tpuntos = arange(0.0,tmax,th)
Tlpuntos = []
Pbpuntos = []
#ciclo principal
for t in tpuntos:
Tlpuntos.append(NTl)
Pbpuntos.append(NPb)
#calculando el numero de nucleos que decayeron
decay=0
for i in range(NTl):
if random()< p:
decay+=1
NTl -= decay
NPb += decay
#La grafica
plot(tpuntos,Tlpuntos)
plot(tpuntos,Pbpuntos)
xlabel("Tiempo")
ylabel("Numero de nucleos")
show()

36
decaa.py Normal file
View File

@ -0,0 +1,36 @@
from random import random
from numpy import arange
from pylab import plot,xlabel,ylabel,show
#Constantes
NTl = 1000 #numero inicial de thalios
NPb = 0 #numero inicial de plomos
tau = 3.053*60 #vida media del thalio en seg.
th = 1.0 #tamano de paso en el tiempo seg.
p = 1-2**(-th/tau) #probabilidad de que decaiga un Tl
tmax = 1000 #tiempo total seg
#los puntos a graficar
tpuntos = arange(0.0,tmax,th)
Tlpuntos = []
Pbpuntos = []
#ciclo principal
for t in tpuntos:
Tlpuntos.append(NTl)
Pbpuntos.append(NPb)
#calculando el numero de nucleos que decayeron
decay=0
for i in range(NTl):
if random()< p:
decay+=1
NTl -= decay
NPb += decay
#La grafica
plot(tpuntos,Tlpuntos)
plot(tpuntos,Pbpuntos)
xlabel("Tiempo")
ylabel("Numero de nucleos")
show()

36
decaa.py~ Normal file
View File

@ -0,0 +1,36 @@
from random import random
from numpy import arange
from pylab import plot,xlabel,ylable,show
#Constantes
NTl = 1000 #numero inicial de thalios
NPb = 0 #numero inicial de plomos
tau = 3.053*60 #vida media del thalio en seg.
th = 1.0 #tamano de paso en el tiempo seg.
p = 1-2**(-th/tau) #probabilidad de que decaiga un Tl
tmax = 1000 #tiempo total seg
#los puntos a graficar
tpuntos = arange(0.0,tmax,h)
Tlpuntos = []
Pbpuntos = []
#ciclo principal
for t in tpuntos:
Tlpuntos.append(NTl)
Pbpuntos.append(NPb)
#calculando el numero de nucleos que decayeron
decay=0
for i in range(NTl):
if random()< p:
decay+=1
NTl -= decay
NPb += decay
#La grafica
plot(tpuntos,Tlpuntos)
plot(tpuntos,Pbpuntos)
xlabel("Tiempo")
ylabel("Numero de nucleos")
show()

BIN
decay2.jpg Normal file

Binary file not shown.

After

Width:  |  Height:  |  Size: 17 KiB

22
inte.py Normal file
View File

@ -0,0 +1,22 @@
import numpy as np
from random import random
#Definimos la funcion
def f(x):
return np.sin(1/(x*(2-x)))**2
#Definimos las inicializaciones
N= 1000
count = 0
#ciclo principal
for i in range(N):
x=2*random()
y= random()
if y<f(x):
count+=1
#Evaluamos la integral
I=2*count/N
print(I)

Binary file not shown.

View File

@ -912,11 +912,17 @@ El programa para hacerlo en python se puede ver en el \emph{listing} \ref{alg:az
\label{alg:azar}
\end{center}
El algoritmo es un generador de números aleatorios linealmente congruentes, y a pesar de que parecen ser aleatorios realmente están correlacionados. Para evitar las correlaciones se usa el algoritmo \textit{Mersenne twistter} que es un ``\textit{generalized feedback shift-register generator}'', algoritmo implementado en python.
El algoritmo es un generador de números aleatorios linealmente congruentes, y a pesar de que parecen ser aleatorios realmente están correlacionados. Para evitar las correlaciones se usa el algoritmo \textit{Mersenne twistter} que es un ``\textit{generalized feedback shift-register generator}'', algoritmo implementado en python. Implementado por Makoto Matsumoto y Takuji Nishimura en 1996, mejorado en 2002, la versión original es implementada en C y tiene las propiedades de:
Como se hizo con el ejemplo del generador linealmente congruente que se inició con un valor $x=1$ y se fueron obteniendo los demás valores de forma recursiva, lo mismo se hace en el algoritmo \textit{Mersenne twistter}, a esos valores se les llaman semillas.
\begin{itemize}
\item Perido: $2^{19937}-1$
\item Propiedad de equidistribución de $623$-dimensiones\footnote{Una secuencia de números reales es equidistribuida o uniformemente distribuida si la proporción de términos que caen en un subintervalo es proporcional a la longitud de tal subintervalo.}.
\item Uso eficiente de la memoria, consume sólo un espacio de trabajo de 624 palabras (la implementación en C).
\end{itemize}
Ahora hemos creado números todos con el mismo peso, pero si se quieren generar números de acuerdo a una distribución de probabilidad distinta, usar unos dados o una moneda cargada, por ejemplo, si tuviéramos una moneda cargada de tal forma que la probabilidad de que salga águila es el $20\%$ y de que caiga Sol $80\%$. Se corre el mismo algoritmo pero se ponen condiciones. Si el generador genera uniformemente números entre $0$ y $1$, se condiciona a que si el número está entre $0$ y $0.2$ elegimos que es águila, de lo contrario Sol, como se ve el el \textit{listing} \ref{alg:azar_pesado}.
Como se hizo con el ejemplo del generador linealmente congruente que se inició con un valor $x=1$ y se fueron obteniendo los demás valores de forma recursiva, lo mismo se hace en el algoritmo \textit{Mersenne twister}, a esos valores se les llaman semillas.
Ahora hemos creado números todos con el mismo peso, pero si se quieren generar números de acuerdo a una distribución de probabilidad distinta, por ejemplo, si tuviéramos una moneda cargada de tal forma que la probabilidad de que salga águila es el $20\%$ y de que caiga Sol $80\%$. Se corre el mismo algoritmo pero se ponen condiciones. Si el generador genera uniformemente números entre $0$ y $1$, se condiciona a que si el número está entre $0$ y $0.2$ elegimos que el resultado sea águila, de lo contrario Sol, como se ve el el \textit{listing} \ref{alg:azar_pesado}.
\begin{center}
\begin{lstlisting}[language=Python, caption=Generando números de acuerdo a una distribución.]
@ -950,7 +956,7 @@ Consideremos $1000$ núcleos de ${}^{208}Tl$, que tienen una vida media de $3.05
\begin{lstlisting}[language=Python, caption=Generando números de acuerdo a una distribución.]
from random import random
from numpy import arange
from pylab import plot,xlabel,ylable,show
from pylab import plot,xlabel,ylabel,show
#Constantes
NTl = 1000 #numero inicial de thalios
@ -961,7 +967,7 @@ Consideremos $1000$ núcleos de ${}^{208}Tl$, que tienen una vida media de $3.05
tmax = 1000 #tiempo total
#los puntos a graficar
tpuntos = arange(0.0,tmax,h)
tpuntos = arange(0.0,tmax,th)
Tlpuntos = []
Pbpuntos = []
@ -1018,7 +1024,7 @@ Conforme pasa el tiempo la probabilidad de decaimiento es menor. Esta distribuci
Para obtener a partir de las funciones aleatorias del lenguaje de programación, se hace uso del \textit{método de transformación}. Veamos en que consiste.
Supongamos que se tiene una distribución de probabilidad conocida, la función aleatoria del lenguaje de programación, que genera números $z$ con la probabilidad $q(z)$, con la probabilidad de generar un número entre $z$ y $z+dz$ dada por $q(z)dz$. Por otro lado se tiene una función $x=x(z)$, de tal forma que $x(z)$ también es un número aleatorio, pero puede tener una distribución de probabilidad distinta, llamémosle $p(x)$, distinta de $q(z)$.
Supongamos que se tiene una distribución de probabilidad conocida, la función aleatoria del lenguaje de programación, que genera números $z$ con la probabilidad $q(z)$, por lo que la probabilidad de generar un número entre $z$ y $z+dz$ está dada por $q(z)dz$. Por otro lado se tiene una función $x=x(z)$, de tal forma que $x(z)$ también es un número aleatorio, pero puede tener una distribución de probabilidad distinta, llamémosle $p(x)$, distinta de $q(z)$.
Por esta definición, la probabilidad de generar un número entre $z$ y $z+dz$ es igual a la probabilidad de generar un número entre $x$ y $x+dx$, es decir:
@ -1049,7 +1055,45 @@ De lo que resulta
t= -\tau ln\left( 1-\frac{z}{ln(2)}\right).
\end{equation}
Generamos números aleatorios uniformemente y evaluamos esa expresión.
Generamos números aleatorios uniformemente y evaluamos esa expresión.
El programa se puede ver en el listing \ref{alg:decay2}, que da como resultado el histograma en la figura \ref{fig:decay2}.
\begin{center}
\begin{lstlisting}[language=Python, caption=Tiempos de decaimiento en histograma.]
from random import random
import matplotlib.pyplot as plt
import numpy as np
#Constantes
NTl = 1000 #numero inicial de thalios
tau = 3.053*60 #vida media del thalio en seg.
#los puntos a graficar
temps = []
#Determinando tiempo de decaimiento
for i in range(NTl):
t = -tau*np.log(1-random()/np.log(2))
temps.append(t)
#La grafica
(counts, bins, patches) = plt.hist(temps)
plt.xlabel("Tiempos")
plt.ylabel("Eventos")
plt.show()
\end{lstlisting}
\label{alg:decay2}
\end{center}
\begin{figure}[ht!]
\begin{center}
\includegraphics[width=0.7\linewidth]{decay2.jpg}
\caption{Simulación del decaimiento del ${}^{208}Tl$, histograma de tiempos.}
\label{fig:decay2}
\end{center}
\end{figure}
\section{Distribución Gaussian}
@ -1110,7 +1154,7 @@ Podemos aproximar la integral sabiendo esto. Hay que generar número aleatorios
I \approx& \frac{\#\text{\{puntos que caen en I\}}}{\#\text{\{puntos en total generados\}}}\times A
\end{align*}
El programa en python que puede calcular esta integral por el método Montecarlo se ve en el \textit{listing} \ref{alg:int}
El programa en python que puede calcular esta integral por el método Montecarlo se ve en el \textit{listing} \ref{alg:inte}
\begin{center}
\begin{lstlisting}[language=Python, caption=Integración por Montecarlo.]
@ -1137,7 +1181,7 @@ El programa en python que puede calcular esta integral por el método Montecarlo
I=2*count/N
print(I)
\end{lstlisting}
\label{alg:decay}
\label{alg:inte}
\end{center}
\begin{thebibliography}{10}
@ -1146,6 +1190,7 @@ El programa en python que puede calcular esta integral por el método Montecarlo
\bibitem{PDG}Zyla, P.A. et al. (Particle Data Group), to be published in Prog. Theor. Exp. Phys. 2020, 083C01 (2020).
\bibitem{Leo}Leo, W.R. ``Techniques for nuclear and particle physics experiments'' Springer Verlag, 1987.
\bibitem{newman2013} Newman, M.E.J. ``Computational Physics''. \emph{Createspace Independent Pub}, 2013.
\bibitem{Marti-Shaw} Martin, B.R., Shaw, G. ``Particle Physics'' 3a edición, John Wiley \& Sons Ltd., 2008
\end{thebibliography}
\end{document}

BIN
pres5.pdf

Binary file not shown.

189
pres5.tex
View File

@ -908,6 +908,195 @@
\lstinputlisting[basicstyle=\ttfamily\scriptsize,firstline=1,lastline=16,language=python]{azarin.py}
\end{frame}
\begin{frame}{Algoritmo Mersenne twister}
\begin{itemize}
\item Makoto Matsumoto y Takuji Nishimura en 1996, mejora en 2002
\item Perido: $2^{19937}-1$
\item Propiedad de equidistribución de $623$-dimensiones\footnote{Una secuencia de números reales es equidistribuida o uniformemente distribuida si la proporción de términos que caen en un subintervalo es proporcional a la longitud de tal subintervalo.}.
\item Uso eficiente de la memoria, consume sólo un espacio de trabajo de 624 palabras (la implementación en C).
\end{itemize}
\end{frame}
\begin{frame}{Uso del algoritmo Mersenne twister}
\begin{itemize}
\item Un valor inicial: semilla ($x=1$ en el ejemplo anterior)
\item Opción con distinto peso: generar uniformemente en $[0,1]$ y dar pesos.
\end{itemize}
\lstinputlisting[basicstyle=\ttfamily\scriptsize,firstline=1,lastline=16,language=python]{axarinp.py}
\end{frame}
\begin{frame}{Decaimiento radiactivo}
Ley de deacímiento
\begin{equation*}
N(t)= N(0)2^{-t/\tau},
\end{equation*}
Probabilidad de decaimiento
\begin{equation*}
p(t)=1-2^{-t/\tau}
\end{equation*}
Consideramos $1000$ núcleos de ${}^{208}Tl$
\begin{itemize}
\item $tau_{{}^{208}Tl}=3.053$ minutos
\item Decae a ${}^{208}Pb$
\end{itemize}
\end{frame}
\begin{frame}{Simulemos}
\begin{itemize}
\item Generamos uniformemente
\item si cae entre 0 y $p(t)$ decae, de lo contrario no
\item Similar al águila o Sol contamos si es ${}^{208}Tl$ o si es ${}^{208}Pb$
\end{itemize}
\end{frame}
\begin{frame}{Programa I}
\lstinputlisting[basicstyle=\ttfamily\scriptsize,firstline=1,lastline=17,language=python]{decaa.py}
\end{frame}
\begin{frame}{Programa II}
\lstinputlisting[basicstyle=\ttfamily\scriptsize,firstline=17,lastline=37,language=python]{decaa.py}
\end{frame}
\begin{frame}{Resultado}
\begin{figure}[ht!]
\begin{center}
\includegraphics[width=0.7\linewidth]{decay.jpg}
\caption{Simulación del decaimiento del ${}^{208}Tl$ (azul) en ${}^{208}Pb$ (amarillo).}
\label{fig:decay}
\end{center}
\end{figure}
\end{frame}
\begin{frame}{Números aleatorios no uniformes}
La probabilidad de decaimiento diferencial
\begin{align*}
dp=& 1-2^{-dt/\tau} \\
=& 1-exp(ln(2^{-dt/\tau})) \\
=& 1-exp(\frac{-dt}{\tau}ln(2)) \\
=& 1-(1-\frac{dt}{\tau}ln(2)) \text{, aproximando a primer orden la exponencial}\\
=& \frac{dt}{\tau}ln(2)
\end{align*}
\end{frame}
\begin{frame}{Distribución de probabilidad}
La probabilidad de que un núcleo se mantenga sin decaer tras un tiempo $t$
\begin{equation}
P(t)=2^{-t/\tau}\frac{ln(2)}{\tau} dt
\end{equation}
Ya no es constante, conforme pasa $t$ cambia
\end{frame}
\begin{frame}{Método de transformación}
\begin{itemize}
\item Distribución de probabilidad genera $z$ con pobabilidad $q(x)$
\item Probabilidad de generar un número entre $z$ y $z+dz$: $q(z)dz$
\item Función $x=x(z)$, con $p(x)$
\end{itemize}
\begin{equation*}
p(x)dx=q(z)dz
\label{ec:igual}
\end{equation*}
\end{frame}
\begin{frame}{Determinar $x(z)$}
\begin{itemize}
\item Sea $q(z)$ la distribución uniforme:
\begin{equation}
q(z)=
\begin{cases}
1 & \text{si } z \in [0,1]\\
0 & \text{en cualquier otro caso }
\end{cases}
\end{equation}
\end{itemize}
\end{frame}
\begin{frame}{Transformando}
Integramos la igualdad entre las distribuciones
\begin{equation}
\int^{x(z)}_{-\infty} p(x')dx' = \int_0^z (1) dz' = z
\end{equation}
Por suerte el lado derecho es integrable
\begin{align*}
z=& \int^{t(z)}_{-\infty} p(t')dt'\\
=& \int^{t(z)}_{-\infty} \frac{ln(2)}{\tau}2^{-t'/\tau}dt'\\
=& ln(2)(1-e^{-t/\tau})
\end{align*}
\end{frame}
\begin{frame}{Resultado transformación}
\begin{equation}
t= -\tau ln\left( 1-\frac{z}{ln(2)}\right).
\end{equation}
\end{frame}
\begin{frame}{El programa}
\lstinputlisting[basicstyle=\ttfamily\scriptsize,firstline=1,lastline=22,language=python]{deca2.py}
\end{frame}
\begin{frame}{Histograma}
\begin{figure}[ht!]
\begin{center}
\includegraphics[width=0.7\linewidth]{decay2.jpg}
\caption{Simulación del decaimiento del ${}^{208}Tl$, histograma de tiempos.}
\label{fig:decay2}
\end{center}
\end{figure}
\end{frame}
\begin{frame}{Distribución Gaussiana}
\begin{equation*}
p(x)=\frac{1}{\sqrt{2\pi \sigma^2}}exp\left( -\frac{x^2}{2\sigma^2}\right)
\end{equation*}
¿La podemos integrar?
\begin{itemize}
\item Hacemos cambio de variable
\end{itemize}
\begin{align*}
p(x)dx \times p(y)dy =& \frac{1}{\sqrt{2\pi \sigma^2}}exp\left( -\frac{x^2}{2\sigma^2}\right)dx \times \frac{1}{\sqrt{2\pi \sigma^2}}exp\left( -\frac{y^2}{2\sigma^2}\right)dy \\
=& \frac{1}{\sqrt{2\pi \sigma^2}}exp\left( -\frac{x^2+y^2}{2\sigma^2}\right)dxdy,
\end{align*}
\end{frame}
\begin{frame}{En coordenadas polares}
\begin{align*}
p(r,\theta)drd\theta =& \frac{1}{2\pi\sigma^2}exp\left(-\frac{r^2}{2\sigma^2}\right)rdrd\theta \\
=& \frac{r}{\sigma^2}exp\left(-\frac{r^2}{2\sigma^2}\right)dr \times\frac{d\theta}{2\pi}
\end{align*}
Se obtienen los generadores
\begin{align*}
r=& \sqrt{-2\sigma^2 ln(1-z)} \\
\theta =& 2\pi z \text{, de ambos se obtiene}\\
x=& rcos\theta \\
y=& rsen\theta
\end{align*}
\end{frame}
\begin{frame}{Integración Montecarlo}
\begin{equation*}
I=\int_0^2 sin^2 \left[ \frac{1}{x(2-x)}\right]dx.
\label{ec:inte}
\end{equation*}
\begin{figure}[ht!]
\begin{center}
\includegraphics[width=0.7\linewidth]{inte.jpg}
\caption{Gráfica de la función a integrar en la ecuación \ref{ec:inte}}
\label{fig:int}
\end{center}
\end{figure}
\end{frame}
\begin{frame}{Esquema de la integración}
\begin{align*}
\frac{I}{A} \approx& \frac{\#\text{\{puntos que caen en I\}}}{\#\text{\{puntos en total generados\}}} \\
I \approx& \frac{\#\text{\{puntos que caen en I\}}}{\#\text{\{puntos en total generados\}}}\times A
\end{align*}
\end{frame}
\begin{frame}{El programa}
\lstinputlisting[basicstyle=\ttfamily\scriptsize,firstline=1,lastline=23,language=python]{inte.py}
\end{frame}
\begin{frame}{Paso de partículas a través de la materia}
\begin{itemize}
\item Los valores calculados por pedazos