Wegner 1d Exact.py

From Werner KRAUTH

Revision as of 23:18, 1 December 2019 Werner (Talk | contribs) ← Previous diff		Current revision Werner (Talk | contribs)
Line 1:		Line 1:
	import math, cmath, random		import math, cmath, random
	random.seed('Werner')		random.seed('Werner')
-		+
	N = 8 # take even		N = 8 # take even
	beta = math.sqrt(2.0)		beta = math.sqrt(2.0)
-		+
	PhiR = []		PhiR = []
	for k in range(N):		for k in range(N):
	PhiR.append(random.uniform(-1.0, 1.0)) # this is an example		PhiR.append(random.uniform(-1.0, 1.0)) # this is an example
-		+
	rvalues = range(N)		rvalues = range(N)
	kvalues = [2.0 * math.pi / N * float(l) for l in range(N)]		kvalues = [2.0 * math.pi / N * float(l) for l in range(N)]
	epskvalues = [4.0 * math.sin( k / 2.0) **2 for k in kvalues]		epskvalues = [4.0 * math.sin( k / 2.0) **2 for k in kvalues]
-		+
	Energy_real = -2.0 * N + (PhiR[N - 1] - PhiR[0]) ** 2		Energy_real = -2.0 * N + (PhiR[N - 1] - PhiR[0]) ** 2
	for r in range(N - 1):		for r in range(N - 1):
	Energy_real += (PhiR[r + 1] - PhiR[r]) ** 2		Energy_real += (PhiR[r + 1] - PhiR[r]) ** 2
-		+
	print Energy_real, 'Energy of the conf. in real representation'		print Energy_real, 'Energy of the conf. in real representation'
	#		#
Line 42:		Line 42:
	for k in range(1, N / 2):		for k in range(1, N / 2):
	Energy_fourier += 2.0 * epskvalues[k] * PhiK[k] * PhiK[N - k]		Energy_fourier += 2.0 * epskvalues[k] * PhiK[k] * PhiK[N - k]
-		+
	print Energy_fourier, 'energy of the conf. in Fourier representation'		print Energy_fourier, 'energy of the conf. in Fourier representation'
	#		#
Line 62:		Line 62:
	dummy += term * math.sqrt(2.0 / N)		dummy += term * math.sqrt(2.0 / N)
	print r, dummy.real, PhiR[0] - PhiR[r], ' r, Fourier term, direct term'		print r, dummy.real, PhiR[0] - PhiR[r], ' r, Fourier term, direct term'
-		+
	xdata = []		xdata = []
	ydata = []		ydata = []
Line 73:		Line 73:
	xdata.append(r)		xdata.append(r)
	print ydata		print ydata
-	~

---

Current revision

import math, cmath, random
random.seed('Werner')

N = 8 # take even
beta = math.sqrt(2.0)

PhiR = []
for k in range(N):
    PhiR.append(random.uniform(-1.0, 1.0)) # this is an example

rvalues = range(N)
kvalues = [2.0 * math.pi / N * float(l) for l in range(N)]
epskvalues = [4.0 * math.sin( k / 2.0) **2  for k in kvalues]

Energy_real = -2.0 * N + (PhiR[N - 1] - PhiR[0]) ** 2
for r in range(N - 1):
    Energy_real += (PhiR[r + 1] - PhiR[r]) ** 2

print Energy_real, 'Energy of the conf. in real representation'
#
# Fourier transform (by hand, no FFT)
#
PhiK = []
for fourierk in kvalues:
    dummy = 0.0
    for i in range(N):
        term = cmath.exp(- 1j * fourierk * rvalues[i]) * PhiR[i] / math.sqrt(N)
        dummy += term
    PhiK.append(dummy)
#
# real-valued Fourier coefficients, don't forget to count PsiK[N/2] = PhiK[N/2]
#
PsiKpos = [0.0]
PsiKneg = [0.0]
for k in range(1,N / 2):
    PsiKpos.append((PhiK[k] + PhiK[N - k]) / math.sqrt(2.0))
    PsiKneg.append((PhiK[k] - PhiK[N - k]) / (math.sqrt(2.0) * 1j))
#
# Express the energy in Fourier modes
#
Energy_fourier = - 2.0 * N +  4.0 * PhiK[N / 2] ** 2
for k in range(1, N / 2):
    Energy_fourier += 2.0 * epskvalues[k] * PhiK[k] * PhiK[N - k]

print Energy_fourier, 'energy of the conf. in Fourier representation'
#
# Energy in Fourier representation with real-valued terms
#
Energy_real = - 2.0 * N + 4.0 * PhiK[N / 2]**2
for k in range(1, N / 2):
    Energy_real += (PsiKneg[k] ** 2 + PsiKpos[k] ** 2) * epskvalues[k]
print Energy_real, ' energy of th conf. with real-valued Fourier coefficients'
#
# Correlation between site 0 and site k
#
print 'check correlation in real and (real Fourier)'
for r in range(N):
    dummy = PhiK[N / 2] * (1.0 - (-1) ** r) / math.sqrt(N)
    for k in range(1, N / 2):
        term = PsiKpos[k] * (1.0 - math.cos(kvalues[k] * r)) + \
               PsiKneg[k] * math.sin(kvalues[k] * r)
        dummy += term * math.sqrt(2.0 / N)
    print r, dummy.real, PhiR[0] - PhiR[r], ' r, Fourier term, direct term'

xdata = []
ydata = []
for r in range(0, N / 2 + 1):
    summe = (1.0 - (-1) ** r) ** 2 / (32.0)
    for k in range(1, N / 2):
        summe += math.sin(math.pi * k * r / N) ** 2 / (4.0 * math.sin(math.pi *
                  k / N) ** 2)
    ydata.append( math.exp(- 2.0 * summe / N / beta))
    xdata.append(r)
print ydata