« Scratchで新しい三体問題の周期解(1223個)を計算してみる、、、計算できたが遅すぎる。8の字解を計算したものを代わりに。。。 | トップページ | Apple Watch Series3(LTE対応)がiFixitで早くも分解。LTE部分はSiPの後ろに並べただけ、、、ちょっとカッコ悪い。 »

2017年9月25日 (月)

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

Lorenz14thparigp

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

で、このブログではルンゲクッタ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次の比較。

Lorenzcmp01

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

Lorenzcmp02_2

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

Lorenzcmp03

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

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

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

Lorenzcmp04_2

« Scratchで新しい三体問題の周期解(1223個)を計算してみる、、、計算できたが遅すぎる。8の字解を計算したものを代わりに。。。 | トップページ | Apple Watch Series3(LTE対応)がiFixitで早くも分解。LTE部分はSiPの後ろに並べただけ、、、ちょっとカッコ悪い。 »

学問・資格」カテゴリの記事

日記・コラム・つぶやき」カテゴリの記事

コメント

コメントを書く

(ウェブ上には掲載しません)

トラックバック

この記事のトラックバックURL:
http://app.cocolog-nifty.com/t/trackback/512682/65836680

この記事へのトラックバック一覧です: 35段14次のルンゲクッタ法をPARI/GPに実装、ローレンツ方程式を計算し、通常の4次、そして8次、オイラー法と比較してみる。(再掲):

« Scratchで新しい三体問題の周期解(1223個)を計算してみる、、、計算できたが遅すぎる。8の字解を計算したものを代わりに。。。 | トップページ | Apple Watch Series3(LTE対応)がiFixitで早くも分解。LTE部分はSiPの後ろに並べただけ、、、ちょっとカッコ悪い。 »

最近のコメント

2017年11月
      1 2 3 4
5 6 7 8 9 10 11
12 13 14 15 16 17 18
19 20 21 22 23 24 25
26 27 28 29 30    
フォト
無料ブログはココログ