ヴェルレ統合

ヴェルレ積分（フランス語発音： [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})\,\Delta t^{2}}$
2. n = 1, 2, ... を繰り返す$${\displaystyle \mathbf {x} _{n+1}=2\mathbf {x} _{n}-\mathbf {x} _{n-1}+\mathbf {A} (\mathbf {x} _{n})\,\Delta 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}}{\Delta 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{aligned}}}$$

離散化エラー

この手法に固有の時間対称性により、離散化によって積分に導入される局所誤差のレベルが低減されます。これは、奇数次項（ここでは3次項）をすべて削除することによって実現されます。局所誤差は、反復計算に正確な値を挿入し、異なる時間方向における位置ベクトルの時刻におけるテイラー展開を計算することで定量化されます。 ここで、は位置、速度、加速度、およびジャーク（位置の時間に関する 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)}\\[1ex]\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}}}$$ ${\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 つの展開を追加すると、次の式が得られます。 テイラー展開の 1 次項と 3 次項が打ち消され、そのため、Verlet 積分器は単純なテイラー展開のみによる積分よりも 1 次正確になります。 $${\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)}.}$$

ここでの加速度は厳密解 から計算されるのに対し、反復計算では中心反復点 で計算される点に注意する必要があります。厳密解と近似値列の間の距離である全体誤差を計算する際、これら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}}$

この微分方程式にシュトラーマー法を適用すると、線形回帰関係 または が得られます 。これは、特性多項式 の根を求めることで解くことができます 。これらは です 。 線形回帰の基底解はとです。これらを厳密な解と比較するために、テイラー展開が計算されます。 この級数と指数関数の商は で始まるため 、 そこから、最初の基底解の誤差は として計算できることがわかります。 つまり、ローカル離散化誤差は 4 次ですが、微分方程式の 2 次により、グローバル誤差は 2 次となり、時間の経過とともに指数関数的に増加する定数になります。 $${\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}}\left(wh\right)^{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}}\left(wh\right)^{2}\pm wh{\sqrt {1+{\tfrac {1}{4}}\left(wh\right)^{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}}{\left(h^{4}\right)}\right)\\&=e^{wt_{n}}+{\mathcal {O}}{\left(h^{2}t_{n}e^{wt_{n}}\right)}.\end{aligned}}}$$

反復を開始する

ヴェルレ反復法の開始時（ステップ、時刻、計算）には、時刻 における位置ベクトルが既に必要であることに注意してください。一見すると、初期条件は初期時刻 においてのみ既知であるため、これは問題を引き起こす可能性があります。しかし、これらから加速度は既知であり、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 {\begin{aligned}\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)}.\end{aligned}}}$$

この場合、最初のタイムステップの誤差は のオーダーである。これは問題とはみなされない。なぜなら、多数のタイムステップにわたるシミュレーションでは、最初のタイムステップの誤差は、位置ベクトルから までの距離についても、差分商から までの距離についても のオーダーである合計誤差のうち無視できるほど小さな量に過ぎないからである。さらに、この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}.}$$

速度の計算 – Størmer-Verlet 法

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

速度ヴェルレの誤差は、基本ヴェルレと同程度であることが示されます。速度アルゴリズムは必ずしもメモリ消費量が多いわけではありません。基本ヴェルレでは位置ベクトルを2つ保持するのに対し、速度ヴェルレでは位置ベクトルと速度ベクトルをそれぞれ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}$ ${\textstyle \gamma ={\tfrac {1}{2}}}$

アルゴリズム表現

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

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

    /**
     * 「Velocity Verlet」統合を使用して pos と vel を更新します
     * @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 = new_vel ;  
        acc = new_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}}}$$

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

しかし、各粒子に複数の拘束力が作用すると問題が発生します。これを解決する一つの方法は、シミュレーション内の各点をループ処理することです。これにより、各点において、直前の拘束緩和が既に適用され、情報の伝播が高速化されます。シミュレーションでは、シミュレーションに小さな時間ステップを使用する、時間ステップごとに拘束解決ステップ数を固定する、あるいは特定の偏差に達するまで拘束を解くといった方法でこれを実装できます。

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

SuperLU ^[7]のような高度なソフトウェアは、疎行列を用いて複雑な問題を解くために存在します。行列（のクラスター）を用いるなどの特定の手法は、例えば布地を伝わる力が音波を形成せずに伝播するといった特定の問題に対処するために用いられる場合があります。 [ ^8]

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

衝突反応

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

後者の場合、ヴェルレ積分は衝突によってもたらされる速度を自動的に処理します。ただし、衝突物理法則と整合した方法で処理される保証はありません（つまり、運動量の変化が現実的であるとは保証されません）。速度項を暗黙的に変更するのではなく、衝突する物体の最終速度を明示的に制御する必要があります（前の時間ステップから記録された位置を変更することによって）。

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

参照

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

文学

1. Verlet, Loup (1967). 「古典流体に関するコンピュータ「実験」．I. レナード−ジョーンズ分子の熱力学的性質」. Physical Review . 159 (1): 98– 103. Bibcode :1967PhRv..159...98V. doi : 10.1103/PhysRev.159.98 .
2. Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). 「Section 17.4. Second-Order Conservative Equations」. Numerical Recipes: The Art of Scientific Computing (第3版). ニューヨーク: Cambridge University Press. ISBN 978-0-521-88068-8。
3. 、Wayback Machineで 2004-08-03 にアーカイブされてお​​り、Størmer メソッドの説明が記載されています。
4. Dummer, Jonathan. 「A Simple Time-Corrected Verlet Integration Method」. 2020年1月6日時点のオリジナルよりアーカイブ。2014年2月26日閲覧。
5. Swope, William C.; HC Andersen; PH Berens; KR Wilson (1982年1月1日). 「分子の物理的クラスター形成における平衡定数の計算のためのコンピュータシミュレーション法：小さな水クラスターへの応用」. The Journal of Chemical Physics . 76 (1): 648 (付録). Bibcode :1982JChPh..76..637S. doi :10.1063/1.442716.
6. エルンスト・ヘアラー;ルビッチ、クリスチャン。ワナー、ゲルハルト (2003)。 「Størmer/Verlet 法によって示される幾何学的数値積分」。アクタ ヌメリカ。12 : 399–450。書誌コード:2003AcNum..12..399H。CiteSeerX 10.1.1.7.7106。土井：10.1017/S0962492902000144。S2CID  122016794。 
7. SuperLU ユーザーズ ガイド、Wayback Machineに 2011-07-20 アーカイブ。
8. Baraff, D.; Witkin, A. (1998). 「布のシミュレーションにおける大きなステップ」(PDF) .コンピュータグラフィックス論文集. 年次会議シリーズ: 43– 54.

外部リンク

• Verlet 統合デモと Java アプレットとしてのコード
• Thomas Jakobsen 著『Advanced Character Physics』
• 分子動力学シミュレーション理論 2006年4月26日アーカイブ-ページの下部
• Verlet統合は最新のJavaScriptで実装されています – ページの下部