ThermoIntegration.py

From Werner KRAUTH

import random
n = 100
NIter = 1000000
AlphaSim = [0, 0.25, 0.5, 0.75, 1.0]
Pi = {}
for Alpha in AlphaSim:
    Pi[(0, Alpha)] = 0.0
    Pi[(n + 1, Alpha)] = 0.0
    for x in range(1, n + 1):
        Pi[(x, Alpha)] =  (abs( (n + 1) / 2.0 - x)) ** Alpha
x = 1
ZProd = n
for AlphaIndex in range(len(AlphaSim) - 1):
    Alpha = AlphaSim[AlphaIndex]
    AlphaPrime = AlphaSim[AlphaIndex + 1]
    ZRatio = 0
    for iter in range(NIter):
        ZRatio += Pi[(x, AlphaPrime)] / Pi[(x, Alpha)]
        Sigma = random.choice([1, -1])
        if random.uniform(0.0, 1.0) < Pi[(x + Sigma, Alpha)] / Pi[(x, Alpha)]:
            x = x + Sigma
    ZRatio /= NIter
    print(ZRatio)
    ZProd *= ZRatio
print(1.0 / ZProd, 4.0 / n ** 2)