« 酔わないウメッシュのCMに出てるのは北乃きいちゃん。 | トップページ | ちふれのCM”SAVE WOMAN”に出てくる女の子がとんでもなく可愛いです(木村文乃さん)/ドコモスマートフォンCMで桑田さんとも共演 »

2011年4月18日 (月)

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

「DOP853.bas」をダウンロード

例としてローレンツ方程式を解いたドライバは、

「Lorenz.bas」をダウンロード

内容の説明は明日やります。

ライセンスは

「licence.txt」をダウンロード

に従います。

« 酔わないウメッシュのCMに出てるのは北乃きいちゃん。 | トップページ | ちふれのCM”SAVE WOMAN”に出てくる女の子がとんでもなく可愛いです(木村文乃さん)/ドコモスマートフォンCMで桑田さんとも共演 »

学問・資格」カテゴリの記事

コメント

このDOP853VBA移植版も含めてですが、ここのブログにおいてあるコードを利用して別のコードを書き、配布してもよろしいのでしょうか?

オリジナルのDOP853のライセンス条件は、
http://www.unige.ch/~hairer/prog/licence.txt
です。本ブログのソースもそれに従います。
なので、このテキストを一緒に配布すればOKです。

ありがとうございます。

もしかしてこのVBA版をさらにC#に移植されようとしている方でしょうか?
実はC版、C++版を別の方が作られています。
そっちの方を移植された方が速いのでは?
C版:
http://www.unige.ch/~hairer/prog/nonstiff/cprog.tar
C++版:
http://www.unige.ch/~hairer/prog/IntegratorT.tgz

あ、はいそうです汗
C版/C++版のご紹介ありがとうございます。

C#が使い慣れてたので移植させていただこうかと思いました。一応機械的に変換して移植は済みました。

それと、VBA版コード中の

FACOLD = Max(Err, 0.0001)

FACOLD = Max(ERR1, 0.0001)

の誤植ではありませんか?ERR1はFORTRANの元コードではERRとなっていますが…どうなのでしょう?無駄な指摘だったら申し訳ありません。

あ、そうですね。私のバグです。
ErrってのはVBAでは予約語(エラーの状態のオブジェクト)なので、使えなくて
ERR1に書き直したはずが、抜け落ちてました…
予約語なのでエラーにならなかった…
ぼつぼつ直してアップロードしなおします。

あと、、、
BZ反応をやられているならぜひ!スパイラル
パターンを観察してください!
http://sci.tea-nifty.com/blog/2011/01/excel-vbacomple.html
こういうパターンが得られます。

>スパイラルパターン
おおお!すごいですね!

ただ、私は現在高2で、冬になったら受験モードに切り替えなければならないので時間がそれを許してくれるかどうかわからないです…でもできる限りやってみます。

少し質問させていただきたいのですが、古典的RKでOregonatorを計算してみたところ、Scholarpediaなどに載っているグラフが再現しませんでした。これはやはり誤差が原因でしょうか?

高二でこんなことに興味を持つなんてすごいね!将来が楽しみです。
で、Oregonatorっていうのは硬い(Stiffな)微分方程式として知られていて、
パラメータによってはめちゃくちゃ時間刻み幅を小さくしないとたぶん、うまく
計算できません。刻み幅をびっくりするほど小さくしてみてはどうでしょう。
時間かかりますが。
実はStiffなものにはStiff専用の解き方があって(陰解法)、これはまた大学で
勉強してもらえればと。

あと、BZのスパイラルパターンの実験は、
http://www.youtube.com/watch?v=3JAqrRnKFHo&feature=related
なんかがあって、ちょっと途中でクリップでいじったりするとターゲットパターン
(同心円)がスパイラルになったりします。こういうのっておもしろいね。

>刻み幅をびっくりするほど小さくしてみてはどうでしょう。
時間かかりますが。

0.001とか0.0001くらいですかね?試してみます。

>実はStiffなものにはStiff専用の解き方があって(陰解法)、これはまた大学で
勉強してもらえればと。

関連ページを覗いて見たのですが難しくてサッパリでした…

>ちょっと途中でクリップでいじったりするとターゲットパターン
(同心円)がスパイラルになったりします。

おもしろいですよね~。この化学波、普通の波ともまた違う動きをするのがすごく面白いです。
一回ペトリ皿で空間振動の実験をしてみましたが、スパイラルはうまくつくれずにぐちゃぐちゃになってしまいました…

現在は一応、反応系の酸化還元電位を長期的に観測して、それをシミュレーションと結びつけることをメインに置いて研究を進めています。どちらかというと時間振動に重きを置いていますが…

とにかく、時間の許す限りいろいろと頑張ってみようと思います。
ありがとうございました。

コメントを書く

(ウェブ上には掲載しません)

トラックバック

この記事のトラックバックURL:
http://app.cocolog-nifty.com/t/trackback/512682/51389502

この記事へのトラックバック一覧です: Dormand-Prince(ルンゲ・クッタ8次)の有名なルーチンDOP853をExcel VBAに移植:

« 酔わないウメッシュのCMに出てるのは北乃きいちゃん。 | トップページ | ちふれのCM”SAVE WOMAN”に出てくる女の子がとんでもなく可愛いです(木村文乃さん)/ドコモスマートフォンCMで桑田さんとも共演 »

最近のコメント

2017年11月
      1 2 3 4
5 6 7 8 9 10 11
12 13 14 15 16 17 18
19 20 21 22 23 24 25
26 27 28 29 30    
フォト
無料ブログはココログ