Derkf Derkf_r DerkfInt

Derkf() Sub Derkf ( N As  Long, F As  LongPtr, T As  Double, Y() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Info As  Long, Optional Mode As  Long = 0, Optional Cont As  LongPtr  ) 常微分方程式の初期値問題 (5(4)次 ルンゲ・クッタ・フェールベルグ法) 目的 本ルーチンは1階の常微分方程式の初期値問題 dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, t0およびy0は既知でそれぞれtおよびyの初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. 本ルーチンは微分方程式ソルバー・パッケージDEPACに収録されているプログラムのひとつである(DEPACは, Derkf, DeabmおよびDebdfからなる). Derkfは5(4)次ルンゲ・クッタ・フェールベルグ法のプログラムである. 3つの中ではアルゴリズムおよび使用法の両面で最もシンプルなルーチンである. 主に, 微分係数の計算に時間がかからない非スティフおよびややスティフな微分方程式を解くために設計されている. 一般的には, 高精度な解を求める場合や非常に多くの点における解を求める場合には適さない. Derkfはオーバーヘッドが非常に小さいため, 中程度の精度が必要で方程式の計算に時間がかからない場合に, 通常は最も短い計算時間で解を求める. Derkfは H. A. Watts および L. F. Shampine によるプログラム RKF45 の修正版のドライバー・ルーチンである. 下記参考文献を基にした密出力機能が提供される(Modeおよび使用例(2)を参照せよ). 引数 [in] N 微分方程式の数. (N >= 1) [in] F 微分方程式の値を求めるユーザー定義サブルーチンで, 次のように定義すること. Sub F(N As Long, T As Double, Y() As Double, Yp() As Double) Yp(i) = TおよびY()における微分係数の計算値 (i = 0 〜 N-1) End Sub ただし, Nは方程式の数, Yp()は与えられたTおよびY()における微分係数の計算値, すなわち, Yp(i) = dyi/dt = fi(T, Y(0), ..., Y(N-1)) (i = 0 〜 N-1). Yp()以外の変数を変更してはならない. [in,out] T 本ルーチンはTからToutまでの積分を行う. 積分を開始する点を与え, 最終ステップの最後の点が返される. 新たなToutにおける解を求めるために積分を継続することができる(インターバル・モードとよぶ). その場合, 2回目以降の呼び出しでは, Tは前回のToutに等しくなければならない. また, Toutまでの各中間ステップにおいて戻ることもできる(中間結果出力モードとよぶ). 本モードは解の挙動をみたい場合に使うとよい. モードはパラメータModeで指定できる. [in] 独立変数Tの初期値. [out] 独立変数Tの最終ステップの最後の点の値 (インターバル・モードでは通常Toutに等しい). 解がこの点まで正常に求められたことを示す. [in,out] Y() 配列 Y(LY - 1) (LY >= N) [in] Tの初期値における従属変数Y()の初期値. [out] 最終のT(インターバル・モードにおいてはToutに等しい)において求められた解(の近似値). [in] Tout 解を求めたい点を設定する. Tout = T とすることができ, その場合はTにおける微分係数を計算して戻る. 積分を行うのはtについて前進方向(Tout > T)でも後退方向(Tout < T)でもよい. ただし, 初期化時でなければ積分方向を変更することはできない. 要求精度を満たすように自動的に選ばれたステップ幅を用いてTからToutに向かって解が求められる. 必要であれば, 各中間ステップにおける解とその微分係数を確認するために戻ることができる(中間結果出力モード). ただし, その場合であっても本来のとおりToutを指定しなければならない. [in] RTol RTol(LRTol - 1) (LRTol >= 1) (RTol()の全要素 >= 0) 求める解の精度を指定する相対誤差許容値. 本パラメータはスカラー(LRTol = 1)または配列(LRTol = N)のどちらでもよい. 配列は誤差テストをより細かく制御したい場合に使うことができる. LRTol = 2, ... または N-1の場合, LRTol = 1 とみなす. LRTol > Nの場合, LRTol = N とみなす. 許容値は各ステップにおける局所的な誤差テストに使われ, Y()の各要素がおおむね次式を満たすようにする (i = 0 〜 N-1).   abs(局所誤差) <= RTol(0)*abs(Y(i)) + ATol(0) (LRTolまたはLATol = 1 の場合)     または   abs(局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i) (LRTolおよびLATol = N の場合) RTol(i) = 0 と設定するとその要素については純粋に絶対誤差テストとなる. 誤差が約1.0e-8以下の相対精度を必要とするならば, 通常はDerkfを使わないほうがよい. DEPACのDeabmルーチンは厳しい精度を必要とする場合により効率的である. [in] ATol ATol(LATol - 1) (LATol >= 1) (ATol()の全要素 >= 0) 求める解の精度を指定する絶対誤差許容値. 本パラメータはスカラー(LATol = 1)または配列(LATol = N)のどちらでもよい. 配列は誤差テストをより細かく制御したい場合に使うことができる. LATol = 2, ... または N-1の場合, LATol = 1 とみなす. LATol > Nの場合, LATol = N とみなす. 許容値は各ステップにおける局所的な誤差テストに使われ, Y()の各要素がおおむね次式を満たすようにする (i = 0 〜 N-1).   abs(局所誤差) <= RTol(0)*abs(Y(i)) + ATol(0) (LRTolまたはLATol = 1 の場合)     または   abs(局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i) (LRTolおよびLATol = N の場合) ATol(i) = 0 と設定するとその要素については純粋に相対誤差テストとなる. 誤差が約1.0e-8以下の相対精度を必要とするならば, 通常はDerkfを使わないほうがよい. DEPACのDeabmルーチンは厳しい精度を必要とする場合により効率的である. [in,out] Info [in] 制御コード. = 0: 問題の最初の呼び出し時(新たに問題を開始する場合)には Info = 0 と設定せよ. 初期化が行われ, 新たな問題についての計算が開始される. = 1, 2, 11〜15: Info = 1, 2, 11〜15 で戻った場合, 計算を継続するために, Infoの値を変えずに再呼び出しすることができる(下記説明を参照せよ). [out] リターンコード. ユーザーはこの値をチェックし, Info = 1, 2, 11〜15 であれば必要により次のアクションとして本ルーチンを再度呼び出すことができる. = -1: パラメータ N の誤り. (N < 1) = -3: パラメータ T の誤り. (継続呼び出しで T = Tout または T <> 前回のTout) = -4: パラメータ Y() の誤り. = -5: パラメータ Tout の誤り. (積分の方向を変更した) = -6: パラメータ Rtol() の誤り. (RTol(i) < 0) = -7: パラメータ Atol() の誤り. (ATol(i) < 0) = -8: パラメータ Info の誤り. (Info <> 0, 1, 2, 11〜15) = 1: 正常終了 (T = Tout). Toutを再設定して再呼び出しすることができる. ToutはTと異なる値であること. = 2: 中間結果出力モードのため戻った. (まだ Toutに達していない). Tout方向の次のステップを続行するために再呼び出しすることができる. = 11: 最大ステップ数(10000)を超えた. 再呼び出しして続行することができる. さらに10000ステップ続けることができる. = 12: 誤差許容値が厳しすぎる. 許容値を変更してから再呼び出しすることができる. = 13: 求められた解が0であるため, 純粋な相対誤差テスト(ATol = 0)ができない. Atolを正の値にして再呼び出しして続行することができる. = 14: スティフな問題の可能性がある. 再呼び出しして続行することができる. しかし, スティフな問題を効率的に扱うことができるDEPACのDebdfのような他のプログラムが推奨される. = 15: 頻繁な出力のためステップサイズが制限された. 再呼び出しして続行することができる. しかし, 密出力機能を持った他のプログラムが推奨される. = 16: 無限ループが検出された. [in] Mode (省略可) 動作モード. (省略時 = 0) = 0: Toutで戻る (インターバル・モード). = 1: ステップごとに戻る (中間結果出力モード) = 2: ステップごとに戻る (中間結果出力モード). 密出力を可能にする. すなわち, Info = 2 (最後のステップでは Info = 1) で戻ったときに DerkfInt ルーチンを使って直近のステップの区間内で補間値を求めることができる. (他の値であれば省略時の既定値とみなす.) [out] Cont (省略可) DerkfInt が使用する密出力のための情報. Info = 2 (最後のステップでは Info = 1) かつ Mode = 2 の場合にのみ返される. 出典 SLATEC (DEPAC) 参考文献 W H Enright et al. "Interpolants for Runge-Kutta Formulas" ACM Transactions on Mathematical Software Vol.12, No.3, 1986, pp.193-218 使用例 (1) 次の常微分方程式の初期値問題を解く. dy1/dt = -2*y1 + y2 - cos(t) dy2/dt = 2*y1 - 3*y2 + 3*cos(t) - sin(t) (t = 0 において y1 = 1, y2 = 2) Sub F1(N As Long, T As Double, Y() As Double, Yp() As Double) Yp(0) = -2 * Y(0) + Y(1) - Cos(T) Yp(1) = 2 * Y(0) - 3 * Y(1) + 3 * Cos(T) - Sin(T) End Sub Sub Ex_Derkf() Const N = 2 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double Dim RTol(0) As Double, ATol(0) As Double, Info As Long RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) T = 0: Tend = 10: Y(0) = 1: Y(1) = 2 Info = 0 Do Tout = T + 1 Call Derkf(N, AddressOf F1, T, Y(), Tout, RTol(), ATol(), Info) If Info <> 1 Then Debug.Print "Error in Derkf: Info =", Info Exit Do End If Debug.Print T, Y(0), Y(1) Loop While Tout < Tend End Sub 実行結果 1 0.367879441136182 0.908181747107717 2 0.135335283245885 -0.280811553333642 3 4.97870683958675E-02 -0.940205428291068 4 1.83156389160862E-02 -0.635327982031259 5 6.73794699054023E-03 0.290400132477777 6 2.47875214648626E-03 0.9626490388872 7 9.11881936224063E-04 0.754814136368097 8 3.35462623309683E-04 -0.145164571170354 9 1.23409831706352E-04 -0.91100685213574 10 4.5399959961099E-05 -0.839026129207695 使用例 (2) 次の常微分方程式の初期値問題を解く(密出力を使用). dy1/dt = -2*y1 + y2 - cos(t) dy2/dt = 2*y1 - 3*y2 + 3*cos(t) - sin(t) (t = 0 において y1 = 1, y2 = 2) Sub F1(N As Long, T As Double, Y() As Double, Yp() As Double) Yp(0) = -2 * Y(0) + Y(1) - Cos(T) Yp(1) = 2 * Y(0) - 3 * Y(1) + 3 * Cos(T) - Sin(T) End Sub Sub Ex_Derkf_2() Const N = 2 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double Dim RTol(0) As Double, ATol(0) As Double, Info As Long Dim Y1(N - 1) As Double, Mode As Long, Cont As LongPtr RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) T = 0: Tend = 10: Y(0) = 1: Y(1) = 2 Mode = 2 Tout = T + 1 Info = 0 Do Call Derkf(N, AddressOf F1, T, Y(), Tend, RTol(), ATol(), Info, Mode, Cont) If Info <> 1 And Info <> 2 Then Debug.Print "Error in Derkf: Info =", Info Exit Do End If While T >= Tout Call DerkfInt(N, Tout, Y1(), Cont) Debug.Print Tout, Y1(0), Y1(1) Tout = Tout + 1 Wend Loop While Tout <= Tend End Sub 実行結果 1 0.367879441127794 0.908181747124133 2 0.135335283265252 -0.280811553374618 3 4.97870684046614E-02 -0.940205428309061 4 1.83156389289378E-02 -0.635327982057511 5 6.73794698644643E-03 0.290400132486541 6 2.4787521426282E-03 0.962649038895081 7 9.11881915710959E-04 0.754814136410194 8 3.35462644956918E-04 -0.145164571216893 9 1.23409854714012E-04 -0.911006852183254 10 4.53999601058471E-05 -0.839026129207973