[Unity] シンプレクティック法でオイラー陽解法の精度改善

はじめに

[Unity] オイラー陽解法から始める物理シミュ入門 の続きの記事です。

1. 概要

前回紹介したオイラー陽解法では徐々に誤差が大きくなってきてしまう問題がありました。これを簡単な変更だけでかなり改善することができるシンプレクティック法について書きたいと思います。

先に比較を載せます。グレー色の球が初期位置です。こうして見てみるとただのオイラー陽解法(青球)がすぐずれていくのがよくわかりますね...。一方で今回紹介する方法(赤球)ではほぼほぼ変わっていません。

2. おさらい

2.1. オイラー陽解法

詳細については[Unity] オイラー陽解法から始める物理シミュ入門に任せて割愛しますが、必要な情報だけ簡単におさらいします。

オイラー陽解法では、少し未来の $t+Δtt+\Delta t$ の値を近似によって求めるものでした。
その近似式は以下になります。

$$x(t+Δt)=x(t)+x′(t)Δt x(t+\Delta t) = x(t) + x'(t)\Delta t$$

2.2. コードで書ける形にする

微分というのは極小時間で見た時の変化量であり、 距離xの時間微分は速度vとなり、速度vの時間微分は加速度aとなります。

距離の微分 $x(t)′x(t)'$ を速度 $v(t)v(t)$ で書き直してみると以下になります。

$$x(t+Δt)=x(t)+v(t)Δt x(t+\Delta t) = x(t) + v(t)\Delta t$$

速度 $v(t)v(t)$ も必要なので、同様にオイラー陽解法で求めます。また速度の微分 $v(t)′v(t)'$ も加速度 $aa$ に置き換えておきます。

$$v(t+Δt)=v(t)+aΔtx(t+Δt)=x(t)+v(t)Δt \begin{aligned} v(t+\Delta t) &= v(t) + a\Delta t \\ x(t+\Delta t) &= x(t) + v(t)\Delta t \end{aligned}$$

さらに 加速度 $aa$ は運動方程式 $ma=Fma = F$ から求めることができました。今回も振り子の運動方程式を例に見ていきます。振り子の運動方程式は以下です。

https://w3e.kanazawa-it.ac.jp/math/physics/simulation/henkan-tex.cgi?target=/math/physics/simulation/simple_pendulum.html

$$F=ma=−mgsin⁡θ≒−mglxa=−glx F = ma = -mg \sin \theta \fallingdotseq - \frac{mg}{l}x \\ a = - \frac{g}{l}x$$

なので、もう一度近似式を書き直してみると以下になります。

$$v(t+Δt)=v(t)+(−glx(t))Δtx(t+Δt)=x(t)+v(t)Δt \begin{aligned} v(t+\Delta t) &= v(t) + (- \frac{g}{l}x(t))\Delta t \\ x(t+\Delta t) &= x(t) + v(t)\Delta t \end{aligned}$$

これで定数と時刻 $tt$ の値から時刻 $t+Δtt+\Delta t$ の値を計算できるようになりました。

2.3. コードに落とし込む

数式が苦手な場合はとっつきにくい気もしますが、コードを見てみると実際の計算は各々1行で数式のままですね。

以上がオイラー陽解法でシミュレーションのおさらいです。

void FixedUpdate()
{
	// [説明変数] 定数
	float dt = Time.deltaTime;
	float g = Physics.gravity.magnitude;
	float l = 5;
		
	// [説明変数] 現在の値
	float x_t = _x;
	float v_t = _v;
		
	// 加速度
	float a = -g / l * x_t;
		
	// dt秒先の速度を計算
	_v = v_t + a * dt;
		
	// dt秒先の座標を計算
	_x = x_t + v_t * dt;
}

3. シンプレクティック法

とても簡単です！！！！！

$$v(t+Δt)=v(t)+aΔtx(t+Δt)=x(t)+v(t)Δt \begin{aligned} v(t+\Delta t) &= v(t) + a\Delta t \\ x(t+\Delta t) &= x(t) + v(t)\Delta t \end{aligned}$$

を

$$v(t+Δt)=v(t)+aΔtx(t+Δt)=x(t)+v(t+Δt)Δt \begin{aligned} v(t+\Delta t) &= v(t) + a\Delta t \\ x(t+\Delta t) &= x(t) + v(t+\Delta t)\Delta t \end{aligned}$$

にするだけです。以上！

オイラー陽解法では、 $x(t+Δt)x(t+\Delta t)$ について計算するときに現在の速度 $v(t)v(t)$ を使っていましたが、計算したばかりの速度 $v(t+Δt)v(t+\Delta t)$ を代わりに使うだけです。

載せるまでもないかもしれませんが、コードも載せておきます。

// dt秒先の速度を計算
_v = v_t + a * dt;
		
// dt秒先の座標を計算
_x = x_t + _v * dt;

コードだけ分かればOK！という方は実装についてのお話はここまでです。
以降は「なんで？？本当に？？」という説明をしていくので、実装の話はでてきません。

4. 本当に？？

理論に関してはかなり数学の知識が必要そうでかなり敷居が高そうでした...。ただ、なぜ誤差が広がっていかないかについてはこちらのサイトに分かりやすい説明がありました。

だいたいこんな感じかと思います。

1. オイラー陽解法/シンプレクティック法の式を行列の演算として考える
2. その行列の行列式を求めると、その行列変換による拡大率が分かる
3. オイラー陽解法の行列式では、エネルギーが拡大され発散することが分かる
4. シンプレクティック法の行列式では、エネルギーが保存されることが分かる

ただ、ちょっと理解度が怪しいので話半分で見ていただいて、間違いなどあれば教えてもらえると非常に助かります。

4.1. 行列で表す

オイラー陽解法の数式を行列として表してみます。

$$v(t+Δt)=v(t)+aΔtx(t+Δt)=x(t)+v(t)Δt \begin{aligned} v(t+\Delta t) &= v(t) + a\Delta t \\ x(t+\Delta t) &= x(t) + v(t)\Delta t \end{aligned}$$

ここで、 $aa$ は振り子の運動方程式から求めることができました。

$$a=−glx(t) a = - \frac{g}{l}x(t)$$

$$v(t+Δt)=v(t)+(−glx(t))Δtx(t+Δt)=x(t)+v(t)Δt \begin{aligned} v(t+\Delta t) &= v(t) + (- \frac{g}{l}x(t))\Delta t \\ x(t+\Delta t) &= x(t) + v(t)\Delta t \end{aligned}$$

これは以下のように行列式の形で書くことができます。

$$(v(t+Δt)x(t+Δt))=(abcd)(v(t)x(t))=(av(t)+bx(t)cv(t)+dx(t)) \begin{aligned} \begin{pmatrix} v(t+\Delta t) \\ x(t+\Delta t) \\ \end{pmatrix} &= \begin{pmatrix} a & b \\ c & d \\ \end{pmatrix} \begin{pmatrix} v(t) \\ x(t) \\ \end{pmatrix} \\ &= \begin{pmatrix} av(t) + bx(t) \\ cv(t) + dx(t) \\ \end{pmatrix} \end{aligned}$$

$v(t+Δt)v(t+\Delta t)$ と $x(t+Δt)x(t+\Delta t)$ の式をベクトル表記すると、上記の行列の形そのままなので分かりやすいですね。

$$(v(t+Δt)x(t+Δt))=((1)⋅v(t)+(−glΔt)⋅x(t)(1)⋅v(t)+(1)⋅x(t)) \begin{pmatrix} v(t+\Delta t) \\ x(t+\Delta t) \\ \end{pmatrix} = \begin{pmatrix} (1)\cdot v(t) + (-\frac{g}{l}\Delta t)\cdot x(t) \\ (1)\cdot v(t) + (1)\cdot x(t) \\ \end{pmatrix}$$

見比べれば $a,b,c,da, b, c, d$ の対応が簡単に見つけられると思います。

$$A=(abcd)=(1−glΔt11) A = \begin{pmatrix} a & b \\ c & d \\ \end{pmatrix} = \begin{pmatrix} 1 & -\frac{g}{l}\Delta t \\ 1 & 1 \\ \end{pmatrix} \\$$

よって $u=Aqu = Aq$ のように行列 $AA$ の積という形で表記できました。

4.2. 行列積による拡大率を見る

この行列変換によって速度と座標が変換されます。速度が変わるということは運動量が変化します。運動量は $p=mvp=mv$ であり、式からも速度 $vv$ に比例して運動量が変化することが分かります。

その行列の性質を見ることで運動量がどう増えるのか減るのか、保存されるのか、というのを見ます。また、座標に関しても時間変化を見た時の変化量(図示したときの面積)が保たれる、ということだと思います。

その行列の性質というのが行列式です。

https://ja.wikipedia.org/wiki/行列式

線型変換に対して線形空間の拡大率ということができる。

A(Xu, Xv) = (ad − bc)A(u, v) が成り立っているが、これは X の定める線型変換によって平面内の図形の面積が (ad − bc)-倍される、と解釈できる。

という感じに、行列の変換によって対象が何倍されるか拡大率を表す値です。

行列式 $∣A∣|A|$ は式といっても以下のようなスカラー値です。

$$∣A∣=∣abcd∣=ad−bc |A| = \begin{vmatrix} a & b \\ c & d \\ \end{vmatrix} = ad - bc$$

オイラー陽解法とシンプレクティック法の変換行列の行列式を調べて拡大率を見ることで発散するのか保持されるのかを見ていきます。図形的に見るならば、変換行列によって四角形が回転したり形が変わっても、その四角形の面積が保持されているということになります。

4.3. オイラー陽解法の拡大率を見てみよう

4.1. 行列で表す でオイラー陽解法を行列表記しました。その行列 $AA$ は以下でした。

$$A=(abcd)=(1−glΔt11) A = \begin{pmatrix} a & b \\ c & d \\ \end{pmatrix} = \begin{pmatrix} 1 & -\frac{g}{l}\Delta t \\ 1 & 1 \\ \end{pmatrix} \\$$

行列式 $∣A∣|A|$ を求めてみます。

$$∣A∣=∣1−glΔt11∣=1⋅1−(−glΔt)⋅1=1+glΔt \begin{aligned} |A| &= \begin{vmatrix} 1 & -\frac{g}{l}\Delta t \\ 1 & 1 \\ \end{vmatrix} \\ &= 1\cdot 1 - (- \frac{g}{l}\Delta t) \cdot 1 \\ &= 1 + \frac{g}{l}\Delta t \end{aligned}$$

$gg$ も $ll$ も非負なので $∣A∣=1+glΔt≧1|A| = 1 + \frac{g}{l}\Delta t \geqq 1$ で、変換のたびに増大していくことが分かります。

4.4. シンプレクティック法の拡大率を見てみよう

まだこちらは行列の形式にしていないのでまずは行列の形にしていきます。

$$v(t+Δt)=v(t)+(−glx(t))Δtx(t+Δt)=x(t)+v(t+Δt)Δt v(t+\Delta t) = v(t) + (- \frac{g}{l}x(t))\Delta t \\ x(t+\Delta t) = x(t) + v(t+\Delta t)\Delta t$$

座標 $x(t+Δt)x(t+\Delta t)$ の式中の $v(t+Δt)v(t+\Delta t)$ を代入して具体的に書いていきます。

$$x(t+Δt)=x(t)+v(t+Δt)Δt=x(t)+(v(t)+(−glx(t))Δt)Δt=x(t)+v(t)Δt−gΔt2lx(t)=Δtv(t)+(1−gΔt2l)x(t) \begin{aligned} x(t+\Delta t) &= x(t) + v(t+\Delta t) \Delta t \\ &= x(t) + ( v(t) + (- \frac{g}{l}x(t))\Delta t )\Delta t \\ &= x(t) + v(t)\Delta t - \frac{g\Delta t^2}{l}x(t) \\ &= \Delta t v(t) + (1 - \frac{g\Delta t^2}{l}) x(t) \\ \end{aligned}$$

では、オイラー陽解法と同様に行列形式で書いてみます。

$$(v(t+Δt)x(t+Δt))=(abcd)(v(t)x(t))=(av(t)+bx(t)cv(t)+dx(t)) \begin{aligned} \begin{pmatrix} v(t+\Delta t) \\ x(t+\Delta t) \\ \end{pmatrix} &= \begin{pmatrix} a & b \\ c & d \\ \end{pmatrix} \begin{pmatrix} v(t) \\ x(t) \\ \end{pmatrix} \\ &= \begin{pmatrix} av(t) + bx(t) \\ cv(t) + dx(t) \\ \end{pmatrix} \end{aligned}$$

$$(v(t+Δt)x(t+Δt))=((1)⋅v(t)+(−lgΔt)⋅x(t)(Δt)⋅v(t)+(1−gΔt2l)⋅x(t)) \begin{pmatrix} v(t+\Delta t) \\ x(t+\Delta t) \\ \end{pmatrix} = \begin{pmatrix} (1)\cdot v(t) + (-\frac{l}{g}\Delta t)\cdot x(t) \\ (\Delta t)\cdot v(t) + (1 - \frac{g\Delta t^2}{l})\cdot x(t) \\ \end{pmatrix}$$

よって

$$A=(abcd)=(1−glΔtΔt1−glΔt2) A = \begin{pmatrix} a & b \\ c & d \\ \end{pmatrix} = \begin{pmatrix} 1 & -\frac{g}{l}\Delta t \\ \Delta t & 1 - \frac{g}{l}\Delta t^2 \\ \end{pmatrix} \\$$

行列式 $∣A∣|A|$ を求めます。

$$∣A∣=∣1−glΔtΔt1−glΔt2∣=1⋅(1−glΔt2)−(−glΔt)Δt=1−glΔt2+glΔt2=1 \begin{aligned} |A| &= \begin{vmatrix} 1 & -\frac{g}{l}\Delta t \\ \Delta t & 1 - \frac{g}{l}\Delta t^2 \\ \end{vmatrix} \\ &= 1\cdot (1 - \frac{g}{l}\Delta t^2) - (-\frac{g}{l}\Delta t)\Delta t \\ &= 1 - \frac{g}{l}\Delta t^2 + \frac{g}{l}\Delta t^2 \\ &= 1 \end{aligned}$$

シンプレクティック法では行列式 $∣A∣|A|$ が $11$ となりました！！
よって、行列によって変換はされてもその拡大率は1でエネルギーが保存され、発散などが生じないと言えます。

雑な絵で厳密ではないですが、拡大率が1の場合は斜線部分の面積が $±0\pm 0$ で周期的に繰り返しても総和が0のままで、一方で拡大率が1以上の場合は赤線のように徐々に拡大して発散してしまうイメージです。

5. まとめ

オイラー陽解法の式を少しだけ変更するだけで誤差がほとんど拡大していかないシンプレクティック法について説明しました。

おわりに

たったこれだけの変更で明らかに差が出るのは面白いですね。
ただなぜそうなるのかなど理論の方はだいぶ数学の知識が必要で難しいです。なぜ発散しないのかについては理解が浅い気がしたので書くか迷いましたが助言など頂ければなと思い書きました。調べている中でだいぶ数学の知識を忘れていたり、そもそも履修していない内容だったりで改めて数学力が欲しいなと感じました。また学んだことなど書いていきたいと思います。