Vander Waals勢能問題不變積分保持算法研究

馬大柱

(湖北民族學院 理學院，湖北 恩施 445000)

可靠的數(shù)值方法是非線性研究的基礎(chǔ).由于傳統(tǒng)低階數(shù)值算法如四階龍格庫塔法(RK4)本身的截斷誤差較大，高階算法引入了人工耗散等因素，而相對熱門的辛算法應用有限，所以尋找合適的數(shù)值計算方法成為當前非線性工作研究的熱點問題.

自1971年Nacozy[1]提出流形改正方法以來，基于流形改正原理相繼提出了一系列的不變積分保持算法.Vander Waals勢能問題[2]存在兩個穩(wěn)定積分，為改正數(shù)值誤差，可同時穩(wěn)定兩個或者其中一個積分.基于這種考慮，采用Ma[3]設(shè)計的速度因子方法穩(wěn)定其能量積分，以及Zhong[4]構(gòu)造的最小二乘法速度因子改正方法同時穩(wěn)定兩個積分.通過數(shù)值實驗比較兩種穩(wěn)定方法優(yōu)劣，推薦適用算法，力圖避免數(shù)值誤差產(chǎn)生的混沌.

1 不變積分保持算法

對一含有m個運動積分的微分動力系統(tǒng)：φi(x)=ci(i=1,2，…,m)，這里x為狀態(tài)矢量，可分別代表坐標r和速度V，ci為運動常數(shù).理想狀態(tài)下有：

εi(x)=φ(x)-ci=0

(1)

由于數(shù)值誤差的存在，往往在實際計算中εi≠0[5].為消除數(shù)值誤差的影響，可選擇兩種方式，一種是構(gòu)造高階算法，但由于其引入人工耗散，不是理想的方案.另外一種就是改造低階算法，如各種流形改正思想，也是目前通用的方案.

1.1 速度因子方法

Ma[3]在Nacozy流形改正原理的基礎(chǔ)上，提出對速度分量可作如下形式的變換：

v*=s·v

(2)

式中v*和v分別代表速度真實值和數(shù)值解，s是標度因子.將式(2)代入到含標準積分值的表達式中，即可求出標度因子，此方法稱之為速度因子方法.此方法對系統(tǒng)只要求含有一個積分即可.當然，還有很多其它的構(gòu)造方式，此處不一一列舉.

1.2 兩個積分保持算法

若系統(tǒng)存在兩個孤立積分，除采用以上的方法穩(wěn)定一個積分外，還可以對兩個積分同時穩(wěn)定.一類仍然利用流形改正方法[6]，即將兩個積分H1和H2對坐標分量和速度分量的偏導數(shù)組成以下矩陣形式：

(3)

真實值與數(shù)值解之間滿足關(guān)系：

r*=r+Δηx，v*=v+Δηv

(4)

而Δη的計算需從下式中獲得：

Δη=-ET(EET)-1ε

(5)

另外一類就是最小二乘法速度因子改正方法[4].該方法利用到了最小二乘法原理，其構(gòu)造方式類似于速度因子改正方法.即要求速度分量都乘以一改正因子λ，有v*=λv.誤差函數(shù)φ(λ)的形式為：

(6)

這里wi是權(quán)重因子.若誤差函數(shù)最小，即要求φ′(λ)=0，可得出改正因子λ的具體表達式.

兩類不同的改正方法以上已經(jīng)詳細介紹，其不同主要體現(xiàn)在構(gòu)造方式上，尤其是速度因子方法構(gòu)造極為簡單，兩個積分同時保持算法又存在兩種不同形式，在實際運用中效果是否一樣？實際模型將給出具體回答.

2 模型與計算結(jié)果分析

以標準的Vander Waals勢能問題模型為例，其哈密頓形式為：

(7)

系統(tǒng)可積時要求B=15A.計算時為方便計，可令B=1.式(7)為該系統(tǒng)的一個積分，記為H1.此外系統(tǒng)還存在另外一個積分H2：

(8)

需注意的是式(8)是一對稱形式，將x和y，px和py交換順序，積分并無改變.獲得這兩個積分之后，將采用以上幾種穩(wěn)定方法投入實算，并進行比較.

數(shù)值比較之前，規(guī)定基本的數(shù)值積分方法為四階龍格庫塔法(RK4)，初始能量為E=1/12，龐加萊截面取(x,px)平面.首先關(guān)注無任何改正的情況，如圖1(a)所示，截面上點的分布雜亂無章，軌道在截面上的投影非常模糊，軌跡極其粗糙，在平面下部尚能看出大致規(guī)律，上方卻不能.這些都是由于RK4算法本身的誤差較大引入的，實際軌道并不如此.再看單個穩(wěn)定的情況，穩(wěn)定之前需將穩(wěn)定化程序代入到實際模型中去，獲得改正因子.圖1(b)即為能量積分穩(wěn)定情況，與圖1(a)相比，情況明顯好很多，可看出軌道穿越情況，在截面上方的分布并非雜亂無章，只是點的分布略為粗糙，仍然不能完全清晰的反映軌道動力學特征.考慮將兩個積分同時穩(wěn)定的情況，計算過程與前相同，只是穩(wěn)定化因子的求解稍顯復雜，計算結(jié)果如圖1(c)所示.從圖中可以清晰的看到軌道在截面上的分布，與真實軌道最為接近，將數(shù)值誤差的影響降低到最小.

圖2 (x,px)的相圖

圖1 (x,px)平面龐加萊截面圖

通過以上數(shù)值實驗可知，最小二乘法速度因子改正方法與速度因子方法相比，更能準確的保持系統(tǒng)的兩個積分，真實反映相空間結(jié)構(gòu)，最低限度的減少數(shù)值誤差.如前所述，穩(wěn)定兩個積分還可以采用Nacozy提出的流形改正方法.數(shù)值實驗的時候，也投入過實算，效果不太理想.究其原因，不難發(fā)現(xiàn)，該模型中第二個積分式(8)是一完全對稱形式，其對坐標兩個分量的導數(shù)完全相同，速度分量亦是如此.正是由于這樣，導致該算法不再適用.

通過龐加萊截面比較了兩種改正方法，對該模型的其它非線性性質(zhì)的研究我們推薦采用最小二乘法速度因子改正方法.在這里一并做出該模型的相圖，如圖2所示.其中圖2(a)是無改正時的情況，圖2(b)是兩個積分同時改正時的情況.通過(p1,q1)的相圖，不難發(fā)現(xiàn)單純的用RK4法計算得到的相圖與實際情況相比誤差極大.

3 結(jié)論

基于Vander Waals勢能問題模型，通過龐加萊截面分析比較了單個積分改正方法與兩個積分同時改正方法.相比較而言，速度因子方法雖然在一定程度上保持了單個積分，但是兩個積分同時改正的最小二乘法速度因子方法卻能更好的同時保持兩個積分的數(shù)值精度，更加準確的反映了相空間幾何動力學.另外，兩個積分同時改正的流形改正方法不能處理含對稱形式的積分，其應用有限.

[1]Nacozy P E.The use of integrals in numerical integtations of the n-body problem[J].Ap&SS，1971,14(11):40-51.

[2]Ganesan K, Lakshmanan M. Dynamics of atomic hydrogen in a generalized van der Waals potential[J].Phys Rev A,1990,42(10):3 490-3 497.

[3]Ma D Z. Velocity Scaling Method to Correct Individual Kepler Energies[J].New Astron,2008,13(5):216-223.

[4]Zhong S Y. Velocity Scale Factor with Least-squares Correction[J].Astrophys space sci,2009,324(1):31-40.

[5]Wu X. Comparison Among Correction Methods of Individual Kepler Energies in n-Body Simulations[J].A J,2007,133(6):2 643-2 653.

[6]Ma D Z. Velocity Corrections to Kepler Energy and Laplace Integral[J].Int J Mod Phys C,2008,19(9):1 411-1 424.