35段14次のルンゲクッタ法をPARI/GPに実装、ローレンツ方程式を計算し、通常の4次、そして8次、オイラー法と比較してみる。（再掲）

ルンゲクッタ8次のDOP853で三体問題の周期解を1223個見つけた話をこの前見て検証してみたりしてましたが、

14次の公式があるのを知ってますか？

これだ！

http://sce.uhcl.edu/rungekutta/

35段14次！

これは倍精度くらいではそもそも使う意味がなさそうなので、任意精度の計算ができるPARI/GPに実装してみる。

まずはいつものようにローレンツ方程式をやってみよう。式はこちら。

dx/dt = -σ*(x-y)

dy/dt=-y-x*z+r*x

dz/dt=x*y-b*z

でプログラムのソースコードはこちら。60桁で計算してますよ。

「Lorenz.gp」をダウンロード

使い方は、PARI/GPを立ち上げて、

\r Lorenz

でロード。

rk14("ファイル名")

で実行されて、ファイル名のところに保存される。散布図がどうも描きにくそうなので、ファイルに落としてExcelで描いてみた。それがこれ。

なるほど、計算できてそうだ。

で、このブログではルンゲクッタ8次のDormand & Princeルーチン（DOP853)をExcel VBAに移植したりまた普通のルンゲクッタ4段4次をカシオの高精度計算サイトkeisan.casio.jpにUPしたり、結構いろいろやったので全部比較してみよう。

お題は、やはりローレンツ方程式でσ=10,r=26, b=8/3, 時間刻み幅0.01, そして初期値はx=y=z=0.1としてやってみる。

まずは4次、8次、14次の比較。

ｔ=25くらいですでに4次がやばくなってきた。次行こう。

もうt=40過ぎると8次もやばい。

t=100くらいまでいくともうぐだぐだ。

ということで、常微分方程式なら全部4段4次のRunge kuttaでいいや、とか思っているといろいろと落とし穴がでたり。注意しよう。

え？オイラー法？t=1くらいでもう終わってるんだけど、、、