Wegner 1d Exact.py
From Werner KRAUTH
(Difference between revisions)
Revision as of 23:18, 1 December 2019
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
~