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桁で計算してますよ。
使い方は、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くらいでもう終わってるんだけど、、、
« Scratchで新しい三体問題の周期解（1223個）を計算してみる、、、計算できたが遅すぎる。8の字解を計算したものを代わりに。。。 | トップページ | Apple Watch Series3(LTE対応）がiFixitで早くも分解。LTE部分はSiPの後ろに並べただけ、、、ちょっとカッコ悪い。 »
「学問・資格」カテゴリの記事
- アメリカの空港で5G C-bandをを使うなという話で揉めてますが、どう考えても電波高度計(4.2-4.4GHz）が無茶苦茶で、衛星通信(3.7-4.2GHz)を邪魔しているのにC-Band(3.7–3.98GHz)が入ってくるとクレームとは、、、ジャイアニズムか！(2022.01.22)
- 不定方程式a*x+b*y=gcd(a,b) （ディオファントス方程式、ベズーの等式）の拡張ユークリッド互除法を使った計算をカシオの高精度計算サイトkeisan.casio.jpにUP!2022年共通テストの数学I・数学Aの問題4に出てても結果に自信ないので機械に頼る。。。(2022.01.19)
- 2022年度共通テストの数学I・数学Aの問題、例の防衛省のイージス・アシュアの縮尺違い地図で角度間違ってる話をそのまま持ってきてる！秋田県男鹿市の本山か。。。(2022.01.17)
- 高周波（RF・マイクロ波・ミリ波・5G）関連ニュース2022年1月18日 IEEE Microwave Magazineは女性高周波研究者特集、Microwave Journalはレーダで月面が見える。ドコモの誘電体線路アンテナ、3.45GHzのオークション結果、hBNを使ったスイッチ、Scikit-rfの使い方など。(2022.01.18)
- 犬や猫などのペットにマイクロチップを埋め込むことが義務化されるが、どんな無線通信技術を使っているか調べた。ISO11784/11785に準拠した134.2kHzを使うRFIDで、フェライトロッドアンテナ、キャパシタ、チップ内蔵。15桁の個体識別番号が読み取れる。(2022.01.10)
「日記・コラム・つぶやき」カテゴリの記事
- IMO(国際数学オリンピック）1998年問題6、というのがメールマガジンで回ってきた。a²+b² が ab +1で割れるとき、a²+b² / ab +1が完全平方数になるというもの。Pythonで何も考えず総当たりで計算する。(2022.01.28)
- 共通テスト試験問題流出で、カメラとSIMカード（セルラー通信用）使えるスマートウォッチ使ったら簡単にできてしまうのでは、と思った。(2022.01.27)
- PREMによる地球内部・外部の重力加速度gの計算をカシオの高精度計算サイトkeisan.casio.jpに自作式としてＵＰ！地球内部のほうが表面よりｇが大きいところがあるのが面白いです。(2022.01.27)
- BALMUDA Phoneの分解見た。近年まれに見るメインボードのスカスカさ。Xiaomi Mi 11、Oppo Reno 6 pro、Google Pixel 6 pro、iPhone 13 Proと比べよう。また中華圏のトレンドのバッテリーをフレキでまたぐ構造もしてない。(2022.01.26)
- PREMによる地球内部の密度分布をカシオの高精度計算サイトkeisan.casio.jpにUP!(2022.01.25)
« Scratchで新しい三体問題の周期解（1223個）を計算してみる、、、計算できたが遅すぎる。8の字解を計算したものを代わりに。。。 | トップページ | Apple Watch Series3(LTE対応）がiFixitで早くも分解。LTE部分はSiPの後ろに並べただけ、、、ちょっとカッコ悪い。 »
コメント