多区域N体问题解析函数的近似计算

刘 超，章社生, 吴永红 (武汉理工大学理学院，湖北 武汉 430063)

多区域N体问题解析函数的近似计算

刘 超，章社生, 吴永红 (武汉理工大学理学院，湖北 武汉 430063)

研究了粒子相互作用的多区域N体问题解析函数的近似计算，利用多重积分和delta函数，求得径向分布函数的解析表达式，推导出双球域的分析解。计算结果表明，当粒子数为2000个时，计算误差小于千分之一，且计算速度快。

粒子模拟；N体问题；径向分布函数；多重积分；分析解；数值解

粒子模拟已经成为车路协同等科学领域的重要研究工具[1]，它将复杂系统(生物大分子、星系、车路协同等)看成是大量经典物理可描述的单个个体(原子、星体等 )的集合，然后对所有这些个体特征数据进行分析。笔者专注于一种在粒子模拟数据分析中占有极其重要地位的双体关联函数查询-空间距离直方图(SDH)的研究[2]，并且关注Barnes2Hut算法(BHA)和快速多极子方法( FMM)[3]。SDH是一种被称为径向分布函数(RDF)的连续统计分布函数的直接估计[4-5]，与该系统的结构因子(structure factor)相关联[6]。当研究区域为单连通时，文献[7]求出了粒子相互作用径向分布函数的解析表达式，该工作能推广到互不相连的多区域N体问题。

1 多区域粒子径向分布

径向分布函数(RDF)[7]的定义为：

(1)

式中，Fh(r)是在粒子之间距离为r，厚度为h的球壳内粒子距离数量，称为距离和函数。设粒子坐标为Xi，i=1,…,N，则距离和函数Fh(r)可表示为：

(2)

式中，r*=|Xi-Xj|为Xi与Xj之间的距离；当|r*-r|≤h/2时函数δh=1，当|r*-r|>h/2时函数δh=0。式(2)中除以2是为了消去重复累加。

用一个立方体Ω包含所有粒子。按粒子分布的稀疏情况，将立方体Ω划分为N′个子区域Ωk(k=1,…,N′),则距离和函数Fh(r)可表示为：

(3)

式中,Fi,j为第i个区域中的粒子与第j个区域中的粒子之间的距离和函数;当j=i时表示第i个区域内部粒子之间的距离和函数,当j≠i时表示第i个区域中的粒子与第j个区域中的粒子之间的距离和函数。

2 双球域分析解

考虑粒子分布在2个球域O和O′，球内分别有N和N′个粒子，对应粒子密度分别为ρ和ρ′(粒子数除球体积)，两球心距离为a,半径分别为R和R′。取粒子距离和函数为:

Fh(r)=F11(r)+F12(r)+F22(r)

(4)

图1 双球域对数示意图

式中,r=mh,h为步长，m=1,2,3,…;F11(r)为球O内部粒子产生的距离和函数;F22(r)为球O′内部粒子产生的距离和函数;F12(r)为球O中的粒与球O′中的粒子相互作用产生的距离和函数;F11(r)和F22(r)用文献[7]中的公式计算。设点Y为球O′表面点，它到点O的距离为L(见图1)，则点Y与球O中粒子距离和函数为：

(5)

根据空间距离直方图(SDH)的定义[2]，将球域O划分为M′个子区域dΩ，第i个子区域dΩi内的粒子数为n(i)，子区域的中点坐标为Yi，第i个子区域中点和区域外点Y的距离为|Yi-Y|，用n(i)|Yi-Y|近似表示第i个子区域内所有粒子与在Y点粒子的距离的和。将下标i从1取到M′，然后求和， 再根据多维积分的定义，当dV非常小时，有:

Y(r)=0rL+R

(6)

式(6)只与距离L和R有关，与粒子坐标Y无关。容易求得，球心在O点，半径为L的球与球O′的中粒子相互距离为r的距离函数：

(7)

式中:

3 算例及误差分析

例1设2球半径R=R′=1，球心相距a=3, 每个球内随机均匀分布M个粒子,研究粒子的分布。

例1中2球之间的距离大于球的半径，它们能放置于半径为2.5的大球中,但用半径为2.5的大球的粒子径向分布解析公式计算将带来较大的误差。下面笔者用双球粒子径向分布解析公式计算。

由于半径R=1， 则有:

(8)

且:

此时，只有2种情况:①r≤a,(L1,L2)=(a-1,r+1); ②r>a;(L1,L2)=(r-1,a+1)。

取粒子数M=500，半径R=R′=1，两球心距离a=3，步长h=0.1；球中粒子点用随机方法产生。对应的计算值绘于图2。图2中实线为理论值，星号为直接计算值(用2点之间距离公式计算)。图2中有2个峰值，第1个峰在r=1附近，由球内粒子相互作用产生。第2个峰值在r=3, 它由球与球中的粒子相互作用产生。由图2可知，当粒子仅为M=500时，理论解与直接法计算结果接近。进一步计算表明，当球中粒子数增加时，能减少理论计算误差。

定义平均误差为Er:

(9)

表1 平均误差Er和直接法CPU计算时间t随粒子数M的变化

图2 2球距离和函数随r的变化

当R=1、R′=1、a=3、h=0.1时，平均误差Er和CPU计算时间t如表1所示。表1中M为球内粒子数，计算平均误差Er乘了1000，计算时间t为直接法的CPU计算时间。由表1可知，当粒子数M=100时，平均误差小于2‰，粒子数增加后，平均误差总体上减少，但发现有波动。原因是粒子坐标点是用随机方法选取的，当粒子总数M改变后，粒子坐标都改变。由表1可知，随着粒子数增加，CPU计算时间增加很快。

4 结 语

笔者推广了已有的理论[7],构造了多区域解析计算模型，得到了近似计算径向距离的积分表达式。在粒子密度为均匀分布的条件下，得到双球域积分解。实例计算结果表明，积分解有很高的精确度，且计算速度快。

[1]胡旻，祝大军，刘大刚，等.粒子模拟软件吸收边界的研究[J].强激光与粒子束，2006，18(8)：1315-1318.

[2] Tu Yi-cheng, Chen Shao-ping, Sagar Pandit.Computing Distance Histograms Efficiently in Scientific Databases[A].In procedings of ICDE[C]. 2009：796-807.

[3]王武，冯仰德，迟学斌.树结构在N体问题中的应用[J].计算应用研究，2008，125(1)：42-44.

[4] Bamdad M, Alavi S, Najafi B,et al.A new expression for radial distribution function and infinite shear modulus of Lennard Jones fluids [J].Chem Phys, 2006, 325：554-562.

[5]陈念贻，徐驰，李通化.LiF-KCl熔盐溶液的Monte Carlo法计算机模拟研究-I:径向分布函数和热力学性质[J].中国科学(B辑)，1987(1)： 21-26.

[6] Filipponi A. The radial distribution function probed by X-ray absorption spectroscopy[J].J Phys Condens Matter, 1994,6：8415-8427.

[7] 陈绍平,章社生.N体问题解析函数近似计算[J].数值计算与计算机应用,2011 (2):143-147.

2012-11-16

国家自然科学基金项目(61240048)；武汉理工大学ITS开放基金(ERCTSF2013A001)。

刘超(1990-)，男，硕士生，现主要从事应用数学方面的研究工作。

吴永红(1976-)，男，博士，讲师，现主要从事应用数学方面的教学与研究工作;E-mail:sheshengz@qq.com。

O242.2

A

1673-1409(2013)04-0004-03

[编辑] 洪云飞