Metropolis harmonic.py
From Werner KRAUTH
(Difference between revisions)
| Revision as of 18:13, 19 February 2025 Werner (Talk | contribs) ← Previous diff |
Revision as of 18:14, 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)] | ||
| Iter = 1000000 | Iter = 1000000 | ||
| for iter in range(Iter): | for iter in range(Iter): | ||
| - | xstart = random.uniform(0.0, L) | + | k = random.randint(0, N - 1) # random slice |
| - | xend = xstart + L | + | k_plus = (k + 1) % N |
| - | x = [xstart] | + | k_minus = (k - 1) % N |
| - | for k in range(1, N): # loop over internal slices | + | x_plus = x[k_plus] |
| - | x_mean = ((N - k) * x[k - 1] + xend) / (1.0 + N - k) | + | if k == N - 1: x_plus += L |
| - | sigma = math.sqrt(1.0 / (1.0 + 1.0 / (N - k) )) | + | if k_plus == N: x_plus = x[0] + L |
| - | x.append(random.gauss(x_mean, sigma)) | + | 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) | x.append(xend) | ||
| - | |||
| ==Further information== | ==Further information== | ||
Revision as of 18:14, 19 February 2025
Contents |
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
