(cache)二重振り子の精度保証付き数値計算

2017/08/08(火)二重振り子の精度保証付き数値計算

少し前(5月頃)、二重振り子のシミュレーション動画がtwitterで流行したことがありました。最初のtweetがこれ。

わずかに異なる初期値をもった50本の二重振り子のシミュレーションを動画にしたもので、途中まで完全に重なっているように見える振り子が一気にバラける様子がとても面白い。

次は、それを三重振り子にしたもの。

gmpを用いて高精度計算をしており、ねとらぼで記事にもなりました。

どちらも常微分方程式の計算にはルンゲクッタ法を用いています。わずかな初期値のずれが後に大きな違いをもたらすことがとてもよく分かる動画ですが、一方で、ルンゲクッタ法で計算された値も真値とはわずかにずれており、当然その誤差も同様に後で大きな違いをもたらすことになります。だとすると、果たして意味のある計算になっているのかという疑問が生じます。

そこで、二重振り子の軌道を精度保証付きで計算し、ルンゲクッタ法とどのくらいずれるのか検証してみました。そのtweetがこれです。

このように、4次のルンゲクッタ法(刻み幅

2^{- 7}

)と精度保証付き計算とを比較してみると、

t = 8

あたりから先は完全に違う軌道になってしまっていることが分かりました。

以下では、twitterに書けなかった細かい情報を備忘録も兼ねて書いておこうと思います。

使った式は、

\begin{aligned} (m_{1} + m_{2}) l_{1} \ddot{θ_{1}} + m_{2} l_{2} \ddot{θ_{2}} \cos (θ_{1} - θ_{2}) + m_{2} l_{2} {\dot{θ_{2}}}^{2} \sin (θ_{1} - θ_{2}) + (m_{1} + m_{2}) g \sin θ_{1} = 0 \\ l_{1} l_{2} \ddot{θ_{1}} \cos (θ_{1} - θ_{2}) + l_{2}^{2} \ddot{θ_{2}} - l_{1} l_{2} {\dot{θ_{1}}}^{2} \sin (θ_{1} - θ_{2}) + g l_{2} \sin θ_{2} = 0 \end{aligned}

のようなものです(Wikipediaの記事と全く同じです)。舞台設定は

の通り。この式を

\ddot{θ_{1}}, \ddot{θ_{2}}

の2変数連立一次方程式と見て解いて正規形に直し、

\dot{θ_{1}} = ϕ_{1}, \dot{θ_{2}} = ϕ_{2}

のように置き直して1階4変数の常微分方程式とします。パラメータは

\begin{aligned} m_{1} & = m_{2} = 1 \\ l_{1} & = l_{2} = 1 \\ g & = 9.8 \end{aligned}

とし、初期値は

\begin{aligned} θ_{1} (0) & = \frac{3}{4} π \\ θ_{2} (0) & = \frac{3}{4} π \\ \dot{θ_{1}} (0) & = 0 \\ \dot{θ_{2}} (0) & = 0 \end{aligned}

としました。

kvライブラリを用いて、この初期値問題を精度保証付きで計算するプログラムと、4次のルンゲクッタ法で計算するプログラムを書きました。両方のプログラムをzipで上げておきます ( doublependulum.zip )。4次のルンゲクッタ法の刻み幅は固定で

2^{- 7}

としました。精度保証付きの方は全ての時刻の解を多項式の形で連続的に計算しているのですが、ルンゲクッタ法と同じ

2^{- 7}

の刻み幅で密出力させています。この計算結果を動画にしたものが、先のtweetというわけです。

プログラムを走らせて得られた生データを上げておきます ( doublependulum-result.zip )。

t = 1

のときの両プログラムのデータを抜き出すと、

	$θ_{1}$	$θ_{2}$
kv	[ 0.14046540255551162,0.14046540255558421]	[ -0.71766989300731332,-0.71766989300724726]
rk4	0.1404652582264235	-0.71766977132932186

となっていました。精度保証付きの方は全てのデータを区間で保持しており、この閉区間の中に真の解が存在することが数学的に保証されています。これを見ると、精度保証付き計算ではこの時点で13桁の精度を保っていますが、ルンゲクッタ法はすでに6桁程度しか精度がないことが分かります。例えば最初のtweetの動画では初期値に

10^{- 12}

の摂動を与えていますが、

t = 1

の時点ですでにルンゲクッタ法そのものの誤差の方が支配的になっているとするならば、本当に初期値の摂動の影響を計算していると言えるのでしょうか。

このデータをグラフにしてみました。

t = 8

付近でずれが目に見えるようになり、そのあとは完全に異なる軌道を描いています。

精度保証付き計算もいつまでも高精度を保っていられるわけではなく、

t = 17

付近を見るとわずかにグラフに幅があるのが分かります。この後一気に区間幅が爆発し、

t = 17.05

付近で計算不能に陥ります。この精度保証付き計算プログラムはgmpのような高精度数を使用しているわけではなく、内部で使っているのは普通のdoubleです。試しに内部の数値型をdd(擬似4倍精度数)に変えてみたところ、

t = 30

くらいまでは普通に破綻すること無く計算出来ました。

二重振り子の動画を見るだけなら、精度保証してもしなくても動きの「印象」はほとんど変わらないのですが、カオス系に対して数値シミュレーションで何らかの数学的な結論を得ようと思うならば、誤差が明確に把握できる精度保証付き数値計算は必須だと思います。みなさま、精度保証付き数値計算を使いましょう！