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で桑田さんとも共演 »
「学問・資格」カテゴリの記事
- 高周波・RFニュース 2024年12月6日 NGMNが無線パフォーマンス評価フレームワーク発行、5GAAがC-V2Xのロードマップ発行、Marvellの3nm 1.6Tbps PAM4インターコネクト、Nokiaの2.4Tbps光伝送、Silicon Labsの低消費電力モジュール、Xiaomi 14T Pro分解動画(2024.12.06)
- 高周波回路シミュレータQucsStudioがuSimmicsに名称変更し、バージョンも4.8.3から5.8にアップデートされた。Qucsと区別するためだそうだ。また、Pythonの高周波用ライブラリscikit-rfもv1.5.0にバージョンアップされていた(2024.12.04)
- 日経サイエンス2025年1月号の特集 和算再発見の佐藤賢一さんの記事「算聖 関孝和の実像」に出てきた矢高に対する円弧の２乗の近似式をカシオの高精度計算サイトkeisan.casio.jpの自作式として作った。ものすごい精度であることがよくわかる。(2024.12.03)
- MATLAB Onlineで高周波用のRF Toolboxを使ってみる。Touchstoneファイルの読み込み、dB表示グラフ、スミスチャートなど簡単にできるし、フィルタ合成やIEEE P370 De-embedding(ZC-2xThru)も使える（MATLABで書かれたものがオリジナル）。(2024.12.05)
コメント
« 酔わないウメッシュの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分