Computerübung 2 - Sydney Stern - 10024300
import numpy as np
import math
import scipy.constants as sc
from scipy.optimize import root
import matplotlib.pyplot as plt
from scipy.optimize import root
from scipy.misc import derivative
from scipy import optimize
def expected_particle_number(n, μ, T):
nsq = n[0]*n[0] + n[1]*n[1] + n[2]*n[2]
return 1/(np.exp((nsq-μ)/(T))-1)
delta = 1e-8
def N_overall(μ, T, tol):
nmax = np.sqrt((np.log(1/tol+1)*T+μ)/3)
if np.isnan(nmax):
nmax = np.nan_to_num(nmax)
nmax = int(nmax.real) +1
nval = [expected_particle_number(np.array([i+1,j+1,k+1]),μ,T) for i in range (nmax) for j in range(nmax) for k in range(nmax)]
return np.sum(nval)
xval=np.linspace(-2,3.1,1000)
temp=[0.1,0.5,1,2,5]
yval=[ [N_overall(x,T,delta) for x in xval] for T in temp]
plt.style.use('fivethirtyeight')
fig = plt.figure(figsize=(15, 7))
ax = fig.add_subplot(1,1,1)
plt.title("Expected Particle Number $ N(\mu) $", fontsize=18)
plt.xlabel("Chemical Potential")
plt.ylabel("N")
plt.ylim(0,4)
plt.xlim(-2,3.1)
for i in range(len(temp)):
plt.plot(xval,yval[i],label=temp[i])
plt.grid(which='major')
plt.grid(which='minor', linestyle = ':')
plt.grid()
plt.legend()
plt.show()
/shared-libs/python3.7/py-core/lib/python3.7/site-packages/ipykernel_launcher.py:4: RuntimeWarning: invalid value encountered in sqrt
after removing the cwd from sys.path.
def μ(n,T):
tol=1e-8
guess=np.linspace(1,2.99,10)
guess=np.flipud(guess)
for x0 in guess:
rs =root(lambda μ: N_overall(μ,T,tol)-n, x0, method="lm")
if rs.success:
return rs.x[0]
print("not found")
def tbe(n):
return (4*n/(2.315*np.pi))**(2/3)
nvalues = np.arange(1,9,1)
tvalues = [np.linspace(0.01*tbe(n), 2*tbe(n) ,30) for n in nvalues]
μvalues = np.zeros((nvalues.size, tvalues[0].size))
for i in range (nvalues.size):
print(f"T_BE(n={nvalues[i]})=",tbe(nvalues[i]))
μvalues[i] = [μ(nvalues[i],T) for T in tvalues[i]]
plt.plot(tvalues[i]/tbe(nvalues[i]), μvalues[i], label=f"n={nvalues[i]}",marker=".")
plt.style.use('ggplot')
plt.title("Chem. Potential", fontsize=18)
plt.xlabel("$ T/T_{BE} $")
plt.ylabel("$\mu(N,T/T_{BE})$")
plt.grid(which='major')
plt.grid(which='minor', linestyle = ':')
plt.grid()
plt.legend()
plt.show()
T_BE(n=1)= 0.6712836696126203
/shared-libs/python3.7/py-core/lib/python3.7/site-packages/ipykernel_launcher.py:3: RuntimeWarning: overflow encountered in exp
This is separate from the ipykernel package so we can avoid doing imports until
T_BE(n=2)= 1.0655964033121468
/shared-libs/python3.7/py-core/lib/python3.7/site-packages/ipykernel_launcher.py:3: RuntimeWarning: overflow encountered in exp
This is separate from the ipykernel package so we can avoid doing imports until
T_BE(n=3)= 1.3963263018401306
T_BE(n=4)= 1.6915288515912315
T_BE(n=5)= 1.962845357319927
T_BE(n=6)= 2.216529840431889
T_BE(n=7)= 2.4564321652586356
T_BE(n=8)= 2.6851346784504813
lowest=[expected_particle_number([1,1,1],μvalues[i],tvalues[i]) for i in range(len(μvalues))]
secondlowest=[expected_particle_number([2,1,1],μvalues[i],tvalues[i]) for i in range(len(μvalues))]
fig = plt.figure(figsize=(15, 7))
plt.style.use('fivethirtyeight')
ax = fig.add_subplot(1,2,1)
plt.title("Lowest State $ N_{1,1,1} $", fontsize=18)
plt.xlabel("$T/T_{BE}$")
plt.ylabel("$ N_{1,1,1} $")
plt.ylim(0,max(nvalues)+0.1)
for i in range(nvalues.size):
plt.plot(tvalues[i]/tbe(nvalues[i]),lowest[i],label=f"n={nvalues[i]}")
plt.grid(which='major')
plt.grid(which='minor', linestyle = ':')
plt.grid()
plt.legend()
ax = fig.add_subplot(1,2,2)
plt.title("Lowest State $ N_{2,1,1} $", fontsize=18)
plt.xlabel("$T/T_{BE}$")
plt.ylabel("$ N_{1,1,1} $")
for i in range(nvalues.size):
plt.plot(tvalues[i]/tbe(nvalues[i]),secondlowest[i],label=f"n={nvalues[i]}")
plt.grid(which='major')
plt.grid(which='minor', linestyle = ':')
plt.grid()
plt.legend()
plt.show()
xvalues = np.linspace(1,100,100)
for T in [0.5, 1.5]:
plt.plot(
xvalues,
[expected_particle_number(n=[1,1,1], μ=μ(n=N,T=T*tbe(N)), T=T*tbe(N))/N for N in xvalues],
label=f"T = {T}"
)
#plt.plot(xvalues, [expected_particle_number()/N for N in xvalues], label=f"N={nvalues[i]})")
plt.vlines(x = 0.5, ymin = 0, ymax = 1)
#plt.vlines(x = 1.5, ymin = 0, ymax = 1)
plt.xlabel("$ T/T_{BE} $")
plt.ylabel("$ N_0/N $")
plt.title("$ N_{(0,0,0)}/N $")
plt.grid(which='major')
plt.grid(which='minor', linestyle = ':')
plt.xlim(1,100)
plt.grid()
plt.legend()
plt.show()
def energy(n):
return np.dot(n,n)
def system_energy(μ, T):
tolerance = 1e-5
N = 0
nmax = 2
N_new = 1
difference = 1
while (difference>tolerance):
n_values = [energy([i,j,k])*expected_particle_number([i,j,k], μ, T) for i in range(1,nmax+1) for j in range(1,nmax+1) for k in range(1,nmax+1)]
N_new = np.sum(n_values)
difference = N_new-N
nmax = nmax+1
N = N_new
return N
system_energy(-3,1)#test
def my_heat_capacity(T, N):
def f(T):
return system_energy(μ(N,T), T)
return derivative(func=f, x0=T, dx=1e-5)
n_values=np.arange(2,6,1)
t_values=np.linspace(0.1,10,50)
heat_capacities = []
for N in n_values:
TBE = tbe(N)
heat_capacities.append([my_heat_capacity(temp * TBE, N) for temp in t_values])
plt.style.use('fivethirtyeight')
plt.figure(figsize=(30, 14))
plt.xlabel('$T/T_{BE}$')
plt.ylabel('$C_N$')
plt.title('Heat Capacity $C_N$')
for i in range(n_values.size):
plt.plot(
t_values,
heat_capacities[i],
label = f"N ={n_values[i]}"
)
plt.legend()
plt.show()
specific_heat_capacities = []
for N in n_values:
TBE = tbe(N)
specific_heat_capacities.append([my_heat_capacity(temp*TBE, N)/N for temp in t_values]