Factor ZigZag X2X4.py

From Werner KRAUTH

Revision as of 06:50, 11 June 2024 Werner (Talk | contribs) ← Previous diff		Revision as of 06:51, 11 June 2024 Werner (Talk | contribs) Next diff →
Line 4:		Line 4:
	def u(x):		def u(x):
	return x ** 2 / 2.0 + x ** 4 / 4.0		return x ** 2 / 2.0 + x ** 4 / 4.0
-		+
	x = 0.0		x = 0.0
	time_ev = 0.0		time_ev = 0.0
	sigma = 1 if x <= 0.0 else -1		sigma = 1 if x <= 0.0 else -1
-		+
	data = []		data = []
	n_samples = 10 ** 6		n_samples = 10 ** 6
Line 23:		Line 23:
	time_ev = new_time_ev		time_ev = new_time_ev
	sigma *= -1		sigma *= -1
-		+
	plt.title('Factor-Zig-Zag algorithm, anharmonic oscillator')		plt.title('Factor-Zig-Zag algorithm, anharmonic oscillator')
	plt.xlabel('$x$')		plt.xlabel('$x$')

---

Revision as of 06:51, 11 June 2024

import math
import random
import matplotlib.pyplot as plt
def u(x):
    return x ** 2 / 2.0 + x ** 4 / 4.0

x = 0.0
time_ev = 0.0
sigma = 1 if x <= 0.0 else -1

data = []
n_samples = 10 ** 6
while len(data) < n_samples:
    delta_u2 = -math.log(random.uniform(0.0, 1.0))
    delta_u4 = -math.log(random.uniform(0.0, 1.0))
    new_x2 = math.sqrt(2.0 * delta_u2)
    new_x4 = (4.0 * delta_u4) ** 0.25
    new_x = sigma * min(abs(new_x2), abs(new_x4))
    new_time_ev = time_ev + abs(new_x - x)
    for t in range(math.ceil(time_ev), math.floor(new_time_ev) + 1):
        data.append(x + sigma * (t - time_ev))
    x = new_x
    time_ev = new_time_ev
    sigma *= -1

plt.title('Factor-Zig-Zag algorithm, anharmonic oscillator')
plt.xlabel('$x$')
plt.ylabel('$\pi(x)$')
plt.hist(data, bins=100, density=True, label='data')
XValues = []
YValues = []
for i in range(-1000,1000):
    x = i / 400.0
    XValues.append(x)
    YValues.append(math.exp(- u(x)) / 1.93525)
plt.plot(XValues, YValues, label='theory')
plt.legend(loc='upper right')
plt.show()