蒙特卡洛模拟XY伊辛模型（python）
前言故事
世界上最早的通用电子计算机之一----ENIAC在发明后即被用于曼哈顿计划，乌拉姆敏锐地意识到在计算机的帮助下，可通过重复数百次模拟过程的方式来对概率变量进行统计估计。冯诺依曼立即认识到这个想法的重要性并给予支持。1947年，乌拉姆提出这种统计方法并应用于计算裂变的连锁反应。由于乌拉姆常说他的叔叔又在蒙特卡洛赌场输钱了，因此他的同事Nicolas Metropolis 戏称该方法为“蒙特卡洛”，不料却流传开去。
算法说明
参考文档：维基百科
1，把握这个算法，首先要体会到随机的意思。这是一个通过随机初始化，然后进行优化的一种算法。
2，直接上算法流程，之后再进行解释。
| 一，初始化：选择任意一个x0，假设我们的选择服从某种概率分布，比如高斯分布或者平均分布。 |
| 二，for i from 0 to N 开始循环迭代 |
| （1）Generate 产生一个候选x 。最简单的方式，随机一个 |
| （2）Calculate 计算接受比 alpha 这里定义的接受比等于 新产生的x在某种分布中的概率/x0在某种分布中的概率。这里的某种分布跟前面的高斯或者平均不是一回事，这里的分布更多的要体现问题的方面。这也是问题在这个算法中唯一的一次介入。 |
| （3）Accepted? |
| (3_1) 产生一个随机数u,其范围在[0,1] |
| (3_2) if u<=alpha 那么候选的x就去掉候选状况，正式替代 |
| （3_3）if u>alpha 那么候选的x无效 |
说明:1,为什么u<alpha才有效替换？
我们可以看到，当p(before)=p(after)时，这时，按照我们的常人观点，此时变不变是一半对一半！因为两者的概率一样，所以发生的概率应该一样的。但是这里并不是，如果发生了这种事，直接替换。
当p(after)<p(before)时，此时变化后的状态存在的概率要低于不变的概率，此时并不是一定不变，仍然以一定的概率变。
说明2，第（2）步的计算概率是最跟问题有关的。前面的都是随机分布变化，跟问题一点关系没有。请问，问题的概率分布一定存在吗？如果存在，一般怎么构造这种概率表达式？
我也没接触太多。目前，只见到过这个伊辛模型。在热力学统计中，确实有一个物理量对应这种概率分布是多少，称为配分函数。这里是配分函数的统计意义解释。
代码实现
2024 copilot 优化：
import numpy as np
import matplotlib.pyplot as plt
def get_coupling_energy(row, col, grid):
"""
Calculate the coupling energy of a single cell in a grid.
Parameters:
row (int): The row index of the cell.
col (int): The column index of the cell.
grid (numpy.ndarray): A 2D array representing the grid, where each element is an angle.
Returns:
float: The coupling energy of the specified cell.
"""
# 计算单个格子的能量
angle = grid[row][col]
neighbors = [(row-1, col), (row+1, col), (row, col-1), (row, col+1)]
energy = 0
for neighbor in neighbors:
neighbor_row, neighbor_col = neighbor
if 0 <= neighbor_row < grid.shape[0] and 0 <= neighbor_col < grid.shape[1]:
energy += -np.cos(angle - grid[neighbor_row][neighbor_col])
return energy
def energy_logger(func):
"""
Decorator to log the energy calculated by the function.
Args:
func (function): The function to be decorated.
Returns:
function: The wrapped function with logging.
"""
def wrapper(*args, **kwargs):
energy = func(*args, **kwargs)
print(f"Calculated energy: {energy}")
return energy
return wrapper
@energy_logger
def calculate_all_energy(grid):
"""
Calculate the average energy of the entire grid.
This function computes the total energy of the grid by summing up the coupling
energy for each cell in the grid. The average energy is then calculated by
dividing the total energy by the number of cells in the grid and further
divided by 2.
Args:
grid (numpy.ndarray): A 2D array representing the grid where each
element is an angle.
Returns:
float: The average energy of the grid.
"""
# 计算整个格子的能量，为了求平均内能
total_energy = 0
for row in range(len(grid)):
for col in range(len(grid[0])):
total_energy += get_coupling_energy(row, col, grid)
average_energy = total_energy / (len(grid[0]) * len(grid))
return average_energy / 2
def metropolis_hastings(grid, temperature, num_iterations):
"""
Perform the Metropolis-Hastings algorithm on a given grid.
Parameters:
grid (numpy.ndarray): A 2D array representing the initial state of the system.
temperature (float): The temperature parameter controlling the acceptance probability.
num_iterations (int): The number of iterations to perform the algorithm.
Returns:
numpy.ndarray: The updated grid after performing the Metropolis-Hastings algorithm.
"""
# Metropolis-Hastings 算法的实现
for _ in range(num_iterations):
row = np.random.randint(0, grid.shape[0])
col = np.random.randint(0, grid.shape[1])
old_energy = get_coupling_energy(row, col, grid)
old_angle = grid[row, col]
grid[row, col] = 2 * np.pi * np.random.random()
new_energy = get_coupling_energy(row, col, grid)
if new_energy > old_energy and np.exp((old_energy - new_energy) / temperature) < np.random.random():
grid[row, col] = old_angle
return grid
def plot_grid(grid, title_suffix = ""):
X, Y = np.meshgrid(np.arange(0, grid.shape[0]), np.arange(0, grid.shape[0]))
U = np.cos(grid)
V = np.sin(grid)
plt.figure()
plt.quiver(X, Y, U, V, units='inches')
plt.title(f'Ising XY Model Simulation_{title_suffix}')
plt.show(block=False) # 确保显示图表
def initialize_grid(grid_size):
# 初始化格子
return 2 * np.pi * np.random.random((grid_size, grid_size))
def simulate_single_sample(grid, grid_size, temperature):
# 模拟单个样本
new_grid = np.array(grid)
for _ in range(100):
new_grid = metropolis_hastings(new_grid, temperature, grid_size**2)
average_energy = calculate_all_energy(new_grid)
return new_grid, average_energy
def simulate_ising_model(grid_size, temperature):
"""
Simulates the Ising model for a given grid size and temperature.
Args:
grid_size (int): The size of the grid (number of rows and columns).
temperature (float): The temperature at which to simulate the Ising model.
Returns:
tuple: A tuple containing:
- new_grid (numpy.ndarray): The grid after simulation.
- grid (numpy.ndarray): The initial grid before simulation.
- initial_energy (float): The average energy of the initial grid.
- final_energy (float): The average energy of the final grid.
"""
# 初始化样本
grid = initialize_grid(grid_size)
initial_energy = calculate_all_energy(grid)
new_grid, final_energy = simulate_single_sample(grid, grid_size, temperature)
return new_grid, grid, initial_energy, final_energy
# 主程序部分
if __name__ == "__main__":
grid_size = 30
temperature = 0.5
final_grid, initial_grid,initial_energy, final_energy = simulate_ising_model(grid_size, temperature)
plot_grid(initial_grid, "Init")
plot_grid(final_grid,"Final") # 绘制最终的格子图原始代码：
import random
import matplotlib.pyplot as plt
import numpy as np
import copy
import math
import time
'''这个就是MCMC模拟，用来模拟某个温度下XY Ising模型分布，
几点注意：
注意1，二维伊辛模型，我们用矩阵来模拟。
注意2，旋转的方向，我们用0到2pi表示吧
算法过程：
一，用一个对称的分布，高斯分布初始化矩阵
二，下面是循环
（1）产生一个候选的自旋方向，为了连续，假设新产生的自旋方向变化最多比原来变化pi/2
也就是旋转90度
（2）计算两个概率，这里热统中的配分函数正比于概率，因此我们用配分函数Z的比值
Z(变化后)/Z(变化前)=exp(-(E后-E前)/kT) ,这里k是玻尔兹曼常数，T是开尔文温度。令
这个值是alpha
(3)判断是都接受*********这个规则有进一步讨论*************
（3_1)产生一个随机数u在【0,1】之间
（3_2）如果u<=alhpa,接受候选的，改变此时的自旋状态
(3_3)如果u>alpha，不接受候选的，不改变此时的自旋状态
inputdata: S :matrix 随机分布,假设已经产生
param: T 温度
delta 最大的变化度数，默认是90度，也可以调整为其他
outputdata:S
'''
def MetropolisHastings(S,T,numsOfItera):
deltamax=0.5
k = 1 #玻尔兹曼常数
for sdw in range(numsOfItera):
# k2 = np.random.randint(0,S.shape[0]**2)
i = np.random.randint(0,S.shape[0])
j = np.random.randint(0,S.shape[0])
# print('产生的随机位置是：',i,j)
#time.sleep(0.1)
for m in range(1):
delta = (2*np.random.random()-1)*deltamax
newAngle = S[i][j]+delta
# print(delta)
energyBefore = getEnergy(i=i,j=j,S=S,angle=None)
energyLater = getEnergy(i,j,S=S,angle=newAngle)
alpha = math.exp(-(energyLater-energyBefore)/(k*T))
#print(alpha)
# if alpha>=1:
# print('大于1的哦')
if alpha >=1:
S[i][j]=newAngle
elif np.random.uniform(0.0,1.0)<=1.0*alpha:
S[i][j]=newAngle
return S
#计算i,j位置的能量 = 与周围四个的相互能之和
def getEnergy(i,j,S,angle=None):
width = S.shape[0]
height = S.shape[1]
# print('矩阵的宽和高是',width,height)
top_i = i-1 if i>0 else width-1
bottom_i = i+1 if i<(width-1) else 0
left_j = j-1 if j>0 else height-1
right_j = j+1 if j<(height-1) else 0
enviroment = [[top_i,j],[bottom_i,j],[i,left_j],[i,right_j]]
# print(i,j,enviroment)
#print(enviroment)
energy = 0
if angle ==None:
# print('angle==None')
for num_i in range(0,4,1):
energy += -np.cos(S[i][j]-S[enviroment[num_i][0]][enviroment[num_i][1]])
else:
# print('Angle')
for num_i in range(0,4,1):
energy += -np.cos(angle-S[enviroment[num_i][0]][enviroment[num_i][1]])
return energy
#S=2*np.pi*np.random.rand(30,30)
#计算整个格子的能力，为了求平均内能
def calculateAllEnergy(S):
energy =0
for i in range(len(S)):
for j in range(len(S[0])):
energy +=getEnergy(i,j,S)
averageEnergy = energy/(len(S[0])*len(S))
return averageEnergy/2
#print(S)
#for j in range(1000):
# print(j)
# MetropolisHastings(S,10)
#这个是输入样本的多少，格子的尺寸，温度。中间那个循环，是随机取迭代的过程
def getWeightValue(numsOfSample,sizeOfSample,temperature):
for i in range(numsOfSample): #产生个数
print('+++++++正在计算第%s个样本++++++++++'%i)
S=2*np.pi*np.random.random((sizeOfSample,sizeOfSample))
initialEnergy = calculateAllEnergy(S)
print('系统的初始能量是:',initialEnergy)
newS = np.array(copy.deepcopy(S))
for nseeps in range(100):
newS = MetropolisHastings(newS,temperature,sizeOfSample**2)
aveEnergy = calculateAllEnergy(newS)
plot(newS)
print('系统的平均能量是:',aveEnergy)
reshaped = np.reshape(newS,(1,sizeOfSample**2))
if i ==0:
s = copy.deepcopy(reshaped)
continue
else:
s = np.row_stack((s,reshaped))
return s
#运行getweightValue函数，中间已经把结果会成图了
res = getWeightValue(1,40,2)
#print(len(res))
#画成箭头图表示出现
def plot(S):
X, Y = np.meshgrid(np.arange(0,S.shape[0]),np.arange(0,S.shape[0]))
U = np.cos(S)
V = np.sin(S)
plt.figure()
#plt.title('Arrows scale with plot width, not view')
Q = plt.quiver(X, Y, U, V, units='inches')
#qk = plt.quiverkey(Q, 0.3, 0.3, 1, r'$2 \frac{m}{s}$', labelpos='E',
# coordinates='figure')
plt.show()结果展示：
系统值是：
+++++++正在计算第0个样本++++++++++
系统的初始能量是: 0.006919117595
系统的平均能量是: -0.494289314266
+++++++正在计算第0个样本++++++++++
系统的初始能量是: 0.154688600481
系统的平均能量是: -0.315447987184
+++++++正在计算第0个样本++++++++++
系统的初始能量是: 0.00142378493144
系统的平均能量是: -0.588001265121跟点击打开链接对比，
我们看到在温度等于2的时候，还是很类似的。
其他温度的结果：
R 语言，不保证正确。
library(ggplot2)
# 计算单个格子的能量
get_coupling_energy <- function(row, col, grid) {
angle <- grid[row, col]
neighbors <- list(c(row-1, col), c(row+1, col), c(row, col-1), c(row, col+1))
energy <- 0
for (neighbor in neighbors) {
neighbor_row <- neighbor[1]
neighbor_col <- neighbor[2]
if (neighbor_row > 0 && neighbor_row <= nrow(grid) && neighbor_col > 0 && neighbor_col <= ncol(grid)) {
energy <- energy - cos(angle - grid[neighbor_row, neighbor_col])
}
}
return(energy)
}
# 计算整个格子的能量，为了求平均内能
calculate_all_energy <- function(grid) {
total_energy <- 0
for (row in 1:nrow(grid)) {
for (col in 1:ncol(grid)) {
total_energy <- total_energy + get_coupling_energy(row, col, grid)
}
}
average_energy <- total_energy / (nrow(grid) * ncol(grid))
return(average_energy / 2)
}
# Metropolis-Hastings 算法的实现
metropolis_hastings <- function(grid, temperature, num_iterations) {
for (i in 1:num_iterations) {
row <- sample(1:nrow(grid), 1)
col <- sample(1:ncol(grid), 1)
old_energy <- get_coupling_energy(row, col, grid)
old_angle <- grid[row, col]
grid[row, col] <- runif(1, 0, 2 * pi)
new_energy <- get_coupling_energy(row, col, grid)
if (new_energy > old_energy && exp((old_energy - new_energy) / temperature) < runif(1)) {
grid[row, col] <- old_angle
}
}
return(grid)
}
# 初始化格子
initialize_grid <- function(grid_size) {
return(matrix(runif(grid_size * grid_size, 0, 2 * pi), nrow = grid_size, ncol = grid_size))
}
# 模拟单个样本
simulate_single_sample <- function(grid, grid_size, temperature) {
new_grid <- grid
for (i in 1:100) {
new_grid <- metropolis_hastings(new_grid, temperature, grid_size^2)
}
average_energy <- calculate_all_energy(new_grid)
return(list(new_grid = new_grid, average_energy = average_energy))
}
# 模拟 Ising 模型
simulate_ising_model <- function(grid_size, temperature) {
grid <- initialize_grid(grid_size)
initial_energy <- calculate_all_energy(grid)
result <- simulate_single_sample(grid, grid_size, temperature)
new_grid <- result$new_grid
final_energy <- result$average_energy
return(list(new_grid = new_grid, grid = grid, initial_energy = initial_energy, final_energy = final_energy))
}
# 绘制格子图
plot_grid <- function(grid, title) {
grid_df <- expand.grid(x = 1:nrow(grid), y = 1:ncol(grid))
grid_df$angle <- as.vector(grid)
grid_df$u <- cos(grid_df$angle)
grid_df$v <- sin(grid_df$angle)
ggplot(grid_df, aes(x = x, y = y, angle = angle)) +
geom_segment(aes(xend = x + u, yend = y + v), arrow = arrow(length = unit(0.2, "cm"))) +
coord_fixed() +
ggtitle(title) +
theme_minimal()
}
# 主程序部分
grid_size <- 30
temperature <- 0.5
result <- simulate_ising_model(grid_size, temperature)
initial_grid <- result$grid
final_grid <- result$new_grid
initial_energy <- result$initial_energy
final_energy <- result$final_energy
print(paste("Initial Energy:", initial_energy))
print(paste("Final Energy:", final_energy))
plot_grid(initial_grid, "Initial Grid")
plot_grid(final_grid, "Final Grid")Fortran语言，不保证正确
module ising_model
implicit none
contains
! 计算单个格子的能量
real function get_coupling_energy(row, col, grid, grid_size)
integer, intent(in) :: row, col, grid_size
real, dimension(grid_size, grid_size), intent(in) :: grid
real :: angle
integer :: neighbor_row, neighbor_col
real :: energy
integer :: i
angle = grid(row, col)
energy = 0.0
do i = 1, 4
select case (i)
case (1)
neighbor_row = row - 1
neighbor_col = col
case (2)
neighbor_row = row + 1
neighbor_col = col
case (3)
neighbor_row = row
neighbor_col = col - 1
case (4)
neighbor_row = row
neighbor_col = col + 1
end select
if (neighbor_row >= 1 .and. neighbor_row <= grid_size .and. neighbor_col >= 1 .and. neighbor_col <= grid_size) then
energy = energy - cos(angle - grid(neighbor_row, neighbor_col))
end if
end do
get_coupling_energy = energy
end function get_coupling_energy
! 计算整个格子的能量，为了求平均内能
real function calculate_all_energy(grid, grid_size)
integer, intent(in) :: grid_size
real, dimension(grid_size, grid_size), intent(in) :: grid
real :: total_energy
integer :: row, col
total_energy = 0.0
do row = 1, grid_size
do col = 1, grid_size
total_energy = total_energy + get_coupling_energy(row, col, grid, grid_size)
end do
end do
calculate_all_energy = total_energy / (grid_size * grid_size) / 2.0
end function calculate_all_energy
! Metropolis-Hastings 算法的实现
subroutine metropolis_hastings(grid, temperature, num_iterations, grid_size)
integer, intent(in) :: num_iterations, grid_size
real, dimension(grid_size, grid_size), intent(inout) :: grid
real, intent(in) :: temperature
integer :: i, row, col
real :: old_energy, new_energy, old_angle, rand_num
do i = 1, num_iterations
call random_number(rand_num)
row = int(rand_num * grid_size) + 1
call random_number(rand_num)
col = int(rand_num * grid_size) + 1
old_energy = get_coupling_energy(row, col, grid, grid_size)
old_angle = grid(row, col)
call random_number(rand_num)
grid(row, col) = 2.0 * 3.141592653589793 * rand_num
new_energy = get_coupling_energy(row, col, grid, grid_size)
if (new_energy > old_energy) then
call random_number(rand_num)
if (exp((old_energy - new_energy) / temperature) < rand_num) then
grid(row, col) = old_angle
end if
end if
end do
end subroutine metropolis_hastings
! 初始化格子
subroutine initialize_grid(grid, grid_size)
integer, intent(in) :: grid_size
real, dimension(grid_size, grid_size), intent(out) :: grid
integer :: row, col
real :: rand_num
do row = 1, grid_size
do col = 1, grid_size
call random_number(rand_num)
grid(row, col) = 2.0 * 3.141592653589793 * rand_num
end do
end do
end subroutine initialize_grid
! 模拟单个样本
subroutine simulate_single_sample(grid, new_grid, average_energy, grid_size, temperature)
integer, intent(in) :: grid_size
real, dimension(grid_size, grid_size), intent(in) :: grid
real, dimension(grid_size, grid_size), intent(out) :: new_grid
real, intent(out) :: average_energy
integer :: i
new_grid = grid
do i = 1, 100
call metropolis_hastings(new_grid, temperature, grid_size**2, grid_size)
end do
average_energy = calculate_all_energy(new_grid, grid_size)
end subroutine simulate_single_sample
! 模拟 Ising 模型
subroutine simulate_ising_model(new_grid, grid, initial_energy, final_energy, grid_size, temperature)
integer, intent(in) :: grid_size
real, intent(in) :: temperature
real, dimension(grid_size, grid_size), intent(out) :: new_grid, grid
real, intent(out) :: initial_energy, final_energy
call initialize_grid(grid, grid_size)
initial_energy = calculate_all_energy(grid, grid_size)
call simulate_single_sample(grid, new_grid, final_energy, grid_size, temperature)
end subroutine simulate_ising_model
end module ising_model
program main
use ising_model
implicit none
integer :: grid_size
real :: temperature
real, dimension(:,:), allocatable :: initial_grid, final_grid
real :: initial_energy, final_energy
grid_size = 30
temperature = 0.5
allocate(initial_grid(grid_size, grid_size))
allocate(final_grid(grid_size, grid_size))
call simulate_ising_model(final_grid, initial_grid, initial_energy, final_energy, grid_size, temperature)
print *, "Initial Energy:", initial_energy
print *, "Final Energy:", final_energy
! 绘制格子图的代码可以使用外部工具，例如 Python 的 matplotlib 或者 gnuplot
! 这里我们只输出数据，可以使用其他工具进行绘图
open(unit=10, file='initial_grid.dat', status='replace')
write(10, *) initial_grid
close(10)
open(unit=20, file='final_grid.dat', status='replace')
write(20, *) final_grid
close(20)
end program main
扫一扫
1701
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
