Metropolis harmonic.py

From Werner KRAUTH

Revision as of 18:14, 19 February 2025 Werner (Talk | contribs) (→Python program) ← Previous diff		Revision as of 18:15, 19 February 2025 Werner (Talk | contribs) (→Python program) Next diff →
Line 4:		Line 4:
	==Python program==		==Python program==
	import math, random		import math, random
-		+
	def U(x):		def U(x):
	U = 0.0		U = 0.0
Line 13:		Line 13:
	U += (x[k] - x_minus) ** 2 / 2.0		U += (x[k] - x_minus) ** 2 / 2.0
	return U		return U
-		+
	N = 8		N = 8
	L = 16		L = 16
	delta = 1.0		delta = 1.0
-		+
	x = [L / N * k for k in range(N)]		x = [L / N * k for k in range(N)]
	Iter = 1000000		Iter = 1000000

---

Revision as of 18:15, 19 February 2025

Contents 1 Context 2 Python program 3 Further information 4 References

Context

This page is part of my 2025 Oxford Lectures

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
    x.append(xend)

Further information

References

Krauth, W. XXX