基本
の形の2階微分方程式を初期条件 および のもとで、 をステップサイズとして、 における数値的近似解 は、下のようなアルゴリズムで求めることができる。
- と置く。
- n = 1、2、...について次を繰り返す。
-
運動方程式
保存系におけるニュートンの運動方程式は
-
もしくは粒子毎に
-
と表わされる。ここで、
- tは時刻、
- はN個の物体の位置ベクトル、
- Vはスカラーポテンシャル関数、
- Fは負のポテンシャル勾配、すなわち粒子にはたらく力の集合、
- M質量行列で、各粒子の質量 を対角要素にもつ行列である。
この等式は、相互作用する分子群や惑星の軌道などさまざまな物理系が、さまざまなポテンシャル関数Vを決めたときどのように時間発展するかを記述するために用いることができる。
右辺に質量を移し、粒子系の構造を忘れることにすると、上の式は以下のように単純化することができる。
-
位置に依存する加速度を表す適切なベクトル値関数Aを使用します。通常、も与えられます。ここで、ベクトル値関数Aは位置依存加速度をあらわす。通常は、初期位置 および初速度 も所与である。
ベレ積分(速度なし)
この初期値問題を離散化し、数値的に解くため、タイムステップ を選び、サンプル点列 を考える。このとき、問題は厳密解 をよく近似する点群列 をもとめることに帰着する。
オイラー法が1階微分方程式中の1次微分を前進差分近似するのに対し、ベレ積分は2次微分を中心差分近似すると見ることができる。
-
「Störmerの方法」[3]で用いられる形式の「ベレ積分」は、 この式を用い、速度を使わずに以前の2つの位置から次の座標を与える。
-
離散化誤差
このアルゴリズムは本質的に時間反転対称性を持ち、離散化の際に奇数次の項が消え、 について3次の項がなくなるため、局所誤差の大きさを低減することができる。局所誤差は厳密解 を代入し、 におけるテイラー展開から をそれぞれ計算することにより定量できる。
-
ここで、 は位置、 は速度、 は加速度、 は躍度である。
これら2つの級数を足し合わせると、以下の式を得る。
-
この式から、テイラー展開の1次および3次の項がうち消しあい、ベレ積分の精度が単なるテイラー展開よりも1次高くなっていることがわかる。
ここで、加速度 は厳密解を用いて計算されているが、逐次計算時には中心座標を用いて のように計算されることに注意が必要である。大域誤差を計算する際には、この厳密解と近似解列との差は消えず、大域誤差に寄与する。
シンプルな例
局所誤差と大域誤差との関係について洞察を得るため、厳密解と近似解が陽に書き下せるシンプルな例を考えてみる。このような例として標準的なものとして指数関数があげられる。
wを定数として、線形微分方程式 を考える。この方程式の厳密解は および である。
この微分方程式にStörmerの方法を適用すると、以下の線形漸化式を得る。
-
もしくは
-
これは特性多項式の根を求めること、すなわち の解を求めることにより解くことができ、以下の根が得られる。
-
上の線形漸化式のbasis solutions[訳語疑問点]は および である。これらを厳密解と比較するため、テイラー展開すると以下を得る。
-
この級数の指数関数 との商は となるため、以下のように書ける。
-
ここから、最初のbasis solutionの誤差は以下のように算出される。
-
したがって、局所離散化誤差は4次以上となるが、微分方程式が2階のため、大域誤差は2次で、時間的に指数的に増加する定数項をもつ。
逐次計算の開始
ベルレ法の最初のステップ、 、 で、 を計算するためには における位置ベクトル が必要となる。初期条件は に対してのみ所与なので、一見これは問題をはらんでいるようにみえる。しかし、加速度 は既知であるため、最初のステップは2次までのテイラー展開を用いて計算することができる。
-
この最初のステップの誤差は である。しかし、シミュレーションは長い時間、何ステップにもわたって行われ、時刻 におけるトータルの誤差は、 と との距離および と との比の両方でオーダー であり、最初のステップにおける誤差は無視できる。さらに言えば、この2次の大域誤差を得るためには初期誤差は最低でも3次である必要がある。
非一様なタイムステップ
Störmer–Verlet法の弱点として、タイムステップ が変化すると微分関数の近似解を与えなくなってしまう点である。この弱点は次の式を使うことにより修正できる[4]。
-
より厳密な導出は、 における2次までのテイラー展開に と を代入し、 を消去して以下を得る。
-
したがって、次の式を得る。
-
速度の計算:Störmer–Verlet法
基本的なStörmer方程式は速度を陽に与えないが、運動エネルギーなど特定の物理量を計算するためには速度が必要となる。このことは分子動力学法において、時刻 における運動エネルギーや瞬時温度を時刻 における位置を計算するまで計算できないという技術的な問題をひきおこす。この欠点は速度ベレ法を用いるか、平均値の定理を用いて以下のように位置から速度を推定することで対処することができる。
-
ここで、この速度項は時刻t + Δtではなく時刻 の速度であり、位置項の1ステップ後ろである。このことは は の2次近似であることを意味する。同じ議論として、タイムステップを半分 にした は の2次近似である。
精度を犠牲にすれば、時刻t + Δtにおける速度を次のように近似することもできる。
-
速度ベルレ法
関連する、より広く用いられているアルゴリズムとして速度ベルレ法があげられる[5]。この手法はリープ・フロッグ法に似ているが、同時刻の位置と速度を計算する点(リープ・フロッグ法では名前の含意するとおり計算しない)で違う。速度ベルレ法とベルレ法は似ているが、陽に速度を扱い、初期ステップにおける速度を明示的に解く点で異なる。
-
-
速度ベルレ法とベルレ法の誤差は同じオーダーであることは示すことができる。ベルレ法では2時刻における位置を記憶しておく必要があるため、1時刻における位置と速度を記憶する速度ベルレ法の方がメモリ消費量が多いとは限らない。
- を計算する。
- を計算する。
- を用いて相互作用ポテンシャルを計算し、 を導出する。
- を計算する。
半ステップ速度計算を省く場合、以下のように簡略化される。
- を計算する。
- を用いて相互作用ポテンシャルを計算し、 を導出する。
- を計算する。
ただし、加速度 が のみに依存し、 に依存しないことを前提としている。
速度ベルレ法は、長期的にはリープ・フロッグ法と同様に半陰的オイラー法(英語版)よりもオーダー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; }};
誤差項
拘束
衝突応答
衝突に応答する一つの方法として、ペナルティに基く方法が挙げられる。この方法では基本的に衝突した点に力を加える。問題は、加える力をどのように決めるかである。力が強すぎると物体は不安定になり、弱すぎると互いに貫通してしまう。衝突に応答する別の方法としては投影法があり、衝突した物体を最小の移動距離で相手の物体の内部から出すように動かすことで衝突に応答する。
ベルレ法では、後者の方法で衝突した場合の速度が自動的に決定する。しかし、このことは衝突の物理が自動的に満たされる(すなわち運動量の変化が現実的なものとなる)ことを意味しない。速度項を暗黙的に変化させる代わりに、衝突した物体の最終速度を陽に(直前のタイムステップにおいて記録された位置を変化させることにより)制御する必要がある。
新しい速度を決める最も単純な方法として、完全弾性衝突もしくは完全非弾性衝突を仮定することが挙げられる。より複雑な制御方法では、反発係数を用いることもある。
関連項目
出典
外部リンク