Maxwell-Boltzmann Distribution

We import the libraries we are going to need

In [1]:
import numpy as np
import matplotlib.pylab as plt
from scipy.optimize import curve_fit

We save the url of our data

In [2]:
url='https://raw.githubusercontent.com/jmsevillam/Herramientas-Computacionales-UniAndes/master/Data/Velocities.dat'

We load our data

In [3]:
data=np.genfromtxt(url)

We plot the data of $v_x$ and $v_y$

In [4]:
v_x=data[:,0]
v_y=data[:,1]
In [5]:
fig=plt.figure(figsize=(10,5))
ax1=fig.add_subplot(121)
ax1.hist(v_x,bins=50)
ax1.set_xlabel("$v_x$")
ax1.set_ylabel("Frequency")
ax2=fig.add_subplot(122)
ax2.hist(v_x,bins=50)
ax2.set_xlabel("$v_x$")
ax2.set_ylabel("Frequency")
Out[5]:
Text(0, 0.5, 'Frequency')

The total velocity $$ v=\sqrt{v_x^2+v_y^2} $$

In [6]:
vel=data[:,2]
In [7]:
histo=plt.hist(vel,bins=20)

We define our model as,

In [8]:
def f(x,a,s):
    return a*x*np.exp(-0.5*(x**2.)/s)

We separate our histogram to do the fit

In [9]:
x_bins=(histo[1][:-1]+histo[1][1:])/2.

We do the fit using the Initial conditions

In [10]:
popt,pcov=curve_fit(f,x_bins,histo[0],p0=[1000,1000**2.])

And finally plot the histogram and the fit function

In [11]:
fig=plt.figure(figsize=(5, 3))
ax=fig.add_subplot(111)
ax.set_xlabel(u'Velocidad '+r'$V_x$')
ax.set_ylabel(u'Frecuencia')
ax.hist(vel,label='Histograma',bins=20)
ax.plot(x_bins,f(x_bins,*popt), label='Ajuste')
plt.legend()
Out[11]:
<matplotlib.legend.Legend at 0x16b14359d68>

Wigner distribution

This problem can be sovlved easily using the function np.linalg.eigvals(), such that

In [12]:
m=1000
n=100
eigen_vals=[]
for i in range(n):
    matrix=np.random.normal(0,1,(m,m))
    matrix=matrix+matrix.T
    eigen_vals=np.append(eigen_vals,np.linalg.eigvals(matrix))

And the results are

In [13]:
fig=plt.figure(figsize=(7, 5))
ax=fig.add_subplot(111)

ax.hist(eigen_vals,bins=100);

Central Limit Theorem

We generate the random walks

In [14]:
data=[]
for i in range(10000):
    data.append(np.random.random(1000).mean())

We define our model

In [15]:
def f(x,a,b,c):
    return a*np.exp(-(x-b)**2./c)

To calculate the histogram,

In [16]:
histogram=plt.hist(data,bins=20)

We separate the bins

In [17]:
x_bins=(histogram[1][:-1]+histogram[1][1:])/2.
In [18]:
print(histogram)

(array([   2.,    2.,    6.,   29.,   55.,  179.,  375.,  690., 1046.,
       1438., 1635., 1534., 1255.,  841.,  498.,  256.,  104.,   41.,
         12.,    2.]), array([0.45932333, 0.46310116, 0.46687898, 0.47065681, 0.47443464,
       0.47821246, 0.48199029, 0.48576812, 0.48954594, 0.49332377,
       0.4971016 , 0.50087942, 0.50465725, 0.50843508, 0.5122129 ,
       0.51599073, 0.51976856, 0.52354638, 0.52732421, 0.53110204,
       0.53487986]), <a list of 20 Patch objects>)

we do the fit,

In [19]:
popt,pcov=curve_fit(f,x_bins,histogram[0])

And at the end, we plot the histogram and the fit function

In [20]:
x_plot=np.linspace(x_bins[0],x_bins[-1],1000)
plt.hist(data,bins=20)
plt.plot(x_plot,f(x_plot,*popt))
Out[20]:
[<matplotlib.lines.Line2D at 0x16b145b49e8>]