Excel VBAで複素数演算（一次方程式・FFT他, Numerical Recipes移植）、フィッティング（非線形含む）、ルンゲクッタ8次（DOP853)などが使えるライブラリ その６：ルンゲクッタ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の説明を参照。

ライブラリ本体：

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

またルンゲクッタ8次のDOP853ルーチンとそのドライバ。

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

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

メルセンヌツイスタ用

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