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))
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()
# 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)
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)
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()