求助：修正R语言实现Ising模型代码并绘制ggplot可视化图

需求说明 #

需要在R语言中实现Ising模型，修正存在逻辑错误的尝试代码（混有伪代码），并使用ggplot2结合scale_fill_viridis(option="D")绘制最终结果图。

原尝试R代码 #

L=4
N=L*L;N
A= matrix(nr= N, nc=L, 0);A
for (i in 0:N-1){
  A[i+1,1] = (i %/% L) * L + (i + 1) %% L
  A[i+1,2] = (i + L) %% N
  A[i+1,3] = (i %/% L) * L + (i - 1) %% L
  A[i+1,4] = (i - L) %% N
}

Temperature = 100
Spin = sample(c(1,-1),N,replace=TRUE)
spin_table = matrix(NA,L,L);spin_table
nsteps = 100 *L*L
for(step in 1:nsteps){
  for(n in A){
    k = sample(1,N-1,replace=TRUE)
    delta_E = 2.0 * Spin[k] * sum(Spin[n])
    if(runif(1) < exp(-delta_E/Temperature)){
      Spin[k] = Spin[k]*-1
    }
    }}
  

for(k in 1:N-1){
  y = k %/% L
  x = k %% L
  spin_table[x,y] = Spin[k]}
spin_table

参考可运行Python代码 #

import numpy as np
import pylab

L = 64
N = L*L

#   // floor division
#   %  modulus
neighbors = {i : ((i // L) * L + (i + 1) % L,   # RIGHT
                 (i + L) % N,                   # DOWN
                 (i // L) * L + (i - 1) % L,    # LEFT
                 (i - L) % N)                   # UP
                 for i in range(N)}

Temperature = 1.0
Spin = [np.random.choice([1, -1]) for k in range(N)]
spin_table = [[None for x in range(L)] for y in range(L)]

nsteps = N * 100
for step in range(nsteps):
    k = np.random.randint(0, N - 1)
    delta_E = 2 * Spin[k] * sum(Spin[n] for n in neighbors[k])
    if np.random.uniform(0.0, 1.0) < np.exp(-delta_E/Temperature):
        Spin[k] *= -1  # 原代码此处笔误，应为-1而非+1

for k in range(N):
    y = k // L
    x = k % L
    spin_table[x][y] = Spin[k]

pylab.close()
pylab.imshow(spin_table,extent=[0,L,0,L],interpolation='nearest')
pylab.title("Temperature: "+str(Temperature)+", Grid Size: "+str(L))
pylab.show()

修正后的R代码及可视化 #

关键修正点

• 索引适配：R采用1-based索引，需将Python的0-based邻居计算逻辑转换为R的索引规则
• 循环逻辑：原代码嵌套遍历邻居矩阵元素，正确逻辑应为每次迭代随机选择一个自旋位点
• 温度参数：原代码设为100（远高于2D Ising模型临界温度≈2.269），改为1.0以呈现有序态效果
• 矩阵填充：修正网格坐标的索引计算错误，确保自旋值正确映射到L×L网格
• 采样逻辑：sample(1:N-1)改为sample(1:N, 1)，保证能选中所有自旋位点

完整修正代码

library(ggplot2)
library(tidyr)

# 参数设置
L <- 64
N <- L * L
Temperature <- 1.0
nsteps <- N * 100

# 构建邻居矩阵（1-based索引）
neighbors <- matrix(0, nrow = N, ncol = 4)
for (i in 0:(N-1)) {
  idx <- i + 1
  # 右侧邻居（同一行循环边界）
  neighbors[idx, 1] <- ((i %/% L) * L + (i + 1) %% L) + 1
  # 下方邻居（同一列循环边界）
  neighbors[idx, 2] <- ((i + L) %% N) + 1
  # 左侧邻居（同一行循环边界）
  neighbors[idx, 3] <- ((i %/% L) * L + (i - 1) %% L) + 1
  # 上方邻居（同一列循环边界）
  neighbors[idx, 4] <- ((i - L) %% N) + 1
}

# 初始化自旋状态
Spin <- sample(c(1, -1), N, replace = TRUE)

# Metropolis算法迭代
for (step in 1:nsteps) {
  k <- sample(1:N, 1)
  delta_E <- 2.0 * Spin[k] * sum(Spin[neighbors[k, ]])
  if (runif(1) < exp(-delta_E / Temperature)) {
    Spin[k] <- -Spin[k]
  }
}

# 转换为L×L网格矩阵
spin_table <- matrix(NA, nrow = L, ncol = L)
for (k in 1:N) {
  i_0 <- k - 1  # 转0-base计算坐标
  y <- (i_0 %/% L) + 1
  x <- (i_0 %% L) + 1
  spin_table[y, x] <- Spin[k]
}

# 转换为ggplot兼容的长数据格式
spin_df <- as.data.frame(spin_table)
colnames(spin_df) <- 1:L
spin_df$y <- 1:L
spin_df <- pivot_longer(spin_df, cols = -y, names_to = "x", values_to = "spin")
spin_df$x <- as.integer(spin_df$x)

# 绘制结果图
ggplot(spin_df, aes(x = x, y = y, fill = spin)) +
  geom_tile() +
  scale_fill_viridis(option = "D", limits = c(-1, 1)) +
  labs(title = paste0("Ising Model: Temperature = ", Temperature, ", Grid Size = ", L),
       x = "", y = "") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text = element_blank())

补充说明 #

• 运行前需确保安装ggplot2和tidyr包，可通过install.packages(c("ggplot2", "tidyr"))完成安装
• 温度低于临界温度时，自旋会呈现明显的有序聚类；高于临界温度则呈随机分布

内容的提问来源于stack exchange，提问作者Homer Jay Simpson