基于PyTorch与四阶龙格-库塔的Kuramoto模型复现问题求助

我正在用四阶龙格-库塔（RK4）方法模拟Kuramoto模型，采用PyTorch做GPU加速。尝试复现论文（DOI: 10.1103/PhysRevE.65.026139）的图1以验证代码正确性，但输出结果存在偏差，虽未完全错误但无法证明代码有效。除论文采用Heun方法外，已设置相同条件。

以下是我的代码：

import numpy as np
import matplotlib.pyplot as plt
import torch

# Initial settings
sample_size = 100

t_max = 4 * 1e3
t_transient = 2 * 1e3 # transient datas are discarded
dt = 0.05
step_max = int(t_max / dt)
step_transient = int(t_transient / dt)

pi = torch.tensor(torch.pi, dtype=torch.float64, device='cuda')
i = torch.tensor(1j, dtype=torch.complex128, device='cuda')

N = 800

k_space = np.linspace(0.0, 5.0, 51)
r_plot = np.zeros_like(k_space)


# Kuramoto Model
def f(theta, w, K, A):

    phi = -theta + theta.reshape(-1, 1)
    return w + K * torch.sum(A * torch.sin(phi), axis = 0) / N


# Runge-Kutta 4th order
def RK4(f, theta, w, K, A):

    F1 = dt * f(theta, w, K, A)
    F2 = dt * f(theta + F1/2, w, K, A)
    F3 = dt * f(theta + F2/2, w, K, A)
    F4 = dt * f(theta + F3, w, K, A)

    return theta + (F1 + 2*F2 + 2*F3 + F4) / 6



r_t = torch.zeros(step_max - step_transient , dtype=torch.float64, device='cuda')
r_s = torch.zeros(sample_size, dtype=torch.float64, device='cuda')


for m, K in enumerate(k_space):
    for sample_num in range(sample_size):

        step = 0

        theta = torch.randn(N, dtype=torch.float64, device='cuda')
        theta = torch.remainder(theta, 2 * pi)
        # Mean field connection
        A = (torch.ones([N, N], dtype=torch.float64) - torch.eye(N, dtype=torch.float64)).to(device = 'cuda')
        # Natural frequencies from Gaussian
        w = torch.randn(N, dtype=torch.float64, device='cuda')


        while step < step_transient:

            theta = RK4(f, theta, w, K, A)
            theta = torch.remainder(theta, 2 * pi)
            step += 1

        while step < step_max:

            order_parameter = torch.sum(torch.exp(i * theta)) / N
            r_t[step - step_transient] = order_parameter.abs()

            theta = RK4(f, theta, w, K, A)
            theta = torch.remainder(theta, 2 * pi)
            step += 1

        r_s[sample_num] = r_t.mean()
        r_plot[m] = r_s.mean()


# drawing graphs and etc
plt.plot(k_space, r_plot, 'bo')
plt.show()

print(k_space)
print(r_plot)

已尝试：

• 调整张量索引顺序
• 使用二阶龙格-库塔（Heun方法）
• 分别尝试行向量与列向量

核心问题分析与修正方案 #

1. 初始相位分布不符合论文要求

论文明确初始相位是**[0, 2π)均匀分布**，但你当前用torch.randn生成正态分布后取余，这会导致初始相位集中在0附近，直接影响同步过程。

修正代码：

# 替换原初始相位生成代码
theta = torch.rand(N, dtype=torch.float64, device='cuda') * 2 * pi

2. 自然频率分布错误（关键差异）

论文中自然频率$\omega_i$服从标准柯西分布，而非你当前使用的高斯分布。柯西分布的重尾特性对Kuramoto模型的同步相变有决定性影响，这是结果偏差的核心原因。

标准柯西分布可通过逆变换采样生成，修正代码：

# 替换原自然频率生成代码
u = torch.rand(N, dtype=torch.float64, device='cuda')
w = torch.tan(pi * (u - 0.5))

3. 其他优化建议

• 提前定义统一设备变量，避免重复写device='cuda'：

  device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
  pi = torch.tensor(torch.pi, dtype=torch.float64, device=device)
  i = torch.tensor(1j, dtype=torch.complex128, device=device)
  

• 测试时可临时减小sample_size（如设为10）和N（如设为100），快速验证趋势是否符合论文图1：K<1时r≈0，K>1时r快速上升，K=5时r接近1。

修正上述两点后，结果应能与论文图1的趋势高度吻合，RK4方法的数值精度不会影响整体相变规律。

内容的提问来源于stack exchange，提问作者Waran