« #仮面ライダービルド 第22話の話数を表す数式はまたもやラマヌジャン!2143/π^4≃22、あるいはπ≃(9^2+19^2/22)^(1/4) | トップページ | Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その7:メルセンヌツイスタ »

2018年2月12日 (月)

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:メルセンヌツイスタ »

パソコン・インターネット」カテゴリの記事

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

日記・コラム・つぶやき」カテゴリの記事

コメント

コメントを書く

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

トラックバック

« #仮面ライダービルド 第22話の話数を表す数式はまたもやラマヌジャン!2143/π^4≃22、あるいはπ≃(9^2+19^2/22)^(1/4) | トップページ | Excel VBAで複素数演算(一次方程式・FFT他, Numerical Recipes移植)、フィッティング(非線形含む)、ルンゲクッタ8次(DOP853)などが使えるライブラリ その7:メルセンヌツイスタ »

最近のコメント

2018年9月
            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            
フォト
無料ブログはココログ