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()
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()
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]