Metropolis harmonic.py

From Werner KRAUTH

1 Context
2 Python program
3 Further information
4 References

Context

This page is part of my public lectures on Algorithms and Computations in theoretical physics at the University of Oxford (see the announcement for a syllabus). For more information, see the main page of the public lectures. It presents the Metropolis Monte Carlo algorithm invented by Metropolis et al. (1953). For an in-depth discussion, see my below reference. Other programs in the same series are the Hamiltonian Monte Carlo algorithm and the Lévy construction, both for the harmonic model.

[edit]

Python program

import math, random

def U(x):
    U = 0.0
    for k in range(N):
        k_minus = (k - 1) % N
        x_minus = x[k_minus]
        if k == 0: x_minus -= L
        U += (x[k] - x_minus) ** 2 / 2.0
    return U

N = 8
L = 16
delta = 1.0
x = [L / N * k for k in range(N)]
Iter = 1000000
for iter in range(Iter):
    k = random.randint(0, N - 1)                   # random slice
    k_plus = (k + 1) % N
    k_minus = (k - 1) % N
    x_plus = x[k_plus]
    if k == N - 1: x_plus += L
    if k_plus == N: x_plus = x[0] + L
    x_minus = x[k_minus]
    if k == 0: x_minus -= L
    x_new = x[k] + random.uniform(-delta, delta)   # new position at slice k
    U_k = ((x_plus - x[k]) ** 2  + (x[k] - x_minus) ** 2) / 2.0
    U_kp = ((x_plus - x_new) ** 2  + (x_new - x_minus) ** 2) / 2.0
    if random.uniform(0.0, 1.0) < math.exp(- (U_kp - U_k)): x[k] = x_new

[edit]