ベレの方法

**ベレの方法**（ベレのほうほう、）は、ニュートンの運動方程式を数値積分する手法の一つ。**ベレのアルゴリズム**、**ベレ法**、**ベルレ積分法**()、**ベルレの方法**などの呼び方もある。分子動力学法における粒子の軌跡（トラジェクトリ）のシミュレーションやコンピュータグラフィックスに頻繁に用いられる。1791年にジャン＝バティスト・ジョゼフ・ドランブルが用いたのが最初で、その後も何回も再発見されているが、1960年代にLoup Verletが分子動力学法に用いてから広く使われるようになった。1909年にハレー彗星の軌道を計算するためにフィリップ・コーウェルとアンドリュー・クロンメリンにより用いられたり、1907年に磁場中の荷電粒子のトラジェクトリを研究するためにCarl Størmerにより用いられたりしている（そのため、**Störmerの方法**という異称もある）。この手法の特徴として、数値的安定性が高く、時間反転対称性を持つことや位相空間上のシンプレクティック形式を保存するなど物理系において重要な性質を持つうえ、オイラー法と比べてもほとんど計算コストが増えないことが挙げられる。

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

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


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

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

もしくは粒子毎に

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

と表わされる。ここで、

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

この等式は、相互作用する分子群や惑星の軌道などさまざまな物理系が、さまざまなポテンシャル関数を決めたときどのように時間発展するかを記述するために用いることができる。

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

: \ddot{\vec x}(t) = A\big(\vec x(t)\big)

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

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

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

: 
\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)


「Störmerの方法」で用いられる形式の「ベレ積分」は、 この式を用い、速度を使わずに以前の2つの位置から次の座標を与える。

: 
\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)


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

: \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}

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

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

: \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).

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

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

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

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

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

: x_{n+1} - 2x_n + x_{n-1} = h^2 w^2 x_n,

もしくは

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

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

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

上の線形漸化式のbasis solutionsはx_n = q_+^nおよびx_n = q_-^nである。これらを厳密解と比較するため、テイラー展開すると以下を得る。

: \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}

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

: \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}

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

: \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}

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

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

: 
\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).


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

### 非一様なタイムステップ
Störmer–Verlet法の弱点として、タイムステップ\Delta tが変化すると微分関数の近似解を与えなくなってしまう点である。この弱点は次の式を使うことにより修正できる。

: 
\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.


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

: 
\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,


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

: 
\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


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

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


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

精度を犠牲にすれば、時刻における速度を次のように近似することもできる。

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

## 速度ベルレ法
関連する、より広く用いられているアルゴリズムとして**速度ベルレ法**があげられる。この手法はリープ・フロッグ法に似ているが、同時刻の位置と速度を計算する（リープ・フロッグ法では名前の含意するとおり計算しない）点で違う。速度ベルレ法とベルレ法は似ているが、陽に速度を扱い、初期ステップにおける速度を明示的に解く点で異なる。

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

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

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

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

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

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

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

オイラー法の大域誤差はオーダーは1であるのに対して、速度ベルレ法の大域誤差のオーダーは中点法と同じく2である。加えて、もし加速度が保存力もしくはから定まる場合、エネルギーの近似値は厳密解における一定値のエネルギーのまわりを振動し、大域誤差は半陰的オイラー法ではオーダー1、速度ベルレ法およびリープ・フロッグ法ではオーダー2の範囲に収まる。系の運動量や角運動量などのシンプレクティック数値積分法で保存されるその他の量についても同じことが言える。

速度ベルレ法は、のとした特殊例である。

### アルゴリズムの表現
速度ベルレ法は3Dアプリケーション一般に有用なアルゴリズムであり、以下のようなC++実装により一般的に解ける。加速度の向きの変化をデモするため、簡略化された抗力を作用させているが、加速度が一定でない場合にのみ必要となる。
```c++
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;
    }
};
```


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

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

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

および

: 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)

より以下が導かれる。

: \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)

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

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

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

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

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

: \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)

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

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

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

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

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

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

: d_1 = x_2^{(t)} - x_1^{(t)},

: d_2 = \|d_1\|,

: d_3 = \frac{d_2 - r}{d_2},

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

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

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

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

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

SuperLUを始めとした、疎行列を用いて複雑な問題を解くための洗練されたも存在する。また、例えば音波を形成することなく力を布を伝わらせるなどの特殊化された問題を取り扱うための特殊な技法もある。

を解くには、を扱えるアルゴリズムを用いる必要がある。

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

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

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

## 関連項目

* CFL条件
* 
* シンプレクティック数値積分法
* リープ・フロッグ法
* 

## 出典




## 外部リンク

* Verlet Integration Demo and Code as a Java Applet
* Advanced Character Physics by Thomas Jakobsen
* Theory of Molecular Dynamics Simulations – bottom of page