<p>我正在寻找一个1D伊辛模型的简单实现,偶然发现了这篇文章。虽然我不是这个领域的专家,但我确实写了一篇相关主题的硕士论文。我在Oriol Cabanas Tirapu的答案中实现了代码,并发现了一些bug(我想)</p>
<p>下面是我的改编版哦他们的代码。希望它对某些人有用</p>
<pre><code>#Coding attempt MCMC 1-Dimensional Ising Model
import numpy as np
import matplotlib.pyplot as plt
#Using Probability Distribution given
def get_probability(delta_energy, Temperature):
return np.exp(-delta_energy / Temperature)
def get_energy(spins):
energy=0
for i in range(len(spins)):
energy=energy+interaction*spins[i-1]*spins[i]
energy= energy-field*sum(spins)
return energy
def delta_energy(spins,random_spin):
#If you do flip one random spin, the change in energy is:
#(By using a reduced formula that only involves the spin
# and its neighbours)
if random_spin==L-1:
PBC=0
else:
PBC=random_spin+1
old = -interaction*(spins[random_spin-1]*spins[random_spin] + spins[random_spin]*spins[PBC]) - field*spins[random_spin]
new = interaction*(spins[random_spin-1]*spins[random_spin] + spins[random_spin]*spins[PBC]) + field*spins[random_spin]
return new-old
def metropolis(L = 100, MC_samples=1000, Temperature = 1, interaction = 1, field = 0):
# intializing
#Spin Configuration
spins = np.random.choice([-1,1],L)
Beta = Temperature**(-1)
#Introducing Metropolis Hastings Algorithim
data = []
magnetization=[]
energy=[]
for i in range(MC_samples):
#Each Monte Carlo step consists in L random spin moves
for j in range(L):
#Choosing a random spin
random_spin=np.random.randint(0,L,size=(1))
#Compuing the change in energy of this spin flip
delta=delta_energy(spins,random_spin)
#Metropolis accept-rejection:
if delta<0:
#Accept the move if its negative
spins[random_spin]=-spins[random_spin]
#print('change')
else:
#If its positive, we compute the probability
probability=get_probability(delta,Temperature)
random=np.random.rand()
if random<=probability:
#Accept de move
spins[random_spin]=-spins[random_spin]
data.append(list(spins))
#Afer the MC step, we measure the system
magnetization.append(sum(spins)/L)
energy.append(get_energy(spins))
return data,magnetization,energy
def record_state_statistics(data,n=4):
ixs = tuple()
sub_sample = [[d[i] for i in range(n)] for d in data]
# get state number
state_nums = [int(sum([((j+1)/2)*2**i for j,i in zip(reversed(d),range(len(d)))])) for d in sub_sample]
return state_nums
# setting up problem
L = 200 # size of system
MC_samples = 1000 # number of samples
Temperature = 1 # "temperature" parameter
interaction = 1 # Strength of interaction between nearest neighbours
field = 0 # external field
# running MCMC
data = metropolis(L = L, MC_samples = MC_samples, Temperature = Temperature, interaction = interaction, field = field)
results = record_state_statistics(data[0],n=4) # I was also interested in the probability of each micro-state in a sub-section of the system
# Plotting
plt.figure(figsize=(20,10))
plt.subplot(2,1,1)
plt.imshow(np.transpose(data[0]))
plt.xticks([])
plt.yticks([])
plt.axis('tight')
plt.ylabel('Space',fontdict={'size':20})
plt.title('Critical dynamics in a 1-D Ising model',fontdict={'size':20})
plt.subplot(2,1,2)
plt.plot(data[2],'r')
plt.xlim((0,MC_samples))
plt.xticks([])
plt.yticks([])
plt.ylabel('Energy',fontdict={'size':20})
plt.xlabel('Time',fontdict={'size':20})
</code></pre>