Levy harmonic.py

From Werner KRAUTH

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

Iter = 1000000
Umean = 0.0
for iter in range(Iter):
    xstart = random.uniform(0.0, L)
    xend = xstart + L
    x = [xstart]
    for k in range(1, N):        # loop over internal slices
        x_mean = ((N - k) * x[k - 1] + xend) / (1.0 + N - k)
        sigma = math.sqrt(1.0 / (1.0  + 1.0 / (N - k) ))
        x.append(random.gauss(x_mean, sigma))
    x.append(xend)
    Umean += U(x)
print(Umean / Iter)

Further information

References

Krauth, W. XXX