ベレの方法

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

基本

${\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+{\frac {1}{2}}{\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 )},

と表わされる。ここで、

tは時刻、

{\vec {x}}(t)={\big (}{\vec {x}}_{1}(t),\ldots ,{\vec {x}}_{N}(t){\big )}

はN個の物体の位置ベクトル、

Vはスカラーポテンシャル関数、

Fは負のポテンシャル勾配、すなわち粒子にはたらく力の集合、

M質量行列で、各粒子の質量

m_{k}

を対角要素にもつ行列である。

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

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

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

位置に依存する加速度を表す適切なベクトル値関数Aを使用します。通常、も与えられます。ここで、ベクトル値関数Aは位置依存加速度をあらわす。通常は、初期位置 ${\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の方法」^[3]で用いられる形式の「ベレ積分」は、この式を用い、速度を使わずに以前の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{aligned}{\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{aligned}}

ここで、 ${\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})$ のように計算されることに注意が必要である。大域誤差を計算する際には、この厳密解と近似解列との差は消えず、大域誤差に寄与する。

シンプルな例

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

wを定数として、線形微分方程式 ${\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+{\tfrac {1}{2}}(wh)^{2}\right)x_{n}+x_{n-1}=0

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

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

上の線形漸化式のbasis solutions^{[訳語疑問点]}は $x_{n}=q_{+}^{n}$ および $x_{n}=q_{-}^{n}$ である。これらを厳密解と比較するため、テイラー展開すると以下を得る。

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

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

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

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

{\begin{aligned}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{aligned}}

したがって、局所離散化誤差は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+{\tfrac {1}{2}}{\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$ が変化すると微分関数の近似解を与えなくなってしまう点である。この弱点は次の式を使うことにより修正できる^[4]。

{\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 + Δtではなく時刻 $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}+{\tfrac {1}{2}}\Delta t$ にした ${\vec {v}}_{n+1/2}={\tfrac {{\vec {x}}_{n+1}-{\vec {x}}_{n}}{\Delta t}}$ は ${\vec {v}}(t_{n+1/2})$ の2次近似である。

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

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

速度ベルレ法

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

{\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+{\tfrac {1}{2}}\,\Delta t\right)={\vec {v}}(t)+{\tfrac {1}{2}}\,{\vec {a}}(t)\,\Delta t$ を計算する。
${\vec {x}}(t+\Delta t)={\vec {x}}(t)+{\vec {v}}\left(t+{\tfrac {1}{2}}\,\Delta t\right)\,\Delta t$ を計算する。
${\vec {x}}(t+\Delta t)$ を用いて相互作用ポテンシャルを計算し、 ${\vec {a}}(t+\Delta t)$ を導出する。
${\vec {v}}(t+\Delta t)={\vec {v}}\left(t+{\tfrac {1}{2}}\,\Delta t\right)+{\tfrac {1}{2}}\,{\vec {a}}(t+\Delta t)\Delta t$ を計算する。

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

${\vec {x}}(t+\Delta t)={\vec {x}}(t)+{\vec {v}}(t)\,\Delta t+{\tfrac {1}{2}}\,{\vec {a}}(t)\,\Delta t^{2}$ を計算する。
${\vec {x}}(t+\Delta t)$ を用いて相互作用ポテンシャルを計算し、 ${\vec {a}}(t+\Delta t)$ を導出する。
${\vec {v}}(t+\Delta t)={\vec {v}}(t)+{\tfrac {1}{2}}\,{\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の範囲に収まる。系の運動量や角運動量などのシンプレクティック数値積分法で保存されるその他の量についても同じことが言える^[6]。

速度ベルレ法は、ニューマークのβ法（英語版）のβ=0, γ=1/2とした特殊例である。

アルゴリズムの表現

速度ベルレ法は3Dアプリケーション一般に有用なアルゴリズムであり、以下のような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の積分アルゴリズムとみなされる。

拘束

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

一次元の場合、i番目の質点の時刻tにおける、拘束されない位置を ${\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を始めとした、疎行列を用いて複雑な問題を解くための洗練されたも存在する。また、例えば音波を形成することなく力を布を伝わらせるなどの特殊化された問題を取り扱うための特殊な技法もある^[7]。

ホロノミック拘束（英語版）を解くには、拘束条件（英語版）を扱えるアルゴリズムを用いる必要がある。

衝突応答

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

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

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

出典

^ Verlet, Loup (1967). “Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard−Jones Molecules”. Physical Review 159 (1): 98–103. Bibcode: 1967PhRv..159...98V. doi:10.1103/PhysRev.159.98.
^ Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. (2007). “Section 17.4. Second-Order Conservative Equations”. Numerical Recipes: The Art of Scientific Computing (3rd ed.). New York: Cambridge University Press. ISBN 978-0-521-88068-8. http://apps.nrbook.com/empanel/index.html#pg=928
^ webpage Archived 2004-08-03 at the Wayback Machine. with a description of the Störmer method.
^ Dummer, Jonathan. “A Simple Time-Corrected Verlet Integration Method”. Template:Cite webの呼び出しエラー：引数 accessdate は必須です。
^ Swope, William C.; H. C. Andersen; P. H. Berens; K. R. Wilson (1 January 1982). “A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters”. The Journal of Chemical Physics 76 (1): 648 (Appendix). Bibcode: 1982JChPh..76..637S. doi:10.1063/1.442716.
^ Hairer, Ernst; Lubich, Christian; Wanner, Gerhard (2003). “Geometric numerical integration illustrated by the Störmer/Verlet method”. Acta Numerica 12: 399–450. Bibcode: 2003AcNum..12..399H. doi:10.1017/S0962492902000144.
^ Baraff, D.; Witkin, A. (1998). “Large Steps in Cloth Simulation”. Computer Graphics Proceedings Annual Conference Series: 43–54. https://www.cs.cmu.edu/~baraff/papers/sig98.pdf.