韋爾萊算法是一種用於求解牛頓運動方程的數值方法,被廣泛應用於分子動力學模擬以及視頻遊戲中。韋爾萊算法的優點在於:數值穩定性比簡單的歐拉方法高很多,並保持了物理系統中的時間可逆性與相空間體積元體積守恆的性質。
Carl Størmer首次應用韋爾萊算法求解磁場中運動粒子的軌跡,因此韋爾萊算法又被稱為Størmer算法。1967年法國物理學家Loup Verlet將其應用於分子動力學計算,從此韋爾萊算法流行起來。
基本韋爾萊算法
牛頓運動方程為
![{\displaystyle a(t)={\frac {f(t)}{m}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/45c923c1a23be3388de529baac8c3e0f8b1b147a)
代入到粒子的坐標關於時間步的Taylor展式中
![{\displaystyle r(t+\Delta t)=r(t)+v(t)\Delta t+{\frac {a(t)}{2}}\Delta t^{2}+{\frac {1}{3!}}{\frac {d^{3}r}{dt^{3}}}\Delta t^{3}+{\mathcal {O}}(\Delta t^{4})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2c39aab557dbb82b168969c8d0196a57a566a448)
得
![{\displaystyle r(t+\Delta t)=r(t)+v(t)\Delta t+{\frac {f(t)}{2m}}\Delta t^{2}+{\frac {1}{3!}}{\frac {d^{3}r}{dt^{3}}}\Delta t^{3}+{\mathcal {O}}(\Delta t^{4})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/06db72e7c673fc984f1b7f5285b31cf5dd7d0f4c)
同理
![{\displaystyle r(t-\Delta t)=r(t)-v(t)\Delta t+{\frac {f(t)}{2m}}\Delta t^{2}-{\frac {1}{3!}}{\frac {d^{3}r}{dt^{3}}}\Delta t^{3}+{\mathcal {O}}(\Delta t^{4}):}](https://wikimedia.org/api/rest_v1/media/math/render/svg/62f1ef6b702d2c91200b598f8960a1093910b5ea)
兩式相加,得
![{\displaystyle r(t+\Delta t)+r(t-\Delta t)=2r(t)+{\frac {f(t)}{m}}\Delta t^{2}+{\mathcal {O}}(\Delta t^{4})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e3e27fbc338255ef56872d4730e3eb34e2f25cbc)
則
![{\displaystyle r(t+\Delta t)\simeq 2r(t)-r(t-\Delta t)+{\frac {f(t)}{m}}\Delta t^{2}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ec25f0020de773f565718af7e556fea473c70134)
新位置的計算誤差為四階,
為時間步。因為韋爾萊算法中不涉及速度,如果希望得到速度,須從軌線中推導速度表達式:
![{\displaystyle v(t)={\frac {r(t+\Delta t)-r(t-\Delta t)+{\mathcal {O}}(\Delta t^{3})}{2\Delta t}}={\frac {r(t+\Delta t)-r(t-\Delta t)}{2\Delta t}}+{\mathcal {O}}(\Delta t^{2})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/bf8b18b61494ca43882ac851fcf6611ec7ffbe18)
速度表示的韋爾萊算法
一般地,速度表示下的韋爾萊算法更為常用,它可以給出同一時間變量下的速度和位置。它實際上與基本的韋爾萊算法等價,精度相同。
首先對位置進行 Taylor 展開
![{\displaystyle r(t+\Delta t)=r(t)+v(t)\Delta t+{\frac {f(t)}{2m}}\Delta t^{2}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7d6e7d0aaa406b23d6f0df2f3638ce7ad8687b87)
![{\displaystyle r(t+2\Delta t)=r(t+\Delta t)+v(t+\Delta t)\Delta t+{\frac {f(t+\Delta t)}{2m}}\Delta t^{2}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/87605390a3c2546317522984a5caaf9d363e000f)
對兩式相減可得
![{\displaystyle r(t+2\Delta t)+r(t)=2r(t+\Delta t)+\left[v(t+\Delta t)-v(t)\right]\Delta t+{\frac {f(t+\Delta t)-f(t)}{2m}}\Delta t^{2}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/76fb51db2075f08d6cd04cf579b1bf2294a7f790)
將最初的 Verlet 公式中的
換為
,
![{\displaystyle r(t+2\Delta t)\simeq 2r(t+\Delta t)-r(t)+{\frac {f(t+\Delta t)}{m}}\Delta t^{2}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fac94b95fc3cab80e7ab41e34a1a62afe35a6e9d)
代入前式,可得
![{\displaystyle v(t+\Delta t)=v(t)+{\frac {f(t+\Delta t)+f(t)}{2m}}\Delta t}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b018f1378f929b8ba3ddf4ab5d6baac47bcd36ae)
此式即為速度表示的韋爾萊算法。實際常用的計算步驟為,
1.首先通過 Taylor 展開
計算得到位置
2.由
和系統的相互作用勢條件(如果相互作用僅依賴位置
)可以求的力場
3.由速度表示的韋爾萊公式求出新的速度
。
參考文獻
- Daan Frenkel. 第四章. "Understanding Molecular Simulation - From Algorithms to Applications" Academic Press. 2002.