Computer Exercise 3
import numpy as np
import matplotlib.pyplot as plt
import numpy.random as random
def create_probability_table(h, T):
def prob(spin,tot):
dE = 2 * spin * tot + 2 * h * spin
if dE > 0:
return np.exp( -dE / T )
else:
return 1
combos = [(-1,4), (-1,2), (-1,0), (-1,-2), (-1,-4), (1,4), (1,2), (1,0), (1,-2), (1,-4)]
return {(spin, tot): prob(spin, tot) for (spin,tot) in combos}
#create_probability_table(-8, 1.5)
# function that gives the positions of neighbors
def positions_of_neighbors(x, y, n):
oben = (x,(y+1) % n)
unten = (x, (y-1) % n)
rechts = ((x+1) % n, y)
links = ((x-1) % n, y)
return { 0: oben, 1: unten, 2: rechts, 3: links}
#positions_of_neighbors(2,2,5)
#Create an random grid for initial state
def initial_state(n):
state = 2*np.random.randint(2, size=(n,n))-1
return state
print(initial_state(3))
[[-1 -1 -1]
[-1 1 -1]
[ 1 -1 1]]
def metropolis_sweep(n, config, h, T):
probs = create_probability_table(h, T)
for process in range(n**2):
x = random.randint(0,n)
y = random.randint(0,n)
s = config[x, y]
neighbors = positions_of_neighbors(x, y, n)
spin_total = 0
for i in range(4):
spin_total += config[neighbors[i]]
if random.random() < probs[(s,spin_total)]:
config[x,y] *= -1
#test the metropolis sweep
grid = np.full(shape=(16,16), fill_value=+1)
metropolis_sweep(n=16, config=grid, h=-8, T=1.5)
grid
fig, ax = plt.subplots(figsize=(10, 10))
X, Y = np.meshgrid(range(17), range(17))
plt.pcolormesh(X, Y, grid, cmap=plt.cm.PiYG)
plt.colorbar()
plt.title('Test script of Metropolis Sweep for h = -8, T = 1.5')
plt.show()
/shared-libs/python3.7/py-core/lib/python3.7/site-packages/ipykernel_launcher.py:3: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
This is separate from the ipykernel package so we can avoid doing imports until
/shared-libs/python3.7/py-core/lib/python3.7/site-packages/ipykernel_launcher.py:4: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
after removing the cwd from sys.path.
# a function that calculates the magnetization, given a certain configuration
def calculate_magnetization(config, n):
return np.sum(config) / (n*n)
#test
calculate_magnetization(grid,16)
random_grid = initial_state(16)
magnetizations = [calculate_magnetization(random_grid, 16)]
while True:
metropolis_sweep(n=16, config=random_grid, h=0, T=1)
magnetizations.append(calculate_magnetization(random_grid, 16))
if np.abs(magnetizations[-1] - magnetizations[-2]) < 10**(-5):
break
#test
magnetizations
# plot the magnetisations
plt.style.use('seaborn')
fig = plt.figure(figsize=(15, 7))
#plt.xlim(0,)
plt.title("Evolution of Magnetization over the Simulation h= 0, T = 1, random initialised grid", fontsize=18)
plt.xlabel("Sweeps")
plt.ylabel("Magnetization")
plt.plot(magnetizations)
plt.show()
random_grid = initial_state(16)
magnetizations = [calculate_magnetization(random_grid, 16)]
while True:
metropolis_sweep(n=16, config=random_grid, h=0, T=5)
magnetizations.append(calculate_magnetization(random_grid, 16))
if np.abs(magnetizations[-1] - magnetizations[-2]) < 10**(-5):
break
# plot the magnetisations
plt.style.use('seaborn')
fig = plt.figure(figsize=(15, 7))
#plt.xlim(0,)
plt.title("Evolution of Magnetization over the Simulation, h= 0, T = 5, randomly initialised grid", fontsize=18)
plt.xlabel("Sweeps")
plt.ylabel("Magnetization")
plt.plot(magnetizations)
plt.show()
ordered_grid = np.full(shape=(16,16), fill_value=+1)
magnetizations = [calculate_magnetization(ordered_grid, 16)]
while True:
metropolis_sweep(n=16, config=random_grid, h=0, T=5)
magnetizations.append(calculate_magnetization(random_grid, 16))
if np.abs(magnetizations[-1] - magnetizations[-2]) < 10**(-5):
break
plt.style.use('seaborn')
fig = plt.figure(figsize=(15, 7))
#plt.xlim(0,)
plt.title("Evolution of Magnetization over the Simulation, h= 0, T = 5, ordered grid", fontsize=18)
plt.xlabel("Sweeps")
plt.ylabel("Magnetization")
plt.plot(magnetizations)
plt.show()
def calculate_equilibrium_time(n, h, T):
ordered_grid = np.full(shape=(n,n), fill_value=+1)
random_grid = initial_state(n)
# create random lattice
# create ordered lattice
magnetizations_ordered = np.empty(shape=10)
magnetizations_random = np.empty(shape=10)
i = 0
while True:
# sweep random lattice
metropolis_sweep(n, random_grid, h, T)
# sweep ordered lattice
metropolis_sweep(n, ordered_grid, h, T)
magnetizations_ordered[i%10] = calculate_magnetization(ordered_grid, n)
magnetizations_random[i%10] = calculate_magnetization(random_grid, n)
ordered_sum = np.sum(np.abs(magnetizations_ordered))
random_sum = np.sum(np.abs(magnetizations_random))
i += 1
if np.abs(random_sum - ordered_sum) < 0.5 and i >= 10:
break
time = i
return time
calculate_equilibrium_time(n=16, h=0, T=1)
def average_equilibrium_time(n, T, h, runs):
times = np.empty(shape=(runs), dtype=np.int32)
for i in range(runs):
times[i] = calculate_equilibrium_time(n, h, T)
return np.mean(times), np.std(times) / np.sqrt(runs)
temperatures = np.linspace(start=1, stop=4, num=61)
equilibrium_times = np.empty(shape=len(temperatures), dtype=np.int32)
standard_error_equilibrium_times = np.empty(shape=len(temperatures))
#create equilibrium times for given temperatures
for i, T in enumerate(temperatures):
equilibrium_times[i], standard_error_equilibrium_times[i] = average_equilibrium_time(n=16, T=T, h=0, runs=30)
#test
print(equilibrium_times)
[110 108 98 239 75 73 93 117 87 76 66 51 57 67 69 69 78 70
60 57 70 68 76 69 78 74 74 75 71 61 74 58 42 36 36 31
36 29 32 27 24 22 24 18 20 19 18 19 16 17 16 16 16 14
16 15 14 15 13 15 13]
plt.style.use('seaborn')
fig = plt.figure(figsize=(30, 14))
plt.title("Estimated Equilibrium Times", fontsize=25)
plt.xlabel("Temperature", fontsize = 20)
plt.ylabel("Equilibrium Time", fontsize = 20)
plt.scatter(temperatures, equilibrium_times)
plt.plot(temperatures, equilibrium_times)
plt.errorbar(temperatures, equilibrium_times, yerr=standard_error_equilibrium_times, fmt='o')
plt.show()
def estimate_magnetization(n, h, T, time):
magnetizations_ordered = np.empty(shape=100)
ordered_grid = np.full(shape=(n,n), fill_value=+1)
# create ordered lattice
for i in range(2*time):
# sweep ordered lattice
metropolis_sweep(n, ordered_grid, h, T)
for i in range(100):
metropolis_sweep(n, ordered_grid, h, T)
magnetizations_ordered[i] = np.abs(calculate_magnetization(ordered_grid, n))
return np.mean(magnetizations_ordered)
estimate_magnetization(n=16, h=0, T=3, time=equilibrium_times[0])
estimated_magnetizations = np.empty(shape=len(temperatures))
for i, T in enumerate(temperatures):
estimated_magnetizations[i] = estimate_magnetization(n=16, h=0, T=T, time=equilibrium_times[i])
print(estimated_magnetizations)
[0.99976562 0.99960938 0.99953125 0.99773437 0.99796875 0.99539062
0.99570313 0.99523438 0.99101562 0.98976562 0.985 0.98023437
0.97789063 0.97335938 0.96367187 0.96671875 0.9584375 0.95460937
0.92445312 0.91484375 0.91804688 0.90539063 0.89046875 0.87289063
0.82507812 0.84820312 0.19578125 0.765 0.31234375 0.48015625
0.56273438 0.47539063 0.3303125 0.34289063 0.24203125 0.194375
0.27375 0.19148438 0.20070313 0.275 0.15898438 0.23625
0.20921875 0.17632813 0.12304688 0.115625 0.1653125 0.11945312
0.136875 0.128125 0.14765625 0.10539063 0.13523437 0.09046875
0.10945312 0.10914063 0.09945313 0.1375 0.10117187 0.11851562
0.09164063]
plt.style.use('seaborn')
fig = plt.figure(figsize=(30, 14))
plt.title("Estimated Absolute Magnetization per Spin", fontsize=25)
plt.xlabel("Temperature", fontsize = 20)
plt.ylabel("Estimated Magnetization", fontsize = 20)
plt.scatter(temperatures, estimated_magnetizations)
plt.plot(temperatures, estimated_magnetizations)
plt.show()