ヴェルレ統合

ヴェルレ積分（フランス語発音: [vɛʁˈlɛ] ）は、ニュートンの運動方程式を積分する ために使用される数値手法である。^[1]分子動力学シミュレーションやコンピュータグラフィックスで粒子の軌道を計算するために頻繁に使用される。このアルゴリズムは、1791年にジャン・バティスト・ドゥランブルによって初めて使用され、それ以来何度も再発見されており、最近では1960年代にルー・ヴェルレによって分子動力学で使用されている。また、 1909年にはPHカウエルとACCクロメリンによってハレー彗星の軌道を計算するために使用され、 1907年にはカール・ストルマーによって磁場内の電気粒子の軌道を調べるために使用された（そのためストルマー法とも呼ばれる）。^[2] ヴェルレ積分法は、単純なオイラー法に比べて大幅な計算コストの増加なしに、優れた数値安定性と、時間可逆性や位相空間上のシンプレクティック形式の保存など、物理システムで重要な他の特性を提供します。

初期条件が および である型の2階微分方程式の場合、ステップサイズが である時刻における近似数​​値解は、次の方法で得ることができます。 ${\displaystyle {\ddot {\mathbf {x} }}(t)=\mathbf {A} {\bigl (}\mathbf {x} (t){\bigr )}}$${\displaystyle \mathbf {x} (t_{0})=\mathbf {x} _{0}}$${\displaystyle {\dot {\mathbf {x} }}(t_{0})=\mathbf {v} _{0}}$${\displaystyle \mathbf {x} _{n}\approx \mathbf {x} (t_{n})}$${\displaystyle t_{n}=t_{0}+n\,\Delta t}$${\displaystyle \Delta t>0}$

1. セット、${\textstyle \mathbf {x} _{1}=\mathbf {x} _{0}+\mathbf {v} _{0}\,\Delta t+{\tfrac {1}{2}}\mathbf { A} (\mathbf {x} _{0})\,\デルタ t^{2}}$
2. n = 1, 2, ...の場合、繰り返し$${\displaystyle \mathbf {x} _{n+1}=2\mathbf {x} _{n}-\mathbf {x} _{n-1}+\mathbf {A} (\mathbf {x} _ {n})\,\デルタ t^{2}.}$$

ニュートンの保存的物理システムに対する運動方程式は

${\displaystyle {\boldsymbol {M}}{\ddot {\mathbf {x} }}(t)=F{\bigl (}\mathbf {x} (t){\bigr )}=-\nabla V{ \bigl (}\mathbf {x} (t){\bigr )},}$

または個別に

${\displaystyle m_{k}{\ddot {\mathbf {x} }}_{k}(t)=F_{k}{\bigl (}\mathbf {x} (t){\bigr )}=- \nabla _{\mathbf {x} _{k}}V\left(\mathbf {x} (t)\right),}$

どこ

• ${\displaystyle t}$時間です、
• ${\displaystyle \mathbf {x} (t)={\bigl (}\mathbf {x} _{1}(t),\ldots ,\mathbf {x} _{N}(t){\bigr )}}$物体の位置ベクトルの集合である。${\displaystyle N}$
• ${\displaystyle V}$はスカラーポテンシャル関数であり、
• ${\displaystyle F}$はポテンシャルの負の勾配であり、粒子に働く力の総和を表す。
• ${\displaystyle {\boldsymbol {M}}}$は質量行列であり、通常は各粒子の質量を持つブロックを含む対角行列です。${\displaystyle m_{k}}$

この方程式は、ポテンシャル関数のさまざまな選択に対して、相互作用する分子の運動から惑星の軌道まで、さまざまな物理システムの進化を記述するために使用できます。 ${\displaystyle V}$

質量を右側に持ってきて、多重粒子の構造を忘れる変換をすると、方程式は次のように簡略化される。

${\displaystyle {\ddot {\mathbf {x} }}(t)=\mathbf {A} {\bigl (}\mathbf {x} (t){\bigr )}}$

位置に依存する加速度を表す適切なベクトル値関数を使用します。通常、初期位置と初期速度も与えられます。 ${\displaystyle \mathbf {A} (\mathbf {x} )}$${\displaystyle \mathbf {x} (0)=\mathbf {x} _{0}}$${\displaystyle \mathbf {v} (0)={\dot {\mathbf {x} }}(0)=\mathbf {v} _{0}}$

この初期値問題を離散化して数値的に解くには、時間ステップを選択し、サンプリング点のシーケンスを考慮します。タスクは、正確な解の軌道上の 点に密接に従う点のシーケンスを構築することです。${\displaystyle \Delta t>0}$${\displaystyle t_{n}=n\,\Delta t}$${\displaystyle \mathbf {x} _{n}}$${\displaystyle \mathbf {x} (t_{n})}$

オイラー法では1 次微分方程式の 1 次導関数に前方差分近似を使用しますが、ヴェルレ積分では2 次導関数に 中心差分近似を使用すると考えられます。

${\displaystyle {\begin{aligned}{\frac {\Delta ^{2}\mathbf {x} _{n}}{\Delta t^{2}}}&={\frac {{\frac {\ mathbf {x} _{n+1}-\mathbf {x} _{n}}{\Delta t}}-{\frac {\mathbf {x} _{n}-\mathbf {x} _{n-1}}{\Delta t}}}{\Delta t}}\\[6pt]&={\frac {\mathbf {x} _{n+ 1}-2\mathbf {x} _{n}+\mathbf {x} _{n-1}}{\デルタ t^{2}}}=\mathbf {a} _{n}=\mathbf {A} (\mathbf {x} _{n}).\end{aligned}}}$

シュテルマー法^[3]で用いられる形式のベルレ積分は、この式を用いて、速度を使わずに前の2つの位置ベクトルから次の位置ベクトルを求める。

${\displaystyle {\begin{aligned}\mathbf {x} _{n+1}&=2\mathbf {x} _{n}-\mathbf {x} _{n-1}+\mathbf {a} _{n}\,\Delta t^{2},\\[6pt]\mathbf {a} _{n}&=\mathbf {A} (\mathbf {x} _{n}).\end{整列}}}$

この方法に固有の時間対称性により、離散化によって積分に導入される局所誤差のレベルが、奇数次項（ここでは 3 次項）をすべて削除することによって低減されます。局所誤差は、反復に正確な値を挿入し、異なる時間方向の位置ベクトルの時点でのテイラー展開を計算することによって定量化されます。 ${\displaystyle \Delta t}$${\displaystyle \mathbf {x} (t_{n-1}),\mathbf {x} (t_{n}),\mathbf {x} (t_{n+1})}$${\displaystyle t=t_{n}}$${\displaystyle \mathbf {x} (t\pm \Delta t)}$

${\displaystyle {\begin{aligned}\mathbf {x} (t+\Delta t)&=\mathbf {x} (t)+\mathbf {v} (t)\Delta t+{\frac {\mathbf {a} (t)\Delta t^{2}}{2}}+{\frac {\mathbf {b} (t)\Delta t^{3}}{6}}+{\mathcal {O}}\left(\Delta t^{4}\right)\\\mathbf {x} (t-\Delta t)&=\mathbf {x} (t)-\mathbf {v} (t)\Delta t+{\frac {\mathbf {a} (t)\Delta t^{2}}{2}}-{\frac {\mathbf {b} (t)\Delta t^{3}}{6}}+{\mathcal {O}}\left(\Delta t^{4}\right),\end{aligned}}}$

ここで、は位置、速度、加速度、ジャーク（時間に対する位置の 3 次微分）です。 ${\displaystyle \mathbf {x} }$${\displaystyle \mathbf {v} ={\dot {\mathbf {x} }}}$${\displaystyle \mathbf {a} ={\ddot {\mathbf {x} }}}$${\displaystyle \mathbf {b} ={\dot {\mathbf {a} }}={\overset {\dots }{\mathbf {x} }}}$

これら2つの展開を加えると

${\displaystyle \mathbf {x} (t+\Delta t)=2\mathbf {x} (t)-\mathbf {x} (t-\Delta t)+\mathbf {a} (t)\Delta t^{2}+{\mathcal {O}}\left(\Delta t^{4}\right).}$

テイラー展開の 1 次項と 3 次項が打ち消され、そのため、Verlet 積分器は単純なテイラー展開のみによる積分よりも 1 次正確になることがわかります。

ここでの加速度は厳密解 から計算されるのに対し、反復では中心反復点 で計算される点に注意する必要があります。厳密解と近似シーケンス間の距離であるグローバル誤差を計算する際、これら 2 つの項は正確には打ち消されず、グローバル誤差の順序に影響を及ぼします。 ${\displaystyle \mathbf {a} (t)=\mathbf {A} {\bigl (}\mathbf {x} (t){\bigr )}}$${\displaystyle \mathbf {a} _{n}=\mathbf {A} (\mathbf {x} _{n})}$

局所誤差と全体誤差の関係を理解するには、正確な解と近似解を明示的な式で表現できる簡単な例を調べると役立ちます。このタスクの標準的な例は指数関数です。

定数 を持つ線形微分方程式を考えます。その正確な基底解はおよびです。 ${\displaystyle {\ddot {x}}(t)=w^{2}x(t)}$${\displaystyle w}$${\displaystyle e^{wt}}$${\displaystyle e^{-wt}}$

この微分方程式にストーマー法を適用すると、線形回帰関係が得られる。

${\displaystyle x_{n+1}-2x_{n}+x_{n-1}=h^{2}w^{2}x_{n},}$

または

${\displaystyle x_{n+1}-2\left(1+{\tfrac {1}{2}}(wh)^{2}\right)x_{n}+x_{n-1}=0.}$

これは、その特性多項式の根を見つけることで解くことができます 。これらは、 ${\displaystyle q^{2}-2\left(1+{\tfrac {1}{2}}(wh)^{2}\right)q+1=0}$

${\displaystyle q_{\pm }=1+{\tfrac {1}{2}}(wh)^{2}\pm wh{\sqrt {1+{\tfrac {1}{4}}(wh)^{2}}}.}$

線形回帰の基本解はとです。これらを厳密解と比較するために、テイラー展開を計算します。 ${\displaystyle x_{n}=q_{+}^{n}}$${\displaystyle x_{n}=q_{-}^{n}}$

${\displaystyle {\begin{aligned}q_{+}&=1+{\tfrac {1}{2}}(wh)^{2}+wh\left(1+{\tfrac {1}{8}}(wh)^{2}-{\tfrac {3}{128}}(wh)^{4}+{\mathcal {O}}\left(h^{6}\right)\right)\\&=1+(wh)+{\tfrac {1}{2}}(wh)^{2}+{\tfrac {1}{8}}(wh)^{3}-{\tfrac {3}{128}}(wh)^{5}+{\mathcal {O}}\left(h^{7}\right).\end{aligned}}}$

この級数と指数級数の商は から始まるので、 ${\displaystyle e^{wh}}$${\displaystyle 1-{\tfrac {1}{24}}(wh)^{3}+{\mathcal {O}}\left(h^{5}\right)}$

${\displaystyle {\begin{aligned}q_{+}&=\left(1-{\tfrac {1}{24}}(wh)^{3}+{\mathcal {O}}\left(h^{5}\right)\right)e^{wh}\\&=e^{-{\frac {1}{24}}(wh)^{3}+{\mathcal {O}}\left(h^{5}\right)}\,e^{wh}.\end{aligned}}}$

そこから、最初の基底解に対する誤差は次のように計算できることがわかる。

${\displaystyle {\begin{aligned}x_{n}=q_{+}^{n}&=e^{-{\frac {1}{24}}(wh)^{2}\,wt_{n}+{\mathcal {O}}\left(h^{4}\right)}\,e^{wt_{n}}\\&=e^{wt_{n}}\left(1-{\tfrac {1}{24}}(wh)^{2}\,wt_{n}+{\mathcal {O}}(h^{4})\right)\\&=e^{wt_{n}}+{\mathcal {O}}\left(h^{2}t_{n}e^{wt_{n}}\right).\end{aligned}}}$

つまり、ローカル離散化誤差は 4 次ですが、微分方程式の 2 次により、グローバル誤差は 2 次となり、時間の経過とともに指数関数的に増加する定数を持ちます。

計算中のステップ、時間での Verlet 反復の開始時に、時間 の位置ベクトルがすでに必要であることに注意してください。一見すると、初期条件は初期時間 でのみわかっているため、これは問題を引き起こす可能性があります。ただし、これらから加速度はわかっており、 2 次 テイラー多項式を使用して最初の時間ステップでの位置の適切な近似値を取得できます。${\displaystyle n=1}$${\displaystyle t=t_{1}=\Delta t}$${\displaystyle \mathbf {x} _{2}}$${\displaystyle \mathbf {x} _{1}}$${\displaystyle t=t_{1}}$${\displaystyle t_{0}=0}$${\displaystyle \mathbf {a} _{0}=\mathbf {A} (\mathbf {x} _{0})}$

${\displaystyle \mathbf {x} _{1}=\mathbf {x} _{0}+\mathbf {v} _{0}\Delta t+{\tfrac {1}{2}}\mathbf {a} _{0}\Delta t^{2}\approx \mathbf {x} (\Delta t)+{\mathcal {O}}\left(\Delta t^{3}\right).}$

この場合、最初のタイム ステップの誤差は のオーダーになります。多数のタイム ステップにわたるシミュレーションでは、最初のタイム ステップの誤差は、位置ベクトルから までの距離についても、差の商から までの距離についても、 の時点での のオーダーである合計誤差のうち無視できるほど小さい量にすぎないため、これは問題とは見なされません。さらに、この 2 次グローバル誤差を取得するには、初期誤差が少なくとも 3 次である必要があります。 ${\displaystyle {\mathcal {O}}\left(\Delta t^{3}\right)}$${\displaystyle t_{n}}$${\displaystyle {\mathcal {O}}\left(e^{Lt_{n}}\Delta t^{2}\right)}$${\displaystyle \mathbf {x} _{n}}$${\displaystyle \mathbf {x} (t_{n})}$${\displaystyle {\tfrac {\mathbf {x} _{n+1}-\mathbf {x} _{n}}{\Delta t}}}$${\displaystyle {\tfrac {\mathbf {x} (t_{n+1})-\mathbf {x} (t_{n})}{\Delta t}}}$

ストルマー・ヴェルレ法の欠点は、時間ステップ（）が変化すると、微分方程式の解を近似できないことである。これは、式^{[4]を使用して修正することができる。}${\displaystyle \Delta t}$

${\displaystyle \mathbf {x} _{i+1}=\mathbf {x} _{i}+\left(\mathbf {x} _{i}-\mathbf {x} _{i-1}\right){\frac {\Delta t_{i}}{\Delta t_{i-1}}}+\mathbf {a} _{i}\Delta t_{i}^{2}.}$

より正確な導出は、時刻と時刻についてテイラー級数（2次）を使い、${\displaystyle t_{i}}$${\displaystyle t_{i+1}=t_{i}+\Delta t_{i}}$${\displaystyle t_{i-1}=t_{i}-\Delta t_{i-1}}$${\displaystyle \mathbf {v} _{i}}$

${\displaystyle {\frac {\mathbf {x} _{i+1}-\mathbf {x} _{i}}{\Delta t_{i}}}+{\frac {\mathbf {x} _{i-1}-\mathbf {x} _{i}}{\Delta t_{i-1}}}=\mathbf {a} _{i}\,{\frac {\Delta t_{i}+\Delta t_{i-1}}{2}},}$

反復式は次のようになる。

${\displaystyle \mathbf {x} _{i+1}=\mathbf {x} _{i}+(\mathbf {x} _{i}-\mathbf {x} _{i-1}){\frac {\Delta t_{i}}{\Delta t_{i-1}}}+\mathbf {a} _{i}\,{\frac {\Delta t_{i}+\Delta t_{i-1}}{2}}\,\Delta t_{i}.}$

速度は基本的なシュテルマー方程式では明示的に与えられていませんが、運動エネルギーなどの特定の物理量の計算にはしばしば必要です。これは分子動力学シミュレーションで技術的な課題を生み出す可能性があります。なぜなら、時間 における位置がわかるまでは、システムの運動エネルギーと瞬間温度を計算できないからです。この欠陥は、速度ベルレアルゴリズムを使用するか、位置項と平均値定理を使用して速度を推定することで対処できます。 ${\displaystyle t}$${\displaystyle t+\Delta t}$

${\displaystyle \mathbf {v} (t)={\frac {\mathbf {x} (t+\Delta t)-\mathbf {x} (t-\Delta t)}{2\Delta t}}+{\mathcal {O}}\left(\Delta t^{2}\right).}$

この速度項は位置項より 1 ステップ遅れていることに注意してください。これは、 ではなく、時間 における速度に対するものであるため、は の 2 次近似であることを意味します。同じ議論ですが、時間ステップを半分にすると、で の2 次近似になります。 ${\displaystyle t}$${\displaystyle t+\Delta t}$${\displaystyle \mathbf {v} _{n}={\tfrac {\mathbf {x} _{n+1}-\mathbf {x} _{n-1}}{2\Delta t}}}$${\displaystyle \mathbf {v} (t_{n})}$${\displaystyle \mathbf {v} _{n+{\frac {1}{2}}}={\tfrac {\mathbf {x} _{n+1}-\mathbf {x} _{n}}{\Delta t}}}$${\displaystyle \mathbf {v} \left(t_{n+{\frac {1}{2}}}\right)}$${\displaystyle t_{n+{\frac {1}{2}}}=t_{n}+{\tfrac {1}{2}}\Delta t}$

精度を犠牲にして、 時間における速度を近似するために間隔を短くすることができます。${\displaystyle t+\Delta t}$

${\displaystyle \mathbf {v} (t+\Delta t)={\frac {\mathbf {x} (t+\Delta t)-\mathbf {x} (t)}{\Delta t}}+{\mathcal {O}}(\Delta t).}$

関連し、より一般的に使用されているアルゴリズムは速度ヴェルレアルゴリズム^{[5]であり、これは}リープフロッグ法に似ていますが、速度と位置が時間変数の同じ値で計算される点が異なります（リープフロッグ法では、名前が示すように、同じ値で計算されません）。これは同様のアプローチを使用しますが、速度を明示的に組み込んで、基本的なヴェルレアルゴリズムの最初の時間ステップの問題を解決します。

${\displaystyle {\begin{aligned}\mathbf {x} (t+\Delta t)&=\mathbf {x} (t)+\mathbf {v} (t)\,\Delta t+{\tfrac {1}{2}}\,\mathbf {a} (t)\Delta t^{2},\\[6pt]\mathbf {v} (t+\Delta t)&=\mathbf {v} (t)+{\frac {\mathbf {a} (t)+\mathbf {a} (t+\Delta t)}{2}}\Delta t.\end{aligned}}}$

速度 Verlet の誤差は、基本 Verlet の誤差と同じ程度であることが示されています。速度アルゴリズムは必ずしもメモリを多く消費するわけではありません。基本 Verlet では 2 つの位置ベクトルを追跡しますが、速度 Verlet では 1 つの位置ベクトルと 1 つの速度ベクトルを追跡します。このアルゴリズムの標準的な実装スキームは次のとおりです。

1. 計算します。${\displaystyle \mathbf {v} \left(t+{\tfrac {1}{2}}\,\Delta t\right)=\mathbf {v} (t)+{\tfrac {1}{2}}\,\mathbf {a} (t)\,\Delta t}$
2. 計算します。${\displaystyle \mathbf {x} (t+\Delta t)=\mathbf {x} (t)+\mathbf {v} \left(t+{\tfrac {1}{2}}\,\Delta t\right)\,\Delta t}$
3. 相互作用ポテンシャルから を使用して導出します。${\displaystyle \mathbf {a} (t+\Delta t)}$${\displaystyle \mathbf {x} (t+\Delta t)}$
4. 計算します。${\displaystyle \mathbf {v} (t+\Delta t)=\mathbf {v} \left(t+{\tfrac {1}{2}}\,\Delta t\right)+{\tfrac {1}{2}}\,\mathbf {a} (t+\Delta t)\Delta t}$

このアルゴリズムは可変時間ステップでも機能し、リープフロッグ法の積分の「キック・ドリフト・キック」形式と同一です。

半音速を排除すると、このアルゴリズムは次のように短縮される。

1. 計算します。${\displaystyle \mathbf {x} (t+\Delta t)=\mathbf {x} (t)+\mathbf {v} (t)\,\Delta t+{\tfrac {1}{2}}\,\mathbf {a} (t)\,\Delta t^{2}}$
2. 相互作用ポテンシャルから を使用して導出します。${\displaystyle \mathbf {a} (t+\Delta t)}$${\displaystyle \mathbf {x} (t+\Delta t)}$
3. 計算します。${\displaystyle \mathbf {v} (t+\Delta t)=\mathbf {v} (t)+{\tfrac {1}{2}}\,{\bigl (}\mathbf {a} (t)+\mathbf {a} (t+\Delta t){\bigr )}\Delta t}$

ただし、このアルゴリズムでは、加速度は位置にのみ依存し、速度には依存しないと想定していることに注意してください。 ${\displaystyle \mathbf {a} (t+\Delta t)}$${\displaystyle \mathbf {x} (t+\Delta t)}$${\displaystyle \mathbf {v} (t+\Delta t)}$

速度ヴェルレ法、および同様にリープフロッグ法の長期的な結果は、半暗黙的オイラー法よりも 1桁優れていることに留意してください。アルゴリズムは、速度の半分の時間ステップのシフトまではほとんど同じです。これは、上記のループを回転させてステップ 3 から開始し、ステップ 2 と 4 を組み合わせることでステップ 1 の加速項を除去できることに気付くことで証明できます。唯一の違いは、速度ヴェルレ法の中間点の速度が半暗黙的オイラー法の最終速度と見なされることです。

すべてのオイラー法のグローバル誤差は 1 次ですが、この方法のグローバル誤差は中点法と同様に 2 次です。さらに、加速度が実際に保存力学系またはハミルトン系の力から生じる場合、近似のエネルギーは厳密に解かれた系の定数エネルギーの周りを本質的に振動し、グローバル誤差は半明示的オイラー法では 1 次、ヴェルレ・リープフロッグ法では 2 次になります。同じことが、線型運動量や角運動量など、常に保存されるか、シンプレクティック積分器でほぼ保存される系の他のすべての保存量にも当てはまります。^[6]

速度ヴェルレ法は、およびを使用したニューマークベータ法の特殊なケースです。 ${\displaystyle \beta =0}$${\displaystyle \gamma ={\tfrac {1}{2}}}$

速度 Verlet は3D アプリケーションで一般的に役立つアルゴリズムである ため、C++ で記述されたソリューションは以下のようになります。このタイプの位置統合により、通常のオイラー法と比較して、3D シミュレーションやゲームの精度が大幅に向上します。

構造体本体 
{
    ベクトル3d位置{ 0.0,0.0,0.0 } ;​​      
    Vec3d vel { 2.0 , 0.0 , 0.0 }; // x軸に沿って2m/s       
    Vec3d acc { 0.0 , 0.0 , 0.0 }; // 最初は加速なし       
    倍質量= 1.0 ; // 1kg    

    /**
     * 「Velocity Verlet」統合を使用して位置と速度を更新します
     * @param dt DeltaTime / タイムステップ [例: 0.01]
     */
    void更新( double dt )  
    {
        Vec3d new_pos = pos + vel * dt + acc * ( dt * dt * 0.5 );       
        Vec3d new_acc = apply_forces ();   
        Vec3d new_vel = vel + ( acc + new_acc ) * ( dt * 0.5 );     
        位置=新しい位置;  
        vel =新しいvel ;  
        acc =新しいacc ;  
    }

    /**
     * オブジェクトに速度を適用するには、代わりに必要な力ベクトルを計算します。
     * ここで蓄積された力を適用します。
     */
    Vec3d apply_forces ()定数  
    {
        Vec3d new_acc = Vec3d { 0.0 , 0.0 , -9.81 }; // Z軸の下方向に9.81 m/s²       
        // ここで他の力を適用します...
        // 注意: Velocity Verlet は加速度が位置に依存すると想定しているため、`vel` に依存しないようにしてください。
        new_accを返します。 
    }
};

ベルレ法のグローバルな打ち切り誤差は、位置と速度の両方で です。これは、位置のローカル誤差が上記のとおりであるという事実とは対照的です。この違いは、すべての反復にわたってローカルな打ち切り誤差が蓄積されるためです。 ${\displaystyle {\mathcal {O}}\left(\Delta t^{2}\right)}$${\displaystyle {\mathcal {O}}\left(\Delta t^{4}\right)}$

グローバル エラーは、次の点に注意して導き出すことができます。

${\displaystyle \operatorname {error} {\bigl (}x(t_{0}+\Delta t){\bigr )}={\mathcal {O}}\left(\Delta t^{4}\right)}$

そして

${\displaystyle x(t_{0}+2\Delta t)=2x(t_{0}+\Delta t)-x(t_{0})+\Delta t^{2}{\ddot {x}}(t_{0}+\Delta t)+{\mathcal {O}}\left(\Delta t^{4}\right).}$

したがって

${\displaystyle \operatorname {error} {\bigl (}x(t_{0}+2\Delta t){\bigr )}=2\cdot \operatorname {error} {\bigl (}x(t_{0}+\Delta t){\bigr )}+{\mathcal {O}}\left(\Delta t^{4}\right)=3\,{\mathcal {O}}\left(\Delta t^{4}\right).}$

同様に:

${\displaystyle {\begin{aligned}\operatorname {error} {\bigl (}x(t_{0}+3\Delta t){\bigl )}&=6\,{\mathcal {O}}\left(\Delta t^{4}\right),\\[6px]\operatorname {error} {\bigl (}x(t_{0}+4\Delta t){\bigl )}&=10\,{\mathcal {O}}\left(\Delta t^{4}\right),\\[6px]\operatorname {error} {\bigl (}x(t_{0}+5\Delta t){\bigl )}&=15\,{\mathcal {O}}\left(\Delta t^{4}\right),\end{aligned}}}$

これは次のように一般化できます (帰納法で示すこともできますが、ここでは証明なしで示します)。

${\displaystyle \operatorname {error} {\bigl (}x(t_{0}+n\Delta t){\bigr )}={\frac {n(n+1)}{2}}\,{\mathcal {O}}\left(\Delta t^{4}\right).}$

との間の位置のグローバル誤差を考慮すると、は明らかである^{[引用が必要]}${\displaystyle x(t)}$${\displaystyle x(t+T)}$${\displaystyle T=n\Delta t}$

${\displaystyle \operatorname {error} {\bigl (}x(t_{0}+T){\bigr )}=\left({\frac {T^{2}}{2\Delta t^{2}}}+{\frac {T}{2\Delta t}}\right){\mathcal {O}}\left(\Delta t^{4}\right),}$

したがって、一定時間間隔における全体的な（累積的な）誤差は次のように表される。

${\displaystyle \operatorname {error} {\bigr (}x(t_{0}+T){\bigl )}={\mathcal {O}}\left(\Delta t^{2}\right).}$

速度はVerlet積分器の位置から非累積的に決定されるため、速度のグローバル誤差も です。 ${\displaystyle {\mathcal {O}}\left(\Delta t^{2}\right)}$

分子動力学シミュレーションでは、通常、グローバル エラーはローカル エラーよりもはるかに重要であるため、Verlet 積分器は 2 次積分器として知られています。

制約のある複数の粒子のシステムは、オイラー法よりもベルレ積分で解く方が簡単です。ポイント間の制約は、たとえば、特定の距離にそれらを制約するポテンシャルや引力です。これらは、粒子を接続するバネとしてモデル化できます。無限の剛性のバネを使用すると、モデルはベルレ アルゴリズムで解くことができます。

1次元では、制約距離が である場合、制約のない位置と時刻における点の実際の位置の関係は、アルゴリズムによって見つけることができます。 ${\displaystyle {\tilde {x}}_{i}^{(t)}}$${\displaystyle x_{i}^{(t)}}$${\displaystyle i}$${\displaystyle t}$${\displaystyle r}$

${\displaystyle {\begin{aligned}d_{1}&=x_{2}^{(t)}-x_{1}^{(t)},\\[6px]d_{2}&=\|d_{1}\|,\\[6px]d_{3}&={\frac {d_{2}-r}{d_{2}}},\\[6px]x_{1}^{(t+\Delta t)}&={\tilde {x}}_{1}^{(t+\Delta t)}+{\tfrac {1}{2}}d_{1}d_{3},\\[6px]x_{2}^{(t+\Delta t)}&={\tilde {x}}_{2}^{(t+\Delta t)}-{\tfrac {1}{2}}d_{1}d_{3}.\end{aligned}}}$

ベルレ積分は、速度を使用して問題を解くのではなく、力と位置を直接関連付けるので便利です。

ただし、複数の拘束力が各粒子に作用すると問題が発生します。これを解決する 1 つの方法は、シミュレーション内のすべてのポイントをループすることです。これにより、すべてのポイントで、最後の拘束緩和がすでに使用されて、情報の拡散が高速化されます。シミュレーションでは、シミュレーションに小さな時間ステップを使用するか、時間ステップごとに固定数の拘束解決ステップを使用するか、特定の偏差に一致するまで拘束を解決することでこれを実装できます。

制約を局所的に一次近似する場合、これはガウス・ザイデル法と同じです。小さな行列の場合、 LU 分解の方が高速であることがわかっています。大規模なシステムはクラスターに分割できます (たとえば、各ラグドール = クラスター)。クラスター内では LU 法が使用され、クラスター間ではガウス・ザイデル法が使用されます。行列コードは再利用できます。位置に対する力の依存性は局所的に一次近似でき、ベルレ積分をより暗黙的にすることができます。

^{SuperLU [7]}のような洗練されたソフトウェアは、疎行列を使用して複雑な問題を解決します。行列（のクラスター）を使用するなどの特定の技術は、音波を形成せずに布のシートを介して伝播する力などの特定の問題に対処するために使用される場合があります。^[8]

ホロノミック制約を解決する別の方法は、制約アルゴリズムを使用することです。

衝突に対処する方法の 1 つは、ペナルティ ベースのシステムを使用することです。これは基本的に、接触時にポイントに一定の力を適用します。これの問題は、与える力を選択するのが非常に難しいことです。力が強すぎると、オブジェクトは不安定になり、力が弱すぎると、オブジェクトが互いに貫通してしまいます。もう 1 つの方法は、投影衝突反応を使用することです。これは、問題のあるポイントを取得し、それを他のオブジェクトからできるだけ短い距離に移動させようとします。

後者の場合、Verlet 積分は衝突によってもたらされる速度を自動的に処理しますが、衝突の物理法則と一致する方法で処理されることは保証されません(つまり、運動量の変化が現実的であるとは保証されません)。速度項を暗黙的に変更する代わりに、衝突するオブジェクトの最終速度を明示的に制御する必要があります (前の時間ステップから記録された位置を変更することによって)。

新しい速度を決定する最も簡単な方法は、完全弾性衝突と非弾性衝突の2 つです。より制御性の高い、もう少し複雑な戦略としては、反発係数の使用があります。

• クーラント・フリードリヒス・レヴィ条件
• エネルギードリフト
• シンプレクティック積分器
• リープフロッグ統合
• ビーマンのアルゴリズム

• Verlet 統合デモと Java アプレットとしてのコード
• Thomas Jakobsen 著「Advanced Character Physics」
• 分子動力学シミュレーションの理論 – ページ下部
• Verlet 統合は最新の JavaScript で実装されています – ページの下部