【Unity】Velocity Verlet（速度ベルレ法）を使って物体の位置を計算する

Velocity Verlet（速度ベルレ法）を使って物体の位置を計算する方法です。

Velocity Verlet?
式(1)を導出してみる
実装
結果
参考

Velocity Verlet?

3DCGで物理シミュレーションを行う際には、漸化式を使った逐次計算が用いられます。
代表的な計算手法にオイラー法やルンゲクッタ法といったものがあります。

light11.hatenadiary.com

これらの計算方法は常微分方程式をプログラムで解くときの手法であるため位置の計算に限らず適用できます。
しかしこれらは時間が経つにつれて誤差がどんどん大きくなる（エネルギーが保存されない）という特徴があるので、
シミュレーションの用途によっては別の手法が必要になります。

Velocity Verletは位置を計算するための手法であり、エネルギーを保存しつつ比較的低負荷な計算できる特徴を持ちます。
ある時点 $t$ における座標 $r(t)$ と速度 $v(t)$ を使って単位時間 $\Delta t$ 後の座標や速度が求めるとき、式は以下のように表されます。

${ \displaystyle r(t + \Delta t) =r(t)+\Delta tv(t)+\frac{\Delta t^{2}}{2}\frac{F(t)}{m} \tag{1} }$

${ \displaystyle v(t + \Delta t) =v(t)+\frac{\Delta t}{2m}\left ( F(t) + F(t + \Delta t) \right ) \tag{2} }$

式(1)を導出してみる

さて次にもう少し理解を深めるため、前節の式(1)を導出してみます。

まず、 $r(t + \Delta t)$ をテイラー展開します。
テイラー展開は以下の式で表されます。

${ \displaystyle f(x)=f(a)+\frac{f(a)'}{1!}(x-a)+\frac{f(a)''}{2!}(x-a)^{2}+\cdot \cdot \cdot \tag{3} }$

ので、 $r(t + \Delta t)$ をテイラー展開すると以下の式が得られます。

${ \displaystyle r(t+ \Delta t)=r(a)+\frac{dr(a)}{dt}(t+ \Delta t -a)+\frac{1}{2}\frac{d^{2}r(a)}{dr^{2}}(t+ \Delta t -a)^{2}+\cdot \cdot \cdot \tag{4} }$

ここで、位置を時間で微分したものは速度 $v$ に、速度を時間で微分したものは加速度 $\alpha=\frac{F}{m}$ になります。
またいま $a=t$ を考えるとすると、式(3)は以下のようになります。

${ \displaystyle r(t+ \Delta t)=r(t)+v(t)\Delta t+\frac{1}{2}\frac{F(t)}{m}\Delta t ^{2}+O(\Delta t^{3}) \tag{5} }$

また同様に $r(t- \Delta t)$ は以下の式で表せます。

${ \displaystyle r(t- \Delta t)=r(t)-v(t)\Delta t+\frac{1}{2}\frac{F(t)}{m}\Delta t ^{2}-O(\Delta t^{3}) \tag{6} }$

$O$ はランダウの記号です。
$\Delta t$ は十分に小さい値であることを想定しているので以降ではこの項は無視します。

次に $r(t+ \Delta t)+r(t- \Delta t)$ と $r(t+ \Delta t)-r(t- \Delta t)$ を求めます。

${ \displaystyle r(t+ \Delta t)+r(t- \Delta t)=2r(t)+\frac{F(t)}{m}\Delta t ^{2} \tag{7} }$

${ \displaystyle r(t+ \Delta t)-r(t- \Delta t)=2v(t)\Delta t \tag{8} }$

式(7)を変形すると以下の式が得られます。

${ \displaystyle r(t+ \Delta t)=2r(t)-r(t-\Delta t)+\frac{F(t)}{m}\Delta t^{2} \tag{9} }$

これをベルレの式と呼びます。
また式(8)を変形すると以下の式が得られます。

${ \displaystyle v(t)=\frac{1}{2 \Delta t}(r(t + \Delta t) - r(t - \Delta t)) \tag{10} }$

ここでさらに式(10)を以下の形に変形します。

${ \displaystyle r(t- \Delta t)=t(t+ \Delta t)-2v(t)\Delta t \tag{11} }$

これを式(9)に代入すると以下の式が得られます。

${ \displaystyle r(t + \Delta t) =r(t)+\Delta tv(t)+\frac{\Delta t^{2}}{2}\frac{F(t)}{m} \tag{12} }$

これは式(1)と一致します。
このようにVelocity Verletは二次精度の解法であることが確認できました。

実装

上述の内容をそのまま実装すると以下のようになります。
単位時間（deltaTime）は固定します。

using UnityEngine;

public class VelocityVerlet : MonoBehaviour
{
    [SerializeField] private float _mass = 1.0f;

    [SerializeField] private float _forceFactor = 500f;

    private Vector3 Gravity { get; } = new Vector3(0f, -9.8f, 0f);
    
    private Vector3 _velocity;
    private Vector3 _force;
    private const int TARGET_FRAME_RATE = 60;
    private const float SIMULATE_DELTA_TIME = 1.0f / TARGET_FRAME_RATE;

    private void Start()
    {
        Application.targetFrameRate = TARGET_FRAME_RATE;
        _force = Gravity * _mass;
    }

    private void Update()
    {
        // このフレームにかかる力を計算
        var force = Gravity * _mass;
        if (Input.GetKeyDown(KeyCode.Space))
        {
            // スペースキーが押されたら上向きの力を加える
            force += Vector3.up * _forceFactor;
        }
        
        transform.position = transform.position + SIMULATE_DELTA_TIME * _velocity + SIMULATE_DELTA_TIME * SIMULATE_DELTA_TIME * _force / _mass / 2.0f;
        
        _velocity = _velocity + SIMULATE_DELTA_TIME * (_force + force) / _mass / 2.0f;

        _force = force;
    }
}

結果

前節のスクリプトを適当なオブジェクトにアタッチして再生します。
スペースキーを押下すると上向きの力が加わります。

f:id:halya_11:20200406214734g:plain — 結果

参考

en.wikipedia.org