如何解决MCMC方法1D铁磁Ising模型
我的问题与使用马尔可夫链蒙特卡洛方法(MCMC)的一维Ising模型的Python编码有关。
我有以下哈密顿量
$$H = - \sum_{i=1}^{L-1}\sigma_{i}sigma_{i+1} - B\sum_{i=1}^{L}\sigma_{i}$$
我想编写一个Python函数,该函数生成一个马尔可夫链,在每个步骤中,它都会计算并保存磁化强度(每个位置)和能量。
能量为(=哈密顿量),我将磁化强度定义为:
$$\frac{1}{L}\sum_{i}\sigma_{i}$$
我的概率分布为:
$$p(x) = e^{-H\beta}$$ where,$T^{-1} = \beta$
对于马尔可夫链,我将实施Metropolis-Hastings算法;
if $$\frac{P(\sigma')}{P(\sigma)} = e^{(H(\sigma)-H(\sigma'))\beta}$$
我的想法是在以下情况下接受过渡
$$H(\sigma') < H(\sigma)$$
并且仅接受转换
$$H(\sigma') > H(\sigma)$$
具有可能性
$$P = e^{(H(\sigma)-H(\sigma'))\beta}$$
所以让我设置一些参数,例如:
$L=20$ - Lattice Size
$T=2$ - Temperature
$B=0$ - Magnetic Field
在计算之后,我需要绘制磁化强度和能量与步长大小的直方图。我对这部分没有问题。
我的python知识不是很好,但是我包括了我的草稿(未完成)。我认为我没有什么进步。任何帮助都会很棒。
#Coding attempt MCMC 1-Dimensional Ising Model
import numpy as np
import matplotlib.pyplot as plt
#Shape of Lattice L
L = 20
Shape = (20,20)
#Spin Configuration
spins = np.random.choice([-1,1],Shape)
#Magnetic moment
moment = 1
#External magnetic field
field = np.full(Shape,0)
#Temperature
Temperature = 2
Beta = Temperature**(-1)
#Interaction (ferromagnetic if positive,antiferromagnetic if negative)
interaction = 1
#Using Probability Distribution given
def get_probability(Energy1,Energy2,Temperature):
return np.exp((Energy1 - Energy2) / Temperature)
def get_energy(spins):
return -np.sum(
interaction * spins * np.roll(spins,1,axis=0) +
interaction * spins * np.roll(spins,-1,axis=1) +
interaction * spins * np.roll(spins,axis=1)
)/2 - moment * np.sum(field * spins)
#Introducing Metropolis Hastings Algorithim
x_now = np.random.uniform(-1,1) #initial value
d = 10**(-1) #delta
y = []
for i in range(L-1):
#generating next value
x_proposed = np.random.uniform(x_now - d,x_now + d)
#accepting or rejecting the value
if np.random.rand() < np.exp(-np.abs(x_proposed))/(np.exp(-np.abs(x_now))):
x_now = x_proposed
if i % 100 == 0:
y.append(x_proposed)
解决方法
在这里,我更改了您的代码以一如既往地解决问题。
请仔细检查代码和公式
#Coding attempt MCMC 1-Dimensional Ising Model
import numpy as np
import matplotlib.pyplot as plt
#Shape of Lattice L
L = 20
#Shape = (20)
#Number of Monte Carlo samples
MC_samples=1000
#Spin Configuration
spins = np.random.choice([-1,1],L)
print(spins)
#Magnetic moment
moment = 1
#External magnetic field
field = 0
#Temperature
Temperature = 2
Beta = Temperature**(-1)
#Interaction (ferromagnetic if positive,antiferromagnetic if negative)
interaction = 1
#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(L):
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:
PBC=0
else:
PBC=random_spin+1
return -2*interaction*(spins[random_spin-1]*spins[random_spin]+
spins[random_spin]*spins[PBC]+field*spins[random_spin])
#Introducing Metropolis Hastings Algorithim
#x_now = np.random.uniform(-1,1) #initial value
#d = 10**(-1) #delta
#y = []
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(L-1,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]
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]
#Afer the MC step,we measure the system
magnetization.append(sum(spins)/L)
energy.append(get_energy(spins))
print(magnetization,energy)
#Do histograms and plots
在模拟结束时,变量磁化强度和能量是包含每个MC步骤的测量值的数组。 您可以直接使用这些数组来计算直方图和图。
注意:能量数组是系统的总能量,而不是能量/ L。
,我一直在寻找1D Ising模型的简单实现,并且遇到了这篇文章。虽然我不是该领域的专家,但是我确实写了一个相关主题的硕士。我在Oriol Cabanas Tirapu的答案中实现了代码,发现了一些错误(我认为)。
下面是我的改编版,哦,他们的代码。希望它对某人有用。
const testString = 'This #here is a #test.';
console.log(testString.match(/[#_]?\w+|[^\w#_]+/g));
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。