Dormand-Prince(ルンゲ・クッタ8次）DOP853をExcel VBAに移植(説明編)

さて、DOP853を移植したものの説明をします。

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

関数の引数はこうなっている。

Sub DOP853(N As Long, X As Double, _
           Y() As Double, _
           XEND As Double, _
           RTOL() As Double, _
           ATOL() As Double, _
           ITOL As Long, _
           IOUT As Long, _
           WORK() As Double, LWORK As Long, _
           IWORK() As Long, LIWORK As Long, _
           RPAR() As Double, IPAR() As Long, _
           IDID As Long)

'     ----------------
'     N           変数の数
'     X           Xの初期値
'     Y(N)        Yの初期値（N個）
'     XEND        Xの最終値
'     RTOL,ATOL   相対/絶対誤差許容値。変数ごとに設定もできるし、一括も可。
'     ITOL        RTOLとATOLの設定
'                   ITOL=0: RTOL(1)とATOL(1)を全体に適用する。
　　　　　　　　　　　誤差は RTOL*ABS(Y(I))+ATOL以下
'                   ITOL=1: 各変数ごとに誤差設定を変える。配列使用。
'                     誤差はRTOL(I)*ABS(Y(I))+ATOL(I).
'     IOUT        出力のオプション
'                    IOUT=0: 出力ルーチンは呼ばれない。
'                    IOUT=1: 成功したステップの後呼ばれる。
'                    IOUT=2: 成功したステップの後、密出力される。
'                            (IWORK(5)設定が必要)
'                    IOUT=3: ユーザが指定した刻み幅で密出力
'     WORK        "LWORK"の大きさの作業用配列
'                 WORK(1),...,WORK(20) が使われる。
'                 通常は0に設定。
'                 "LWORK" は11*N+8*NRDENS+21の大きさ。ここでNRDENS = IWORK(5)
'     IWORK       "LIWORK"の大きさの作業用配列
'                 IWORK(1),...,IWORK(20)が使われる。
'                 通常は0に設定
'                 "LIWORK"は NRDENS+21 以上
'     RPAR, IPAR  FCNで使われるパラメータ（FCNは後述）

でこの他に、FCNとSOLOUTというサブルーチンが必要。それがこちらに書いてある。

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

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

のように設定する。これは見てもらえばわかると思います。パラメータをRPARとIPARで与える。この例はローレンツ方程式。

具体的な計算例は明日から。

ライセンスは

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

に従います。