Wegner 1d Exact.py

From Werner KRAUTH

Jump to: navigation, search
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
Personal tools