劉 超，章社生, 吳永紅 (武漢理工大學(xué)理學(xué)院，湖北 武漢 430063)

多區(qū)域N體問題解析函數(shù)的近似計(jì)算

劉 超，章社生, 吳永紅 (武漢理工大學(xué)理學(xué)院，湖北 武漢 430063)

研究了粒子相互作用的多區(qū)域N體問題解析函數(shù)的近似計(jì)算，利用多重積分和delta函數(shù)，求得徑向分布函數(shù)的解析表達(dá)式，推導(dǎo)出雙球域的分析解。計(jì)算結(jié)果表明，當(dāng)粒子數(shù)為2000個(gè)時(shí)，計(jì)算誤差小于千分之一，且計(jì)算速度快。

粒子模擬；N體問題；徑向分布函數(shù)；多重積分；分析解；數(shù)值解

粒子模擬已經(jīng)成為車路協(xié)同等科學(xué)領(lǐng)域的重要研究工具[1]，它將復(fù)雜系統(tǒng)(生物大分子、星系、車路協(xié)同等)看成是大量經(jīng)典物理可描述的單個(gè)個(gè)體(原子、星體等 )的集合，然后對所有這些個(gè)體特征數(shù)據(jù)進(jìn)行分析。筆者專注于一種在粒子模擬數(shù)據(jù)分析中占有極其重要地位的雙體關(guān)聯(lián)函數(shù)查詢-空間距離直方圖(SDH)的研究[2]，并且關(guān)注Barnes2Hut算法(BHA)和快速多極子方法( FMM)[3]。SDH是一種被稱為徑向分布函數(shù)(RDF)的連續(xù)統(tǒng)計(jì)分布函數(shù)的直接估計(jì)[4-5]，與該系統(tǒng)的結(jié)構(gòu)因子(structure factor)相關(guān)聯(lián)[6]。當(dāng)研究區(qū)域?yàn)閱芜B通時(shí)，文獻(xiàn)[7]求出了粒子相互作用徑向分布函數(shù)的解析表達(dá)式，該工作能推廣到互不相連的多區(qū)域N體問題。

1 多區(qū)域粒子徑向分布

徑向分布函數(shù)(RDF)[7]的定義為：

(1)

式中，F(xiàn)h(r)是在粒子之間距離為r，厚度為h的球殼內(nèi)粒子距離數(shù)量，稱為距離和函數(shù)。設(shè)粒子坐標(biāo)為Xi，i=1,…,N，則距離和函數(shù)Fh(r)可表示為：

(2)

式中，r*=|Xi-Xj|為Xi與Xj之間的距離；當(dāng)|r*-r|≤h/2時(shí)函數(shù)δh=1，當(dāng)|r*-r|>h/2時(shí)函數(shù)δh=0。式(2)中除以2是為了消去重復(fù)累加。

用一個(gè)立方體Ω包含所有粒子。按粒子分布的稀疏情況，將立方體Ω劃分為N′個(gè)子區(qū)域Ωk(k=1,…,N′),則距離和函數(shù)Fh(r)可表示為：

(3)

式中,Fi,j為第i個(gè)區(qū)域中的粒子與第j個(gè)區(qū)域中的粒子之間的距離和函數(shù);當(dāng)j=i時(shí)表示第i個(gè)區(qū)域內(nèi)部粒子之間的距離和函數(shù),當(dāng)j≠i時(shí)表示第i個(gè)區(qū)域中的粒子與第j個(gè)區(qū)域中的粒子之間的距離和函數(shù)。

2 雙球域分析解

考慮粒子分布在2個(gè)球域O和O′，球內(nèi)分別有N和N′個(gè)粒子，對應(yīng)粒子密度分別為ρ和ρ′(粒子數(shù)除球體積)，兩球心距離為a,半徑分別為R和R′。取粒子距離和函數(shù)為:

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

(4)

圖1 雙球域?qū)?shù)示意圖

式中,r=mh,h為步長，m=1,2,3,…;F11(r)為球O內(nèi)部粒子產(chǎn)生的距離和函數(shù);F22(r)為球O′內(nèi)部粒子產(chǎn)生的距離和函數(shù);F12(r)為球O中的粒與球O′中的粒子相互作用產(chǎn)生的距離和函數(shù);F11(r)和F22(r)用文獻(xiàn)[7]中的公式計(jì)算。設(shè)點(diǎn)Y為球O′表面點(diǎn)，它到點(diǎn)O的距離為L(見圖1)，則點(diǎn)Y與球O中粒子距離和函數(shù)為：

(5)

根據(jù)空間距離直方圖(SDH)的定義[2]，將球域O劃分為M′個(gè)子區(qū)域dΩ，第i個(gè)子區(qū)域dΩi內(nèi)的粒子數(shù)為n(i)，子區(qū)域的中點(diǎn)坐標(biāo)為Yi，第i個(gè)子區(qū)域中點(diǎn)和區(qū)域外點(diǎn)Y的距離為|Yi-Y|，用n(i)|Yi-Y|近似表示第i個(gè)子區(qū)域內(nèi)所有粒子與在Y點(diǎn)粒子的距離的和。將下標(biāo)i從1取到M′，然后求和， 再根據(jù)多維積分的定義，當(dāng)dV非常小時(shí)，有:

Y(r)=0rL+R

(6)

式(6)只與距離L和R有關(guān)，與粒子坐標(biāo)Y無關(guān)。容易求得，球心在O點(diǎn)，半徑為L的球與球O′的中粒子相互距離為r的距離函數(shù)：

(7)

式中:

3 算例及誤差分析

例1設(shè)2球半徑R=R′=1，球心相距a=3, 每個(gè)球內(nèi)隨機(jī)均勻分布M個(gè)粒子,研究粒子的分布。

例1中2球之間的距離大于球的半徑，它們能放置于半徑為2.5的大球中,但用半徑為2.5的大球的粒子徑向分布解析公式計(jì)算將帶來較大的誤差。下面筆者用雙球粒子徑向分布解析公式計(jì)算。

由于半徑R=1， 則有:

(8)

且:

此時(shí)，只有2種情況:①r≤a,(L1,L2)=(a-1,r+1); ②r>a;(L1,L2)=(r-1,a+1)。

取粒子數(shù)M=500，半徑R=R′=1，兩球心距離a=3，步長h=0.1；球中粒子點(diǎn)用隨機(jī)方法產(chǎn)生。對應(yīng)的計(jì)算值繪于圖2。圖2中實(shí)線為理論值，星號為直接計(jì)算值(用2點(diǎn)之間距離公式計(jì)算)。圖2中有2個(gè)峰值，第1個(gè)峰在r=1附近，由球內(nèi)粒子相互作用產(chǎn)生。第2個(gè)峰值在r=3, 它由球與球中的粒子相互作用產(chǎn)生。由圖2可知，當(dāng)粒子僅為M=500時(shí)，理論解與直接法計(jì)算結(jié)果接近。進(jìn)一步計(jì)算表明，當(dāng)球中粒子數(shù)增加時(shí)，能減少理論計(jì)算誤差。

定義平均誤差為Er:

(9)

表1 平均誤差Er和直接法CPU計(jì)算時(shí)間t隨粒子數(shù)M的變化

圖2 2球距離和函數(shù)隨r的變化

當(dāng)R=1、R′=1、a=3、h=0.1時(shí)，平均誤差Er和CPU計(jì)算時(shí)間t如表1所示。表1中M為球內(nèi)粒子數(shù)，計(jì)算平均誤差Er乘了1000，計(jì)算時(shí)間t為直接法的CPU計(jì)算時(shí)間。由表1可知，當(dāng)粒子數(shù)M=100時(shí)，平均誤差小于2‰，粒子數(shù)增加后，平均誤差總體上減少，但發(fā)現(xiàn)有波動(dòng)。原因是粒子坐標(biāo)點(diǎn)是用隨機(jī)方法選取的，當(dāng)粒子總數(shù)M改變后，粒子坐標(biāo)都改變。由表1可知，隨著粒子數(shù)增加，CPU計(jì)算時(shí)間增加很快。

4 結(jié) 語

筆者推廣了已有的理論[7],構(gòu)造了多區(qū)域解析計(jì)算模型，得到了近似計(jì)算徑向距離的積分表達(dá)式。在粒子密度為均勻分布的條件下，得到雙球域積分解。實(shí)例計(jì)算結(jié)果表明，積分解有很高的精確度，且計(jì)算速度快。

[1]胡旻，祝大軍，劉大剛，等.粒子模擬軟件吸收邊界的研究[J].強(qiáng)激光與粒子束，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]王武，馮仰德，遲學(xué)斌.樹結(jié)構(gòu)在N體問題中的應(yīng)用[J].計(jì)算應(yīng)用研究，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法計(jì)算機(jī)模擬研究-I:徑向分布函數(shù)和熱力學(xué)性質(zhì)[J].中國科學(xué)(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體問題解析函數(shù)近似計(jì)算[J].數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用,2011 (2):143-147.

2012-11-16

國家自然科學(xué)基金項(xiàng)目(61240048)；武漢理工大學(xué)ITS開放基金(ERCTSF2013A001)。

劉超(1990-)，男，碩士生，現(xiàn)主要從事應(yīng)用數(shù)學(xué)方面的研究工作。

吳永紅(1976-)，男，博士，講師，現(xiàn)主要從事應(yīng)用數(shù)學(xué)方面的教學(xué)與研究工作;E-mail:sheshengz@qq.com。

O242.2

A

1673-1409(2013)04-0004-03

[編輯] 洪云飛