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で桑田さんとも共演 »
「学問・資格」カテゴリの記事
- Qwen3.6-35B-A3Bが発表され、Ollamaでも使える。そこで電子レンジの動作原理（2.45GHzは水分子の共振周波数でない）と隕石が大気圏突入で燃える原理（摩擦熱ではない）を聞くと、誘電緩和と断熱圧縮について正しく答えられた。今までのローカルLLMで一番賢い回答と思う。(2026.04.17)
- 高周波・RFニュース 2026年4月17日 atisの3GPP Rel.20ウェビナー動画公開、MWCバルセロナ2026でのGSMA Device Enablement Summit資料公開、ハリファ大学が無線周波数AI言語モデルRF-GPT発表、レドームの解説など(2026.04.17)
- ExcelのOfficeスクリプト(TypeScript)で数値計算ライブラリmath.jsを使う(1) Officeスクリプトは外部API呼び出せるし、math.jsは RESTful APIで呼び出せることがわかった。まずは選択したセルのデータを読み、行列演算。LU分解で一次方程式を解き、逆行列と行列式を求める。(2026.04.17)
- 高周波・RFニュース 2026年4月16日 AmazonがGlobalstarを買収、GSMAが日本のデジタル化をレポート、Mini-Circuitsがケーブルアセンブリを動画で解説、Kymetaが米国海軍研究局と衛星通信で契約、PerasoがドローンIFF向け60GHzモジュール出荷、SEMCOが1500V耐圧MLCC発表(2026.04.16)
- 高周波・RFニュース 2026年4月15日 Microwave Journalはアンプと発振器特集、Signal Integrity Journalは100GHz越えのインターコネクトのAIを使うHFSSモデル化、ローデ・シュワルツが潜水艦通信をUDT2026で発表、Xiaomi Poco X8 Pro分解動画、atisの5Gポリシーレポート(2026.04.15)
コメント
« 酔わないウメッシュの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分