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次の比較。
t=25くらいですでに4次がやばくなってきた。次行こう。
もうt=40過ぎると8次もやばい。
t=100くらいまでいくともうぐだぐだ。
ということで、常微分方程式なら全部4段4次のRunge kuttaでいいや、とか思っているといろいろと落とし穴がでたり。注意しよう。
え?オイラー法?t=1くらいでもう終わってるんだけど、、、
« Scratchで新しい三体問題の周期解(1223個)を計算してみる、、、計算できたが遅すぎる。8の字解を計算したものを代わりに。。。 | トップページ | Apple Watch Series3(LTE対応)がiFixitで早くも分解。LTE部分はSiPの後ろに並べただけ、、、ちょっとカッコ悪い。 »
「学問・資格」カテゴリの記事
- 高周波・RFニュース 2025年7月17日 IEEE Microwave Magazineで地上監視レーダ解説など、Microwave Journalでミリ波AiPのOTA測定解説など、5G Americasが5G Advancedホワイトペーパー発行、TDKのAI 微細欠陥検出edgeRX Vision、IEEE Trans. on THz の今月号発行(2025.07.17)
- Xのタイムラインが何かランダムウォークのガウス分布じゃなくて逆正弦法則のように見えてきている…仕方ないのでJulia言語でモンテカルロシミュレーションで試行回数を変えて逆正弦法則とガウス分布を見る。(2025.07.16)
- 高周波・RFニュース 2025年7月13日 Pythonの高周波ライブラリscikit-rfがv1.8.0に、SamsungがGalaxy Z Fold7など発表→QualcommがSnapdragon 8 Eliteが使われていると発表、NGMNが基地局アンテナの推奨事項をまとめる、STMicroとMetalenzがメタサーフェス光学のライセンス締結(2025.07.14)
- Interface2025年8月号Pythonで体験!はじめての暗号を買った。上杉暗号からRSA、AES、DHなど、特に楕円曲線暗号についてはコードも実際に動かすところまで詳しくかかれていた。耐量子暗号や聞いたことなかったY-00暗号や関数型暗号も記載。(2025.07.10)
- 高周波・RFニュース 2025年7月8日 NordicとSercommのセルラーIoTモジュール、iFixitがFairphone 6を分解、スコアは10/10、RCR wireless newsのウェビナー2件(6GとIndustry4.0)、SEMCOが高耐圧C0G MLCCを車載急速充電に提案(2025.07.09)
「日記・コラム・つぶやき」カテゴリの記事
- 高周波・RFニュース 2025年7月17日 IEEE Microwave Magazineで地上監視レーダ解説など、Microwave Journalでミリ波AiPのOTA測定解説など、5G Americasが5G Advancedホワイトペーパー発行、TDKのAI 微細欠陥検出edgeRX Vision、IEEE Trans. on THz の今月号発行(2025.07.17)
- Xのタイムラインが何かランダムウォークのガウス分布じゃなくて逆正弦法則のように見えてきている…仕方ないのでJulia言語でモンテカルロシミュレーションで試行回数を変えて逆正弦法則とガウス分布を見る。(2025.07.16)
- 高周波・RFニュース 2025年7月13日 Pythonの高周波ライブラリscikit-rfがv1.8.0に、SamsungがGalaxy Z Fold7など発表→QualcommがSnapdragon 8 Eliteが使われていると発表、NGMNが基地局アンテナの推奨事項をまとめる、STMicroとMetalenzがメタサーフェス光学のライセンス締結(2025.07.14)
- 高周波・RFニュース 2025年7月8日 NordicとSercommのセルラーIoTモジュール、iFixitがFairphone 6を分解、スコアは10/10、RCR wireless newsのウェビナー2件(6GとIndustry4.0)、SEMCOが高耐圧C0G MLCCを車載急速充電に提案(2025.07.09)
- 高周波・RFニュース 2025年7月2日 5G Americasが6Gに向けセンシングと通信ホワイトペーパー発行、KYOCERA AVXが3dBハイブリッドカプラ発表、TDKが車載薄膜インダクタ発表、Nordicが1次電池向けPMIC発表、ローデ・シュワルツの6GとAI/ML解説記事(2025.07.02)
« Scratchで新しい三体問題の周期解(1223個)を計算してみる、、、計算できたが遅すぎる。8の字解を計算したものを代わりに。。。 | トップページ | Apple Watch Series3(LTE対応)がiFixitで早くも分解。LTE部分はSiPの後ろに並べただけ、、、ちょっとカッコ悪い。 »
コメント