Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その6:ルンゲクッタ8次 DOP853
今回は第6回。ルンゲクッタ8次のDormand & Prince法の有名なFortranのルーチンをVBAに移植しています。
使い方は
DOP853Lib.bas
DriverDOP853Lib.bas
という2つのファイルをインポートします。
DOP853Libのほうが本体で、DriverDOP853Libというのはドライバで問題によって書き換えることになります。
関数はドライバの中のFCNに書きます。
例えばローレンツ方程式なら、
Sub FCN(N As Long, x As Double, y() As Double, F() As Double, RPAR() As Double, IPAR() As Long)
Dim sigma As Double, r As Double, b As Double
sigma = RPAR(1)
r = RPAR(2)
b = RPAR(3)
F(1) = -sigma * (y(1) - y(2))
F(2) = -y(2) - y(1) * y(3) + r * y(1)
F(3) = y(1) * y(2) - b * y(3)
End Sub
密出力ルーチンはSOLOUTで、ここを書き換えて自由なところに出力できます。
Sub SOLOUT(NR As Long, XOLD As Double, x As Double, y() As Double, _
N As Long, CON() As Double, ICOMP() As Long, ND As Long, _
RPAR() As Double, IPAR() As Long, IRTRN As Long, XOUT As Double)
' --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS
' --- BY USING "CONTD8", THE CONTINUOUS COLLOCATION SOLUTION
Dim K As Long
If (NR = 1) Then
Worksheets("Sheet2").Cells(JJ + 1, 3) = x
For K = 1 To N
Worksheets("Sheet2").Cells(JJ + 1, 3 + K) = y(K)
Next K
XOUT = HH
JJ = JJ + 1
Else
Label10:
If (x >= XOUT) Then
Worksheets("Sheet2").Cells(JJ + 1, 3) = XOUT
For K = 1 To N
Worksheets("Sheet2").Cells(JJ + 1, 3 + K) = CONTD8(K, XOUT, CON, ICOMP, ND)
Next K
XOUT = CDbl(JJ) * HH
JJ = JJ + 1
GoTo Label10
End If
End If
End Sub
ドライバ本体の誤差に関わるパラメータは、本家のDOP853の説明を参照。
ライブラリ本体:
またルンゲクッタ8次のDOP853ルーチンとそのドライバ。
メルセンヌツイスタ用
« #仮面ライダービルド 第22話の話数を表す数式はまたもやラマヌジャン!2143/π^4≃22、あるいはπ≃(9^2+19^2/22)^(1/4) | トップページ | Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その7:メルセンヌツイスタ »
「パソコン・インターネット」カテゴリの記事
- 高周波回路シミュレータQucsStudioを使ってみる(その3)Mixed Mode S parameterを計算(2018.10.12)
- 高周波回路シミュレータQucsStudioを使ってみる(その2)SパラメータのTouchStoneフォーマットで出力するには?(2018.10.11)
- 高周波回路シミュレータQucsStudioを使ってみる(その1)まずは何をさておきμの文字化けだけには注意。(2018.10.10)
- 円の弧長,弦長,矢高,半径のどれか2つを与えて残りを計算(カシオの高精度計算サイト自作式)で180°以上、複数解に対応。(2018.10.09)
- macbook proをmacOS Mojaveにアップデート。せっかくなんでダークモードにしてみる。(2018.09.26)
「学問・資格」カテゴリの記事
- iPhoneで使える数学ソフトMaple Calculatorで手書き入力でMathematicaの常微分方程式の例を試す(その1)(2021.03.29)
- 高周波(RF・マイクロ波・ミリ波・5G)関連ニュース2021年3月16日 IEEE Microwave magazineの特集は中国のマイクロ波技術、Microwave Journalの特集はCATRでミリ波ビームスキャンを測定、SIJの特集はなんと450ドルのネットアナ、フリーの高周波シミュレータ(Qucsなど)。(2021.03.16)
- 高周波(RF・マイクロ波・ミリ波・5G)関連ニュース2021年2月16日 IEEE Microwave Magazineの特集はオールデジタルのRFID、Microwave JournalはEバンド ミリ波通信に衛星や気球を使う話、アメリカの半導体企業がバイデンに投資を迫る、(2021.02.17)
- カオスを生じる電気回路、Chua’s circuitをLTspiceで回路シミュレーションしてみる。(2021.02.19)
- Labyrinth Chaos(迷宮カオス)を生むThomas-Rössler方程式のパラメータbを色々変えて、Python+Scipyでルンゲクッタ8次のDOP853(Dormand&Prince)を使って計算してGIFアニメ(2021.02.16)
「日記・コラム・つぶやき」カテゴリの記事
- 新型コロナウイルス、中国、日本、韓国、アメリカ、ドイツ、フランス、イギリスでの感染者数を指数関数&ロジスティック関数&Log-Logプロットでべき関数フィッティングした(4/11更新)日本は陽性者50万人を超え第四波の始まり、というかイギリス以外は全部増えだしてる。(2021.04.12)
- 3乗すると(奇数回べき乗すると)元に戻る数59161727622001114846846461792218008213239954784512519836425781249(2乗すると元に戻る数はペレリマン数列)(2021.04.11)
- 円と接線を使った1枚の絵でsin/cos/tan/sec/csc/cotを表すもの、複素数を使って接線のパラメータ表示(z=exp(iθ)+i exp(iθ) t )で計算してみる。(2021.04.05)
- 新型コロナウイルス、中国、日本、韓国、アメリカ、ドイツ、フランス、イギリスでの感染者数を指数関数&ロジスティック関数&Log-Logプロットでべき関数フィッティングした(4/4更新) フランス・ドイツの増え方がエグイが日本・韓国もやはり増加中。イギリスの減速がすごいな。(2021.04.05)
- 分数階微分(Fractional Calculus)の計算をPythonのdifferint(Riemann-Liouvilleの定義で計算するもの)でやってみる。(2021.04.02)
トラックバック
この記事へのトラックバック一覧です: Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その6:ルンゲクッタ8次 DOP853:
« #仮面ライダービルド 第22話の話数を表す数式はまたもやラマヌジャン!2143/π^4≃22、あるいはπ≃(9^2+19^2/22)^(1/4) | トップページ | Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その7:メルセンヌツイスタ »
コメント