ベレの方法のソースを表示

{{翻訳直後|[[:en:Special:Permalink/1024685434|Verlet integration]]|date=2023年5月}}

'''ベレの方法'''（ベレのほうほう、{{lang-en-short|Verlet algorithm}}）は、[[ニュートンの運動方程式]]を[[常微分方程式の数値解法|数値積分]]する手法の一つ<ref name="Verlet1967">{{cite journal|last=Verlet|first=Loup|year=1967|title=Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard−Jones Molecules|journal=Physical Review|volume=159|issue=1|pages=98–103|bibcode=1967PhRv..159...98V|doi=10.1103/PhysRev.159.98|doi-access=free}}</ref>。'''ベレのアルゴリズム'''、'''ベレ法'''、'''ベルレ積分法'''({{Lang|en|Verlet integration}})、'''ベルレの方法'''などの呼び方もある。[[分子動力学法]]における粒子の軌跡（トラジェクトリ）のシミュレーションや[[コンピュータグラフィックス]]に頻繁に用いられる。1791年に[[ジャン＝バティスト・ジョゼフ・ドランブル]]が用いたのが最初で、その後も何回も再発見されているが、1960年代に[[:en:Loup_Verlet|Loup Verlet]]が分子動力学法に用いてから広く使われるようになった。1909年に[[ハレー彗星]]の[[軌道 (力学)|軌道]]を計算するために[[フィリップ・コーウェル]]と[[アンドリュー・クロンメリン]]により用いられたり、1907年に磁場中の[[荷電粒子]]のトラジェクトリを研究するために[[:en:Carl_Størmer|Carl Størmer]]により用いられたりしている（そのため、'''Störmerの方法'''という異称もある）。<ref>{{Cite book |last1=Press |first1=W. H. |last2=Teukolsky |first2=S. A. |last3=Vetterling |first3=W. T. |last4=Flannery |first4=B. P. |year=2007 |title=Numerical Recipes: The Art of Scientific Computing |edition=3rd |publisher=Cambridge University Press |location=New York |isbn=978-0-521-88068-8 |chapter=Section 17.4. Second-Order Conservative Equations |chapter-url=http://apps.nrbook.com/empanel/index.html#pg=928}}</ref>この手法の特徴として、[[数値的安定性]]が高く、[[T対称性|時間反転対称性]]を持つことや[[位相空間 (物理学)|位相空間]]上の[[シンプレクティック数値積分法|シンプレクティック形式を保存]]するなど[[系 (自然科学)|物理系]]において重要な性質を持つうえ、[[オイラー法]]と比べてもほとんど計算コストが増えないことが挙げられる。

== 基本 ==
<math>\ddot{\vec x}(t) = \vec A\big(\vec x(t)\big)</math>の形の2階[[微分方程式]]を初期条件<math>\vec x(t_0) = \vec x_0</math>および<math>\dot{\vec x}(t_0) = \vec v_0</math>のもとで、<math>\Delta t > 0</math>をステップサイズとして、<math>t_n = t_0 + n\,\Delta t</math>における数値的近似解<math>\vec x_n \approx \vec x(t_n)</math>は、下のようなアルゴリズムで求めることができる。

* <math>\vec x_1 = \vec x_0 + \vec v_0\,\Delta t + \frac12 \vec A(\vec x_0)\,\Delta t^2</math> と置く。
* ''n'' = 1、2、...について次を繰り返す。
*: <math>
\vec x_{n+1} = 2 \vec x_n - \vec x_{n-1} + \vec A(\vec x_n)\,\Delta t^2.
</math>

=== 運動方程式 ===
[[保存系]]におけるニュートンの運動方程式は

: <math>M \ddot{\vec x}(t) = F\big(\vec x(t)\big) = -\nabla V\big(\vec x(t)\big)</math>

もしくは粒子毎に

: <math>m_k \ddot{\vec x}_k(t) = F_k\big(\vec x(t)\big) = -\nabla_{\vec x_k} V\big(\vec x(t)\big),</math>

と表わされる。ここで、

: {{Mvar|t}}は時刻、
: <math>\vec x(t) = \big(\vec x_1(t), \ldots, \vec x_N(t)\big)</math> は{{Mvar|N}}個の物体の位置ベクトル、
: {{Mvar|V}}はスカラーポテンシャル関数、
: {{mvar|F}}は負のポテンシャル勾配、すなわち粒子にはたらく力の集合、
: {{mvar|M}}[[質量行列]]で、各粒子の質量<math>m_k</math>を対角要素にもつ行列である。

この等式は、[[分子動力学法|相互作用する分子群]]や[[多体問題|惑星の軌道]]などさまざまな物理系が、さまざまなポテンシャル関数{{Mvar|V}}を決めたときどのように時間発展するかを記述するために用いることができる。

右辺に質量を移し、粒子系の構造を忘れることにすると、上の式は以下のように単純化することができる。

: <math>\ddot{\vec x}(t) = A\big(\vec x(t)\big)</math>

位置に依存する加速度を表す適切なベクトル値関数''A''を使用します。通常、も与えられます。ここで、ベクトル値関数{{mvar|A}}は位置依存加速度をあらわす。通常は、初期位置<math>\vec x(0) = \vec x_0</math>および初速度<math>\vec v(0) = \dot{\vec x}(0) = \vec v_0</math>も所与である。

=== ベレ積分（速度なし） ===
この[[初期値問題]]を離散化し、数値的に解くため、タイムステップ<math>\Delta t > 0</math>を選び、サンプル点列<math>t_n = n\,\Delta t</math>を考える。このとき、問題は厳密解<math>\vec x(t_n)</math>をよく近似する点群列<math>\vec x_n</math>をもとめることに帰着する。

[[オイラー法]]が1階微分方程式中の1次微分を[[有限差分|前進差分]]近似するのに対し、ベレ積分は2次微分を中心差分近似すると見ることができる。

: <math>
\frac{\Delta^2\vec x_n}{\Delta t^2}
= \frac{\frac{\vec x_{n+1} - \vec x_n}{\Delta t} - \frac{\vec x_n - \vec x_{n-1}}{\Delta t}}{\Delta t}
= \frac{\vec x_{n+1} - 2 \vec x_n + \vec x_{n-1}}{\Delta t^2} = \vec a_n = \vec A(\vec x_n)
</math>

「Störmerの方法」<ref>[http://www.fisica.uniud.it/~ercolessi/md/md/node21.html webpage] {{Webarchive|url=https://web.archive.org/web/20040803212552/http://www.fisica.uniud.it/~ercolessi/md/md/node21.html|date=2004-08-03}} with a description of the Störmer method.</ref>で用いられる形式の「ベレ積分」は、 この式を用い、速度を使わずに以前の2つの位置から次の座標を与える。

: <math>
\vec x_{n+1} = 2 \vec x_n - \vec x_{n-1} + \vec a_n\,\Delta t^2,
\qquad \vec a_n = \vec A(\vec x_n)
</math>

=== 離散化誤差 ===
このアルゴリズムは本質的に時間反転対称性を持ち、離散化の際に奇数次の項が消え、<math>\Delta t</math>について3次の項がなくなるため、局所誤差の大きさを低減することができる。局所誤差は厳密解<math>\vec x(t_{n-1}), \vec x(t_n), \vec x(t_{n+1})</math>を代入し、<math>t = t_n</math>における[[テイラー展開]]から<math>\vec{x}(t \pm \Delta t)</math>をそれぞれ計算することにより定量できる。

: <math>\begin{align}
\vec{x}(t + \Delta t)
&= \vec{x}(t) + \vec{v}(t)\Delta t + \frac{\vec{a}(t) \Delta t^2}{2}
+ \frac{\vec{b}(t) \Delta t^3}{6} + \mathcal{O}(\Delta t^4)\\
\vec{x}(t - \Delta t)
&= \vec{x}(t) - \vec{v}(t)\Delta t + \frac{\vec{a}(t) \Delta t^2}{2}
- \frac{\vec{b}(t) \Delta t^3}{6} + \mathcal{O}(\Delta t^4)
\end{align}</math>

ここで、<math>\vec{x}</math>は位置、<math>\vec{v} = \dot{\vec x}</math>は速度、<math>\vec{a} = \ddot{\vec x}</math>は加速度、<math>\vec{b}</math>は[[躍度]]である。

これら2つの級数を足し合わせると、以下の式を得る。

: <math>\vec{x}(t + \Delta t) = 2\vec{x}(t) - \vec{x}(t - \Delta t) + \vec{a}(t) \Delta t^2 + \mathcal{O}(\Delta t^4).</math>

この式から、テイラー展開の1次および3次の項がうち消しあい、ベレ積分の精度が単なるテイラー展開よりも1次高くなっていることがわかる。

ここで、加速度<math>\vec a(t) = A\big(\vec x(t)\big)</math>は厳密解を用いて計算されているが、逐次計算時には中心座標を用いて<math>\vec a_n = A(\vec x_n)</math>のように計算されることに注意が必要である。大域誤差を計算する際には、この厳密解と近似解列との差は消えず、大域誤差に寄与する。

=== シンプルな例 ===
局所誤差と大域誤差との関係について洞察を得るため、厳密解と近似解が陽に書き下せるシンプルな例を考えてみる。このような例として標準的なものとして[[指数関数]]があげられる。

{{Mvar|w}}を定数として、線形微分方程式<math>\ddot x(t) = w^2 x(t)</math>を考える。この方程式の厳密解は<math>e^{wt}</math>および<math>e^{-wt}</math>である。

この微分方程式にStörmerの方法を適用すると、以下の線形[[漸化式]]を得る。

: <math>x_{n+1} - 2x_n + x_{n-1} = h^2 w^2 x_n,</math>

もしくは

: <math>x_{n+1} - 2\left(1 + \tfrac12(wh)^2\right) x_n + x_{n-1} = 0</math>

これは特性[[多項式の根]]を求めること、すなわち<math>q^2 - 2\left(1 + \tfrac12(wh)^2\right)q + 1 = 0</math>の解を求めることにより解くことができ、以下の根が得られる。

: <math>q_\pm = 1 + \tfrac12(wh)^2 \pm wh \sqrt{1 + \tfrac14(wh)^2}</math>

上の線形漸化式のbasis solutions{{訳語疑問点|date=2023年5月|基底解}}は<math>x_n = q_+^n</math>および<math>x_n = q_-^n</math>である。これらを厳密解と比較するため、テイラー展開すると以下を得る。

: <math>\begin{align}
q_+ &= 1 + \tfrac12(wh)^2 + wh\left(1 + \tfrac18(wh)^2 - \tfrac{3}{128}(wh)^4 + \mathcal O(h^6)\right)\\
    &= 1 + (wh) + \tfrac12(wh)^2 + \tfrac18(wh)^3 - \tfrac{3}{128}(wh)^5 + \mathcal O(h^7).
\end{align}</math>

この級数の指数関数<math>e^{wh}</math>との商は<math>1 - \tfrac1{24}(wh)^3 + \mathcal O(h^5)</math>となるため、以下のように書ける。

: <math>\begin{align}
q_+ &= \left(1 - \tfrac1{24}(wh)^3 + \mathcal O(h^5)\right)e^{wh}\\
    &= e^{-\frac{1}{24}(wh)^3 + \mathcal O(h^5)}\,e^{wh}
\end{align}</math>

ここから、最初のbasis solutionの誤差は以下のように算出される。

: <math>\begin{align}
x_n = q_+^{\;n}
 &= e^{-\frac{1}{24}(wh)^2\,wt_n + \mathcal O(h^4)}\,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(h^2 t_n e^{wt_n})
\end{align}</math>

したがって、局所離散化誤差は4次以上となるが、微分方程式が2階のため、大域誤差は2次で、時間的に指数的に増加する定数項をもつ。

=== 逐次計算の開始 ===
ベルレ法の最初のステップ、<math>n = 1</math>、<math>t = t_1 = \Delta t</math>で、<math>\vec x_2</math>を計算するためには<math>t = t_1</math>における位置ベクトル<math>\vec x_1</math>が必要となる。初期条件は<math>t_0 = 0</math>に対してのみ所与なので、一見これは問題をはらんでいるようにみえる。しかし、加速度<math>\vec a_0 = \vec A(\vec x_0)</math>は既知であるため、最初のステップは2次までのテイラー展開を用いて計算することができる。

: <math>
\vec x_1 = \vec{x}_0 + \vec{v}_0 \Delta t + \tfrac12 \vec a_0 \Delta t^2
\approx \vec{x}(\Delta t) + \mathcal{O}(\Delta t^3).
</math>

この最初のステップの誤差は<math>\mathcal O(\Delta t^3)</math>である。しかし、シミュレーションは長い時間、何ステップにもわたって行われ、時刻<math>t_n</math>におけるトータルの誤差は、<math>\vec x_n</math>と<math>\vec x(t_n)</math>との距離および<math>\tfrac{\vec x_{n+1} - \vec x_n}{\Delta t}</math>と<math>\tfrac{\vec x(t_{n+1}) - \vec x(t_n)}{\Delta t}</math>との比の両方でオーダー<math>\mathcal O(e^{Lt_n} \Delta t^2)</math>であり、最初のステップにおける誤差は無視できる。さらに言えば、この2次の大域誤差を得るためには初期誤差は最低でも3次である必要がある。

=== 非一様なタイムステップ ===
Störmer–Verlet法の弱点として、タイムステップ<math>\Delta t</math>が変化すると微分関数の近似解を与えなくなってしまう点である。この弱点は次の式を使うことにより修正できる<ref>{{cite web |last=Dummer |first=Jonathan |title=A Simple Time-Corrected Verlet Integration Method |url=http://lonesock.net/article/verlet.html}}</ref>。

: <math>
\vec x_{i+1}
= \vec x_i
+ (\vec x_i - \vec x_{i-1}) (\Delta t_i / \Delta t_{i-1})
+ \vec a_i \Delta t_i^2.
</math>

より厳密な導出は、<math>t_i</math>における2次までのテイラー展開に<math>t_{i+1} = t_i + \Delta t_i</math>と<math>t_{i-1} = t_i - \Delta t_{i-1}</math>を代入し、<math>\vec v_i</math>を消去して以下を得る。

: <math>
\frac{\vec x_{i+1} - \vec x_i}{\Delta t_i}
+ \frac{\vec x_{i-1} - \vec x_i}{\Delta t_{i-1}}
= \vec a_i\,\frac{\Delta t_{i} + \Delta t_{i-1}}2,
</math>

したがって、次の式を得る。

: <math>
\vec x_{i+1}
= \vec x_i
+ (\vec x_i - \vec x_{i-1}) \frac{\Delta t_i}{\Delta t_{i-1}}
+ \vec a_i\,\frac{\Delta t_{i} + \Delta t_{i-1}}2\,\Delta t_i
</math>

=== 速度の計算：Störmer–Verlet法 ===
基本的なStörmer方程式は速度を陽に与えないが、運動エネルギーなど特定の物理量を計算するためには速度が必要となる。このことは分子動力学法において、時刻<math>t</math>における運動エネルギーや瞬時温度を時刻<math>t + \Delta t</math>における位置を計算するまで計算できないという技術的な問題をひきおこす。この欠点は速度ベレ法を用いるか、[[平均値の定理]]を用いて以下のように位置から速度を推定することで対処することができる。

: <math>
\vec{v}(t)
=
\frac{\vec{x}(t + \Delta t) - \vec{x}(t - \Delta t)}{2\Delta t}
+ \mathcal{O}(\Delta t^2).
</math>

ここで、この速度項は時刻{{Math|''t'' + &Delta;''t''}}ではなく時刻<math>t</math>の速度であり、位置項の1ステップ後ろである。このことは<math>\vec v_n = \tfrac{\vec x_{n+1} - \vec x_{n-1}}{2\Delta t}</math>は<math>\vec{v}(t_n)</math>の2次近似であることを意味する。同じ議論として、タイムステップを半分<math>t_{n+1/2} = t_n + \tfrac12 \Delta t</math>にした<math>\vec v_{n+1/2} = \tfrac{\vec x_{n+1} - \vec x_n}{\Delta t}</math>は<math>\vec{v}(t_{n+1/2})</math>の2次近似である。

精度を犠牲にすれば、時刻{{Math|''t'' + &Delta;''t''}}における速度を次のように近似することもできる。

: <math>\vec{v}(t + \Delta t) = \frac{\vec{x}(t + \Delta t) - \vec{x}(t)}{\Delta t} + \mathcal{O}(\Delta t)</math>

== 速度ベルレ法 ==
関連する、より広く用いられているアルゴリズムとして'''速度ベルレ法'''があげられる<ref>{{Cite journal|last=Swope|first=William C.|last2=H. C. Andersen|last3=P. H. Berens|last4=K. R. Wilson|date=1 January 1982|title=A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters|journal=The Journal of Chemical Physics|volume=76|issue=1|pages=648 (Appendix)|bibcode=1982JChPh..76..637S|DOI=10.1063/1.442716}}</ref>。この手法は[[リープ・フロッグ法]]に似ているが、同時刻の位置と速度を計算する点（リープ・フロッグ法では名前の含意するとおり計算しない）で違う。速度ベルレ法とベルレ法は似ているが、陽に速度を扱い、初期ステップにおける速度を明示的に解く点で異なる。

: <math>\vec{x}(t + \Delta t) = \vec{x}(t) + \vec{v}(t)\, \Delta t + \frac{1}{2} \,\vec{a}(t) \Delta t^2,</math>
: <math>\vec{v}(t + \Delta t) = \vec{v}(t) + \frac{\vec{a}(t) + \vec{a}(t + \Delta t)}{2} \Delta t</math>

速度ベルレ法とベルレ法の誤差は同じオーダーであることは示すことができる。ベルレ法では2時刻における位置を記憶しておく必要があるため、1時刻における位置と速度を記憶する速度ベルレ法の方がメモリ消費量が多いとは限らない。

# <math>\vec{v}\left(t + \tfrac12\,\Delta t\right) = \vec{v}(t) + \tfrac12\,\vec{a}(t)\,\Delta t</math>を計算する。
# <math>\vec{x}(t + \Delta t) = \vec{x}(t) + \vec{v}\left(t + \tfrac12\,\Delta t\right)\, \Delta t</math>を計算する。
# <math>\vec{x}(t + \Delta t)</math>を用いて相互作用ポテンシャルを計算し、<math>\vec{a}(t + \Delta t)</math> を導出する。
# <math>\vec{v}(t + \Delta t) = \vec{v}\left(t + \tfrac12\,\Delta t\right) + \tfrac12\,\vec{a}(t + \Delta t)\Delta t</math>を計算する。

半ステップ速度計算を省く場合、以下のように簡略化される。

# <math>\vec{x}(t + \Delta t) = \vec{x}(t) + \vec{v}(t)\,\Delta t + \tfrac12 \,\vec{a}(t)\,\Delta t^2</math> を計算する。
# <math>\vec{x}(t + \Delta t)</math>を用いて相互作用ポテンシャルを計算し、<math>\vec{a}(t + \Delta t)</math>を導出する。
# <math>\vec{v}(t + \Delta t) = \vec{v}(t) + \tfrac12\,\big(\vec{a}(t) + \vec{a}(t + \Delta t)\big)\Delta t</math> を計算する。

ただし、加速度<math>\vec{a}(t + \Delta t)</math>が<math>\vec{x}(t + \Delta t)</math>のみに依存し、<math>\vec{v}(t + \Delta t)</math>に依存しないことを前提としている。

速度ベルレ法は、長期的にはリープ・フロッグ法と同様に{{仮リンク|半陰的オイラー法|en|Semi-implicit Euler method}}よりもオーダー1つ分良い近似である。速度の時刻が半ステップずれる以外は殆ど同じである。このことは上述のループのステップ3から初めて、ステップ2と4を組み合わせればステップ1における加速度項は取り除けることに注意すればすぐに証明することができる。唯一の違いは、速度ベルレ法のける中間速度が半陰的オイラー法における最終速度として扱われることである。

オイラー法の大域誤差はオーダーは1であるのに対して、速度ベルレ法の大域誤差のオーダーは[[中点法]]と同じく2である。加えて、もし加速度が[[ポテンシャル|保存力]]もしくは{{仮リンク|ハミルトニアン系|en|Hamiltonian system}}から定まる場合、エネルギーの近似値は厳密解における一定値のエネルギーのまわりを振動し、大域誤差は半陰的オイラー法ではオーダー1、速度ベルレ法およびリープ・フロッグ法ではオーダー2の範囲に収まる。系の運動量や角運動量などの[[シンプレクティック数値積分法]]で保存されるその他の量についても同じことが言える<ref name="Hairer2003">{{cite journal|last=Hairer|first=Ernst|last2=Lubich|first2=Christian|last3=Wanner|first3=Gerhard|year=2003|title=Geometric numerical integration illustrated by the Störmer/Verlet method|journal=Acta Numerica|volume=12|pages=399–450|bibcode=2003AcNum..12..399H|doi=10.1017/S0962492902000144|citeseerx=10.1.1.7.7106}}</ref>。

速度ベルレ法は、{{仮リンク|ニューマークのβ法|en|Newmark-beta method}}の{{Math|1=''&beta;''=0, ''&gamma;''=1/2}}とした特殊例である。

=== アルゴリズムの表現 ===
速度ベルレ法は3Dアプリケーション一般に有用なアルゴリズムであり、以下のようなC++実装により一般的に解ける。加速度の向きの変化をデモするため、簡略化された抗力を作用させているが、加速度が一定でない場合にのみ必要となる。<syntaxhighlight lang="c++" line="1">
struct Body
{
    Vec3d pos { 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 }; // 加速度の初期値は0
    double mass = 1.0; // 1kg
    double drag = 0.1; // rho*C*Area - 例示のための簡略化された抗力

    /**
     * Update pos and vel using "Velocity Verlet" integration
     * @param dt DeltaTime / time step [eg: 0.01]
     */
    void update(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);
        pos = new_pos;
        vel = new_vel;
        acc = new_acc;
    }

    Vec3d apply_forces() const
    {
        Vec3d grav_acc = Vec3d{0.0, 0.0, -9.81 }; // Z軸負方向に 9.81m/s^2
        Vec3d drag_force = 0.5 * drag * (vel * abs(vel)); // D = 0.5 * (rho * C * Area * vel^2)
        Vec3d drag_acc = drag_force / mass; // a = F/m
        return grav_acc - drag_acc;
    }
};
</syntaxhighlight>

== 誤差項 ==
ベルレ法の局所位置誤差は上述のとおり<math>O(\Delta t^4)</math>、局所速度誤差は<math>O(\Delta t^2)</math>である。

しかし、大域位置誤差は<math>O(\Delta t^2)</math>、大域速度誤差は<math>O(\Delta t^2)</math>である。これは次のように導出できる。

: <math>\mathrm{error}\bigl(x(t_0 + \Delta t)\bigr) = O(\Delta t^4)</math>

および

: <math>x(t_0 + 2\Delta t) = 2x(t_0 + \Delta t) - x(t_0) + \Delta t^2 x''(t_0 + \Delta t) + O(\Delta t^4)</math>

より以下が導かれる。

: <math>\mathrm{error}\bigl(x(t_0 + 2\Delta t)\bigr) = 2\mathrm{error}\bigl(x(t_0 + \Delta t)\bigr) + O(\Delta t^4) = 3\,O(\Delta t^4)</math>

同様に繰り返せば、以下を得る。

: <math>\mathrm{error}\bigl(x(t_0 + 3\Delta t)\bigl) = 6\,O(\Delta t^4)</math>
: <math>\mathrm{error}\bigl(x(t_0 + 4\Delta t)\bigl) = 10\,O(\Delta t^4)</math>
: <math>\mathrm{error}\bigl(x(t_0 + 5\Delta t)\bigl) = 15\,O(\Delta t^4)</math>

この数列の一般項は以下のように得られる。

: <math>\mathrm{error}\bigl(x(t_0 + n\Delta t)\bigr) = \frac{n(n+1)}{2}\,O(\Delta t^4)</math>

大域位置誤差は<math>T = n\Delta t</math>としたときの位置<math>x(t + T)</math>の誤差であるから、次を得る。

: <math>\mathrm{error}\bigl(x(t_0 + T)\bigr) = \left(\frac{T^2}{2\Delta t^2} + \frac{T}{2\Delta t}\right) O(\Delta t^4)</math>{{要出典|date=July 2018}}

したがって、時間間隔を一定とした場合の大域（累積）誤差は以下のオーダーとなる。

: <math>\mathrm{error}\bigr(x(t_0 + T)\bigl) = O(\Delta t^2)</math>

ベルレ法では速度は位置に依存して非累積的に算出されるため、大域速度誤差のオーダーも<math>O(\Delta t^2)</math>となる。

分子動力学シミュレーションでは通常、局所誤差よりも大域誤差のほうがはるかに重要であるため、ベルレ法はオーダー2の積分アルゴリズムとみなされる。

== 拘束 ==
拘束条件下にある多粒子系は、ベルレ法を用いるほうがオイラー法よりも単純に解くことができる。質点間の拘束条件はたとえば距離を特定値に拘束したり引力を印加したりする条件が考えられるが、粒子間を[[ばね]]で繋いだと考えることができる。無限に硬いばねを用いれば、このモデルはベルレ法で解くことができる。

一次元の場合、{{Mvar|i}}番目の質点の時刻{{Mvar|t}}における、拘束されない位置を<math>\tilde{x}_i^{(t)}</math>、実際の位置を<math>x_i^{(t)}</math>とすると、これらの間の関係は以下のようなアルゴリズムで得ることができる。

: <math>d_1 = x_2^{(t)} - x_1^{(t)},</math>

: <math>d_2 = \|d_1\|,</math>

: <math>d_3 = \frac{d_2 - r}{d_2},</math>

: <math>x_1^{(t + \Delta t)} = \tilde{x}_1^{(t + \Delta t)} + \frac{1}{2} d_1 d_3,</math>

: <math>x_2^{(t + \Delta t)} = \tilde{x}_2^{(t + \Delta t)} - \frac{1}{2} d_1 d_3.</math>

ベルレ法は位置を直接力に関係づけるため、速度を用いて問題を解くよりも使いやすい。

しかし、各粒子に複数の拘束条件が課される場合、問題が起きる。この問題を解決するには、シミュレーションのタイムステップを小さくしたり、タイムステップごとに決まったステップ数の拘束緩和を行ったり、残差が特定の値を下回るまで拘束を緩和する、などの方法でシミュレーションを実装する。

拘束条件を1次までで局所近似する場合、これは[[ガウス＝ザイデル法]]と同等となる。[[行列]]が小さい場合、[[LU分解]]の方が速いことが知られている。大きな系もクラスター（例:[[ラグドール物理|ラグドール]]=クラスター）に分解できることがある。クラスター内ではLU分解を用い、クラスター間にはガウス＝ザイデル法を用いればよい。行列のコードを再利用し、力の位置への依存性を1次局所近似することにより、ベルレ積分をより暗黙的にすることができる。

SuperLUを始めとした、疎行列を用いて複雑な問題を解くための洗練されたも存在する。また、例えば[[音波]]を形成することなく力を布を伝わらせるなどの特殊化された問題を取り扱うための特殊な技法もある<ref>{{Cite journal|last=Baraff|first=D.|last2=Witkin|first2=A.|date=1998|title=Large Steps in Cloth Simulation|url=https://www.cs.cmu.edu/~baraff/papers/sig98.pdf|journal=Computer Graphics Proceedings|volume=Annual Conference Series|pages=43–54}}</ref>。

{{仮リンク|ホロノミック拘束|en|Holonomic constraints}}を解くには、{{仮リンク|拘束条件|en|Constraint (computational chemistry)}}を扱えるアルゴリズムを用いる必要がある。

== 衝突応答 ==
[[衝突]]に応答する一つの方法として、ペナルティに基く方法が挙げられる。この方法では基本的に衝突した点に力を加える。問題は、加える力をどのように決めるかである。力が強すぎると物体は不安定になり、弱すぎると互いに貫通してしまう。衝突に応答する別の方法としては投影法があり、衝突した物体を最小の移動距離で相手の物体の内部から出すように動かすことで衝突に応答する。

ベルレ法では、後者の方法で衝突した場合の速度が自動的に決定する。しかし、このことは衝突の物理が自動的に満たされる(すなわち運動量の変化が現実的なものとなる)ことを意味しない。速度項を暗黙的に変化させる代わりに、衝突した物体の最終速度を陽に（直前のタイムステップにおいて記録された位置を変化させることにより）制御する必要がある。

新しい速度を決める最も単純な方法として、[[弾性衝突|完全弾性衝突]]もしくは[[非弾性衝突|完全非弾性衝突]]を仮定することが挙げられる。より複雑な制御方法では、[[反発係数]]を用いることもある。

== 関連項目 ==

* [[CFL条件]]
* {{仮リンク|エネルギードリフト|en|Energy drift}}
* [[シンプレクティック数値積分法]]
* [[リープ・フロッグ法]]
* {{仮リンク|ビーマンのアルゴリズム|en|Beeman's algorithm}}

== 出典 ==
<references>
<ref name="Hairer2003">{{cite journal|last=Hairer|first=Ernst|last2=Lubich|first2=Christian|last3=Wanner|first3=Gerhard|year=2003|title=Geometric numerical integration illustrated by the Störmer/Verlet method|journal=Acta Numerica|volume=12|pages=399–450|bibcode=2003AcNum..12..399H|doi=10.1017/S0962492902000144|citeseerx=10.1.1.7.7106}}</ref>
<ref name="Verlet1967">{{cite journal|last=Verlet|first=Loup|year=1967|title=Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard−Jones Molecules|journal=Physical Review|volume=159|issue=1|pages=98–103|bibcode=1967PhRv..159...98V|doi=10.1103/PhysRev.159.98|doi-access=free}}</ref>
</references>

== 外部リンク ==

* [https://bitbucket.org/craigmit/verlet Verlet Integration Demo and Code as a Java Applet]
* [https://www.cs.cmu.edu/afs/cs/academic/class/15462-s13/www/lec_slides/Jakobsen.pdf Advanced Character Physics by Thomas Jakobsen]
* [http://www.ch.embnet.org/MD_tutorial/pages/MD.Part1.html Theory of Molecular Dynamics Simulations] – bottom of page


{{DEFAULTSORT:へれのほうほう}}
[[Category:計算物理学]]
[[Category:分子動力学]]
[[Category:数値微分方程式]]
[[Category:エポニム]]