Metropolis harmonic.py

From Werner KRAUTH

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.

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

Further information

References

Krauth, W. Hamiltonian Monte Carlo vs. event-chain Monte Carlo: an appraisal of sampling strategies beyond the diffusive regime. ArXiv: 2411.11690