如何用scipy.integrate.RK45模拟四粒子向最近邻相向运动？

核心逻辑：运动方程的编写 #

首先，状态向量需要包含所有四个粒子的位置（因为每个粒子的运动依赖于其他粒子的位置），而非单个粒子的位置+速度。我们可以将状态向量y整理为[x0, y0, x1, y1, x2, y2, x3, y3]，其中(xi, yi)是第i个粒子的坐标。

运动方程motion_eq的任务是计算每个粒子位置的变化率（即速度），步骤如下：

1. 将状态向量拆分为四个粒子的坐标矩阵；
2. 对每个粒子，找到当前距离最近的邻居；
3. 计算指向最近邻居的单位方向向量，乘以恒定速度V得到速度分量；
4. 将所有速度分量展平后返回。

具体实现代码：

import numpy as np
from scipy.integrate import RK45
import matplotlib.pyplot as plt

V = 1.0  # 恒定速度大小
a = 2.0  # 正方形边长
meet_threshold = 1e-3  # 粒子相遇的距离阈值

def motion_eq(t, y):
    # 将状态向量拆分为4个粒子的坐标 (4, 2)
    positions = y.reshape(-1, 2)
    velocity = np.zeros_like(positions)
    
    for i in range(4):
        # 计算当前粒子到其他所有粒子的距离
        distances = np.linalg.norm(positions - positions[i], axis=1)
        # 排除自身（距离为0的情况）
        distances[i] = np.inf
        # 找到最近邻居的索引
        nearest_idx = np.argmin(distances)
        # 计算指向最近邻居的单位方向向量
        dir_vec = positions[nearest_idx] - positions[i]
        dir_vec /= np.linalg.norm(dir_vec)
        # 速度 = 恒定速度大小 * 单位方向向量
        velocity[i] = V * dir_vec
    
    # 展平速度向量作为返回值（对应状态向量的导数）
    return velocity.flatten()

初始化与积分设置 #

需要修改初始化代码，设置四个粒子的初始位置（正方形四个角落），并初始化RK45求解器：

# 初始位置：正方形四个角落（中心在原点）
initial_pos = np.array([
    [-a/2, -a/2],
    [a/2, -a/2],
    [a/2, a/2],
    [-a/2, a/2]
])
# 状态向量初始值：展平的初始位置
initial_y = initial_pos.flatten()

dt = 0.1
motion_solver = RK45(
    motion_eq,
    t0=0,
    y0=initial_y,
    t_bound=np.inf,
    first_step=dt,
    max_step=dt
)

可视化与停止逻辑 #

更新可视化函数，添加粒子相遇的停止判断，并同时更新四个粒子的位置：

# 初始化绘图
fig, ax = plt.subplots()
ax.set_aspect('equal')
ax.set_xlim([-a, a])
ax.set_ylim([-a, a])
# 绘制四个粒子，保存每个粒子的Line2D对象
particles = []
for pos in initial_pos:
    p, = ax.plot(pos[0], pos[1], 'o')
    particles.append(p)

def update():
    global motion_solver
    # 执行一步积分
    motion_solver.step()
    # 获取当前所有粒子的位置
    current_pos = motion_solver.y.reshape(-1, 2)
    # 更新每个粒子的显示位置
    for p, pos in zip(particles, current_pos):
        p.set_data(pos[0], pos[1])
    # 检查是否有粒子相遇
    for i in range(4):
        for j in range(i+1, 4):
            dist = np.linalg.norm(current_pos[i] - current_pos[j])
            if dist < meet_threshold:
                timer.stop()
                print(f"粒子{i}和{j}在t={motion_solver.t:.2f}时相遇")
                return
    plt.draw()

# 设置定时器更新画面
timer = fig.canvas.new_timer(interval=50)
timer.add_callback(update)
timer.start()

plt.show()

关键说明 #

• 对称性简化：如果利用四个粒子始终保持正方形对称的特性，可以不用遍历找最近邻居，直接让每个粒子追顺时针/逆时针的下一个粒子（比如粒子0追1，1追2，2追3，3追0），这样计算效率更高，结果和“找最近邻居”一致（因为对称情况下每个粒子的最近邻居就是下一个粒子）。
• 停止条件：用距离阈值判断相遇，避免数值计算中无法完全等于0的问题。
• 状态向量设计：这里没有把速度存入状态向量，因为速度完全由当前位置决定，只需计算位置的变化率即可，符合ODE求解的要求。

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