Canonic bosons.py

From Werner KRAUTH

Revision as of 00:36, 9 December 2016 Werner (Talk | contribs) ← Previous diff		Revision as of 23:06, 16 December 2018 Werner (Talk | contribs) Next diff →
Line 1:		Line 1:
	This is the python program Canonic_bosons.py useful for the homework session		This is the python program Canonic_bosons.py useful for the homework session
	of [[ICFP_Stat_Physics_2016|week 13 of my ICFP Lectures on statistical mechanics]].		of [[ICFP_Stat_Physics_2016|week 13 of my ICFP Lectures on statistical mechanics]].
-		+	Here, the density of state is 1,3,6,10,15, which corresponds to the three-dimensional isotropic harmonic trap. In the below program, we integrate in the complex plane from (-pi, lambda) to (pi, lambda) (neglecting the integral from (-pi, 0) to (-pi, lambda)), and always get the same result for the partition function. At T=1, for example, we get Z=17.373..., a result we can also obtain from the naive sum over the states, as implemented in Naive_bosons.py.
	import math, cmath, numpy, pylab		import math, cmath, numpy, pylab

---

Revision as of 23:06, 16 December 2018

This is the python program Canonic_bosons.py useful for the homework session of week 13 of my ICFP Lectures on statistical mechanics. Here, the density of state is 1,3,6,10,15, which corresponds to the three-dimensional isotropic harmonic trap. In the below program, we integrate in the complex plane from (-pi, lambda) to (pi, lambda) (neglecting the integral from (-pi, 0) to (-pi, lambda)), and always get the same result for the partition function. At T=1, for example, we get Z=17.373..., a result we can also obtain from the naive sum over the states, as implemented in Naive_bosons.py.

import math, cmath, numpy, pylab
 
N0vec = []
dos = [1, 3, 6, 10, 15]
T = 1.0
beta = 1.0 / T
Zint = complex(0.0, 0.0)
eps = 0.001
dellambda = 0.01
oldl = complex(-math.pi, dellambda)
complexi = complex(0.0, 1.0)
for RL in numpy.arange(-math.pi, math.pi, 0.00001):
    newl = RL + eps * complexi
    integrand = cmath.exp(-complexi * 5 * newl)
    for E in range(5):
        integrand /= (1.0 - cmath.exp( - beta * E + complexi * newl)) ** dos[E]
    Zint += integrand * (newl - oldl) / (2.0 * math.pi)
    oldl = newl
print Zint