Dormand-Prince(ルンゲ・クッタ8次）の有名なルーチンDOP853をExcel VBAに移植
このブログも3周年ということで、前々から気になっていたことを修正したい。今までうちのブログでRunge-Kutta8次のドルマン・プリンス法でいろいろ計算してきた（と言っていた）けれど、実は使っていた公式は8次は8次だけれど別の人が作った公式だったorz（コメントいただいた方に指摘されて見直したら間違っていたことに気付いた。ありがとうございます。信頼できないサイトを引用したのがミス。。。）。
そこで、ちゃんと本物のDormand & Prince 8次の公式に基づいて計算しなおそう。もっとも信頼できるのはこのサイトにある
http://www.unige.ch/~hairer/software.html
DOP853ルーチン。ちゃんと適応刻み幅と密出力ができるものとなっている。ただしFORTRAN。
そこでFORTRANからExcel VBAに移植してみた。
メインルーチンは、このDOP853.bas
例としてローレンツ方程式を解いたドライバは、
内容の説明は明日やります。
ライセンスは
「licence.txt」をダウンロード
に従います。
« 酔わないウメッシュのCMに出てるのは北乃きいちゃん。 | トップページ | ちふれのCM”SAVE WOMAN”に出てくる女の子がとんでもなく可愛いです（木村文乃さん）/ドコモスマートフォンCMで桑田さんとも共演 »
「学問・資格」カテゴリの記事
- 「数学がゲームを動かす！ ゲームデザインから人工知能まで」を読んだ。面白い！パックマンのアルゴリズムやドラクエの計算式、ドンキーコングはベルレ法でジャンプ、カルマンフィルタ、遺伝的アルゴリズム、セガの線形代数本を書かれた方は理論物理出身など話題が豊富。(2025.05.13)
- 高周波・RFニュース 2025年5月12日 IEEE Microwave Magazineでミリ波ガラス基板・超伝導・液晶などの記事、Pythonの高周波ライブラリscikit-rfがv1.7.0に、Megamagneticsの希土類を使わないミリ波サーキュレータ、CoilcraftのRFインダクタホワイトペーパー(2025.05.12)
- 高周波・RFニュース 2025年5月9日 QorvoがスマートファクトリーにUWB利用の解説、KYOCERA AVXがUHF,VHF向け高方向性カプラ発表、Wireless Broadband Allianceが企業向けWi-Fi 7トライアル結果報告、NXPが第三世代イメージングレーダプロセッサ発表(2025.05.09)
- 高周波・RFニュース 2025年5月8日 6Gにはテラヘルツどころか7-14GHzも向いてないという意見、フィンランド企業が6G推進コンソーシアムRF ECO3発表、Infineonが7P7T内蔵７ch-LNA発表、SpirentがOctoboxにWi-Fi 6/7自動テスト追加(2025.05.08)
コメント
« 酔わないウメッシュのCMに出てるのは北乃きいちゃん。 | トップページ | ちふれのCM”SAVE WOMAN”に出てくる女の子がとんでもなく可愛いです（木村文乃さん）/ドコモスマートフォンCMで桑田さんとも共演 »
このDOP853VBA移植版も含めてですが、ここのブログにおいてあるコードを利用して別のコードを書き、配布してもよろしいのでしょうか？
投稿: brv314 | 2011年9月25日 (日) 23時27分
オリジナルのDOP853のライセンス条件は、
http://www.unige.ch/~hairer/prog/licence.txt
です。本ブログのソースもそれに従います。
なので、このテキストを一緒に配布すればOKです。
投稿: tonagai | 2011年9月26日 (月) 20時51分
ありがとうございます。
投稿: brv314 | 2011年9月26日 (月) 22時20分
もしかしてこのVBA版をさらにC#に移植されようとしている方でしょうか？
実はC版、C++版を別の方が作られています。
そっちの方を移植された方が速いのでは？
C版：
http://www.unige.ch/~hairer/prog/nonstiff/cprog.tar
C++版：
http://www.unige.ch/~hairer/prog/IntegratorT.tgz
投稿: tonagai | 2011年9月27日 (火) 07時16分
あ、はいそうです汗
C版/C++版のご紹介ありがとうございます。
C＃が使い慣れてたので移植させていただこうかと思いました。一応機械的に変換して移植は済みました。
それと、VBA版コード中の
FACOLD = Max(Err, 0.0001)
は
FACOLD = Max(ERR1, 0.0001)
の誤植ではありませんか？ERR1はFORTRANの元コードではERRとなっていますが…どうなのでしょう？無駄な指摘だったら申し訳ありません。
投稿: brv314 | 2011年9月27日 (火) 22時43分
あ、そうですね。私のバグです。
ErrってのはVBAでは予約語（エラーの状態のオブジェクト）なので、使えなくて
ERR1に書き直したはずが、抜け落ちてました…
予約語なのでエラーにならなかった…
ぼつぼつ直してアップロードしなおします。
あと、、、
BZ反応をやられているならぜひ！スパイラル
パターンを観察してください！
https://sci.tea-nifty.com/blog/2011/01/excel-vbacomple.html
こういうパターンが得られます。
投稿: tonagai | 2011年9月27日 (火) 22時53分
＞スパイラルパターン
おおお！すごいですね！
ただ、私は現在高2で、冬になったら受験モードに切り替えなければならないので時間がそれを許してくれるかどうかわからないです…でもできる限りやってみます。
少し質問させていただきたいのですが、古典的RKでOregonatorを計算してみたところ、Scholarpediaなどに載っているグラフが再現しませんでした。これはやはり誤差が原因でしょうか？
投稿: brv314 | 2011年9月27日 (火) 23時12分
高二でこんなことに興味を持つなんてすごいね！将来が楽しみです。
で、Oregonatorっていうのは硬い(Stiffな）微分方程式として知られていて、
パラメータによってはめちゃくちゃ時間刻み幅を小さくしないとたぶん、うまく
計算できません。刻み幅をびっくりするほど小さくしてみてはどうでしょう。
時間かかりますが。
実はStiffなものにはStiff専用の解き方があって（陰解法）、これはまた大学で
勉強してもらえればと。
あと、BZのスパイラルパターンの実験は、
http://www.youtube.com/watch?v=3JAqrRnKFHo&feature=related
なんかがあって、ちょっと途中でクリップでいじったりするとターゲットパターン
（同心円）がスパイラルになったりします。こういうのっておもしろいね。
投稿: tonagai | 2011年9月27日 (火) 23時27分
＞刻み幅をびっくりするほど小さくしてみてはどうでしょう。
時間かかりますが。
0.001とか0.0001くらいですかね？試してみます。
＞実はStiffなものにはStiff専用の解き方があって（陰解法）、これはまた大学で
勉強してもらえればと。
関連ページを覗いて見たのですが難しくてサッパリでした…
＞ちょっと途中でクリップでいじったりするとターゲットパターン
（同心円）がスパイラルになったりします。
おもしろいですよね～。この化学波、普通の波ともまた違う動きをするのがすごく面白いです。
一回ペトリ皿で空間振動の実験をしてみましたが、スパイラルはうまくつくれずにぐちゃぐちゃになってしまいました…
現在は一応、反応系の酸化還元電位を長期的に観測して、それをシミュレーションと結びつけることをメインに置いて研究を進めています。どちらかというと時間振動に重きを置いていますが…
とにかく、時間の許す限りいろいろと頑張ってみようと思います。
ありがとうございました。
投稿: brv314 | 2011年9月27日 (火) 23時58分