两版粒子正弦演化程序collchk结果不匹配原因排查

我编写了两个基于正弦函数演化粒子位置与动量的模拟程序，二者实现逻辑存在差异：

• 第一版程序在不同时间步更新对应时刻的三角函数系数，始终与初始粒子状态相乘完成计算
• 第二版程序采用固定时间步保持三角函数系数恒定，逐次更新参与计算的粒子位置、动量值

理论上两种实现得到的粒子轨迹应当完全一致，但实际输出的碰撞检测值collchk完全不匹配。

注：以下代码均为精简后的最小可复现版本。

第一版程序实现 #

import numpy as np
import matplotlib.pyplot as plt
from math import *
from numba import jit
@jit(nopython=True)
def f(SP, alf, dt, n):
    "Time"
    counter = 0;
    np.random.seed(0);
    Random=np.random.rand(n-1);
    C=np.array([cos(k*dt) for k in range(0,iter+1)])
    S=np.array([sin(k*dt) for k in range(0,iter+1)])
    for i in range(1, iter + 1):
        t = i * dt;
        Z = [];
        Up = [];
        Down = [];
        c,s=C[i],S[i]
        c1,s1=C[i-1],S[i-1]
        for j in range(n - 1):
            collchk=((c*(SP[j,0])+s*(SP[j,1]))-(c*(SP[j+1,0])+s*(SP[j+1,1])))*(c1*(SP[j,0])+s1*(SP[j,1])-(c1*(SP[j+1,0])+s1*(SP[j+1,1])));
            print(collchk)
    return  counter ,t
if __name__ == '__main__':
    n=5;
    dt=1/(10**(2));
    iter=(4);
    alf=sqrt(n);
    SP = np.array(sorted(np.array([  np.array([i,j,k]) for i, j,k in zip(Zinitial, Pinitial,SPIN)]),
                    key=lambda x: x[0]))
    counter,t = f(SP, alf, dt, n)
    print("Rate of collision per particle = ",(counter/(n*t)))

第一版程序输出 #

0.03103751937297854
0.08109145781934951
0.22901365345468785
2.5308437990270805
0.03078100154076592
0.08070470111986175
0.2263443458399873
2.538728781711264
0.030519417686482468
0.08030276710131726
0.22364568369773868
2.546117509808945
0.03025287244018221
0.07988581653196439
0.22091874645681717
2.553007027927399
Rate of collision per particle =  0.0

第二版程序实现 #

import numpy as np
import matplotlib.pyplot as plt
from math import *
from numba import jit
@jit(nopython=True)
def f(SP, SP1,alf, dt, n):
    "Time"
    counter = 0
    np.random.seed(0)
    Random=np.random.rand(n-1)
    c,s=cos(dt),sin(dt)
    for i in range(1, iter + 1):
      
        for j in range(n): 
          SP1[j,0]=SP[j,0]
          SP1[j,1]=SP[j,1]
          SP[j,0]=c*SP[j,0]+s*SP[j,1]
          SP[j,1]=c*SP[j,1]-s*SP[j,0]
        for j in range(n - 1):
            collchk=(SP[j,0]*SP[j+1,0])*(SP1[j,0]*SP[j+1,0])
            print(collchk)
    return counter ,t

if __name__ == '__main__':
    n=5;
    dt=1/(10**(2));
    iter=(4);
    alf=sqrt(n);
    Zinitial=[-0.9559708728305755, 1.5751906780581504,-0.014394517940369715, -0.7794354593838172, -0.49433681998221324]
    Pinitial=[0.4066095975477199,0.25889154830799843, -0.0046433361444705125, 0.33541895349394885, 0.270278187797722]
    SPIN=np.array([0, 1, 1, 0, 0])
    SP = np.array(sorted(np.array([  np.array([i,j,k]) for i, j,k in zip(Zinitial, Pinitial,SPIN)]),
                    key=lambda x: x[0]))
    SP1= np.array(sorted(np.array([  np.array([i,j,k]) for i, j,k in zip(Zinitial, Pinitial,SPIN)]),
                    key=lambda x: x[0]))
    counter,t = f(SP,SP1, alf, dt, n)
    print("Rate of collision per particle = ",(counter/(n*t)))
    print(t)   

第二版程序输出 #

0.5480084289400645
0.14618603324575455
5.067472238058949e-05
0.0005173929980597125
0.5383898895552246
0.14326676614418565
5.041820847142852e-05
0.0005221806114288456
0.5286869901202106
0.14033507300727344
5.013883366780392e-05
0.0005267892208272126
0.5189082113895205
0.13739341445190348
4.9836908449848145e-05
0.0005312143456136292
Rate of collision per particle =  0.0
0.04

第二版程序的核心错误 #

两版输出不匹配完全是第二版的实现逻辑错误导致，共有两处硬伤：

1. 粒子状态更新顺序错误
   正弦/余弦旋转更新位置和动量的正确逻辑是：新位置、新动量都必须基于更新前的旧位置、旧动量计算。但第二版代码里先更新了SP[j,0]（位置），紧接着计算新动量SP[j,1]时用的是已经更新过的新位置，而非原始旧位置，每一步迭代都会引入计算误差，迭代后粒子状态会完全偏离理论轨迹。
2. 碰撞检测公式完全不对等
   第一版的collchk计算逻辑是：相邻两个粒子t时刻的位置差，乘以二者t-1时刻的位置差，通过乘积符号判断时间步内是否发生粒子穿越碰撞。但第二版的collchk既没有计算相邻粒子的位置差，还混用了新旧状态的粒子索引，公式逻辑和第一版完全不一致，输出值自然不可能匹配。

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