簡諧振子模擬中能量的變化

董亞琴　王圣軍

摘要：本文以簡諧振子為例討論在計算機模擬中保守系統(tǒng)的能量變化。模擬和解析分析表明常微分方程通用解法使能量隨著時間指數(shù)地變化，而Verlet算法只帶來振子能量的漲落，能量不會單調(diào)變化。這些結(jié)果為理解分子動力學(xué)模擬算法的特性提供了直觀的參考。

關(guān)鍵詞：分子動力學(xué)；誤差；能量

中圖分類號：O411 文獻標(biāo)志碼：A 文章編號：1674-9324（2018）50-0193-02

在分子動力學(xué)教學(xué)中一個重要的內(nèi)容是選擇求微分方程的數(shù)值解法。盡管存在多種通用解法，能夠求解各種形式的常微分方程，但是它們在求解力學(xué)系統(tǒng)的運動方程時不能保證能量守恒。分子動力學(xué)模擬需要避免這個問題。雖然能量誤差在分子動力學(xué)中非常重要，但是各種教材只有簡單的描述。例如參考文獻[1]指出：歐拉方法的時間步長太大時，簡諧振子的相空間軌道不是橢圓，能量被粒子吸收。參考文獻[2]指出：用二階龍格庫塔方法模擬簡諧振子的運動時其相空間軌道收縮。參考文獻[3]提到：模擬簡諧振子運動時步長取得太大得到的結(jié)果偏離也大，甚至連能量都不守恒。這些說法都不容易理解，甚至不夠準(zhǔn)確。為幫助學(xué)生深入理解相關(guān)知識，本文討論簡諧振子模擬中通用解法的能量誤差。

模擬得到的能量時間關(guān)系如圖2。它不引起能量的單調(diào)變化，只帶來能量漲落，步長越小漲落的范圍越小。

本文的結(jié)果表明歐拉方法、泰勒方法、二階龍格庫塔方法導(dǎo)致諧振子能量隨著時間單調(diào)遞增，相空間軌道螺旋增大。用四階龍格庫塔方法模擬時能量隨時間單調(diào)遞減，減小步長或者改變參數(shù)不能改變能量的變化趨勢。使用Verlet算法能量只發(fā)生漲落。這些結(jié)果詳細說明了通用算法在處理牛頓運動方程時的誤差特性，為解釋分子動力學(xué)中算法的選擇提供了詳細的實例和分析，是教材的有益補充。

參考文獻：

[1]D.W.Heermann.理論物理學(xué)中的計算機模擬方法[M].北京：北京大學(xué)出版社，1996：24-25，5-6.

[2]顧昌鑫.計算物理學(xué)[M].上海：復(fù)旦大學(xué)出版社，2010：396-397.

[3]馬文淦.計算物理學(xué)[M].北京：科學(xué)出版社，2005：92-93.

[4]S.E.Koonin.計算物理學(xué)[M].北京：高等教育出版社，1993：19-20.

[5]N.J.Giordano，H.Nakanishi.計算物理[M].北京：清華大學(xué)出版社，2007：273.