35段14次のルンゲクッタ法をPARI/GPで計算(Lorenz方程式)

先日導入したPARI/GPに35段14次のルンゲクッタ法を実装してみた。

係数はこちら↓60桁分あります。

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

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

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で描いてみた。それがこれ。

すばらしく安定感のある計算ができるなあ。これからはこれでしばらくいろいろ計算していってみよう（続く）。