Contd5 Contd5_r Contd8 Contd8_r Contx1 Contx1_r Contx2 Contx2_r Deabm Deabm_r Derkf Derkf_r DerkfInt Dop853 Dop853_r Dopri5 Dopri5_r Doprin Doprin_r Dverk Dverk_r DverkInt Odex Odex2 Odex2_r Odex_r Retard Retard_r Ylag Ylag_r

Dverk_r() Sub Dverk_r ( N As  Long, T As  Double, Y() As  Double, Tout As  Double, Tol As  Double, Info As  Long, TT As  Double, YY() As  Double, YYp() As  Double, IRev As  Long, Optional Neval As  Long, Optional Naccept As  Long, Optional MaxEval As  Long = 0, Optional Int1 As  Long = 0, Optional Int2 As  Long = 0, Optional Hmin As  Double = 0, Optional Hmax As  Double = 0, Optional Scal As  Double = 0, Optional Hstart As  Double = 0, Optional Weight As  Long = 0, Optional Floor As  Double, Optional Cont As  LongPtr  ) 常微分方程式の初期値問題 (6(5)次 ルンゲ・クッタ・ヴァーナー法) (リバースコミュニケーション版) 目的 本ルーチンは1階の常微分方程式の初期値問題 dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, t0およびy0は既知でそれぞれtおよびyの初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. 本ルーチンは6(5)次ルンゲ・クッタ・ヴァーナー法のプログラムである. 微分係数の計算に時間がかからない非スティフな問題に適している. 下記参考文献を基にした密出力機能が提供される(Int2および使用例(2)を参照せよ). Dverk_rはDverkのリバースコミュニケーション版である. 引数 [in] N 微分方程式の数. (N >= 1) [in,out] T 本ルーチンはTからToutまでの積分を行う. 積分を開始する点を与え, 最終ステップの最後の点が返される. [in] 独立変数Tの初期値. [out] 独立変数Tの最終ステップの最後の点の値 (通常Toutに等しい). 解がこの点まで正常に求められたことを示す. Toutの値を変えて本ルーチンを再呼び出しすることにより新たな点まで積分を継続することができる. [in,out] Y() 配列 Y(LY - 1) (LY >= N) [in] Tの初期値における従属変数Y()の初期値. [out] 最終のT(通常Toutに等しい)において求められた解(の近似値). [in] Tout 解を求めたい点を設定する. 積分を行うのはtについて前進方向(Tout > T)でも後退方向(Tout < T)でもよい. 要求精度を満たすように自動的に選ばれたステップ幅を用いてTからToutに向かって解が求められる. [in] Tol 誤差許容値. (Tol > 0) 本ルーチンは, 大域誤差がTolに比例するように, 局所誤差のノルムを制御する. ノルムは重み付きの最大ノルムである. デフォルトの重みは, k番目の要素について, 1/max(1, abs(y(k))) である. これは絶対および相対誤差制御を混合したものになる. [in,out] Info [in] 制御コード. 最初の呼び出し時には Info = 0 としなければならない. = 0: 初期化を行う. = 1: 正常終了後に, Toutを新しい値に再設定して計算を継続する. = 2, 3, 4: 割り込み1または2からの再開を行う. [out] リターンコード. ユーザーはこの値をチェックし, info = 1〜4 であれば必要により次のアクションとして本ルーチンを再度呼び出して計算を続行することができる. = -1: パラメータ N の誤り. (N < 1) = -3: パラメータ T の誤り. (継続呼び出しで T = Tout または T <> 前回のTout) = -6: パラメータ Tol の誤り. (Tol <= 0) = -9: パラメータ lc の誤り. (lcが小さすぎる) = -11: パラメータ lwork の誤り. (lworkが小さすぎる) = -12: パラメータ Info の誤り. (Info < 0 または Info > 4) = 1: 正常終了 (T = Tout). Toutを再設定して再呼び出しすることができる. 他の変数は変更してはならない. = 2: 割り込み 1 (Int1参照). = 3: 割り込み 2 (Int2参照, 試行ステップが採用された). = 4: 割り込み 2 (Int2参照, 試行ステップが採用されなかった). = 11: 関数評価回数が最大値を超えた. = 12: ステップ幅の最小値 > 最大値 になった. (Tolが小さすぎる可能性がある) = 13: ステップ幅の最小値に等しいあるいは小さいステップ幅でも収束しなかった. (Tolが小さすぎる可能性がある) [out] TT IRev = 1: 再呼び出し時に微分係数値を求めてYYp()に入れるべき点のTの値を返す. [out] YY() 配列 YY(LYY - 1) (LYY >= N) IRev = 1: 再呼び出し時に微分係数値を求めてYYp()に入れるべき点のYの値を返す. [in] YYp() 配列 YYp(LYYp - 1) (LYYp >= N) IRev = 1: T(= TT)およびY(= YY())における微分係数, すなわち, YYp(i) = dyi/dt = fi(TT, YY(0), ..., YY(N-1)) (i = 0 〜 N-1) を計算して, 再呼び出し時に入力する. [in,out] IRev リバースコミュニケーションの制御変数 [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の場合, 下記処理を行いIRevを変更せずに再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかはInfoをチェックすること. = 1: TTおよびYY()における微分係数の計算値をYYp()に設定する. YYp()以外の変数を変更してはならない. [out] Neval (省略可) 関数評価回数. [out] Naccept (省略可) 採用されたステップ数. [in] MaxEval (省略可) 関数評価回数の最大値. (省略時 = 5000) (MaxEval < 0 であれば絶対値を使用する) (MaxEval = 0 であれば省略時の既定値とみなす) [in] Int1 (省略可) 割り込みモード 1 の指定. (省略時 = 0) = 0: Toutに達するまで戻らない. (通常モード) = 1: 試行ステップを行う前に計算を中断し戻る. (割り込みモード 1) (上記以外の値であれば省略時の既定値とみなす) [in] Int2 (省略可) 割り込みモード 2 の指定. (省略時 = 0) = 0: Toutに達するまで戻らない. (通常モード) = 1: 試行ステップを行った後に計算を中断し戻る. ステップを採用する場合には Info = 3, 採用しない場合には Info = 4 で戻る. (割り込みモード 2) = 2: 密出力が可能な割り込みモード 2. すなわち, Info = 3 (最後のステップでは Info = 1) のときには DverkInt ルーチンを使って直近のステップの区間内で補間値を求めることができる. 採用されたステップあたり1回の追加関数評価が必要になる. (上記以外の値であれば省略時の既定値とみなす) [in] Hmin (省略可) ステップ幅の最小値. (省略時 = 10*max(1.0e-50, eps*max((yの重み付きノルム)/tol, abs(x))), ただしepsはマシンイプシロン) (Hmin = 0 であれば省略時の既定値とみなす) [in] Hmax (省略可) ステップ幅の最大値を次のようにHmaxとScalから決定する. (省略時 = 0) Hmax = 0 かつ Scal = 0: 2 Hmax <> 0 かつ Scal <> 0: min(abs(Hmax), 2/abs(Scal)) Hmax <> 0 かつ Scal = 0: abs(Hmax) Hmax = 0 かつ Scal <> 0: 2/abs(Scal) [in] Scal (省略可) 信頼度パラメータ. (問題のスケーリング因子) (省略時 = 1) スケーリングおよびステップ幅の最大値の決定に使われる. Scalの値を大きくすると結果の信頼性が上がるが関数評価回数は増える. (Scal = 0 であればスケーリングには省略時の既定値が使われる) [in] Hstart (省略可) ステップ幅の初期値. (省略時 = (ステップ幅の最大値)*Tol^(1/6)) (Hstart < 0 であれば絶対値を使用する) (Hstart = 0 であれば省略時の既定値とみなす) [in] Weight (省略可) 誤差制御には重み付き推定誤差ベクトルの最大ノルムを使うが, その重みを指定する. (省略時 = 0) = 0: 重み = 1/max(1, abs(Y(k))) (下限 = 1 の相対誤差制御) = 1: 重み = 1 (絶対誤差制御) = 2: 重み = 1/abs(Y(k)) (相対誤差制御) = 3: 重み = 1/max(abs(Floor), abs(Y(k))) (下限値を指定した相対誤差制御) (上記以外の値であれば省略時の既定値とみなす) [in] Floor (Weight = 3 以外のとき省略可) Weight = 3 の場合の下限値. [out] Cont (省略可) DverkInt が使用する密出力のための情報. Info = 3 (最後のステップでは Info = 1) かつ Int2 = 2 の場合にのみ返される. 出典 netlib/ode 参考文献 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_Dverk_r() Const N = 2 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double Dim Tol As Double, Info As Long Dim TT As Double, YY(N - 1) As Double, YYp(N - 1) As Double, IRev As Long Tol = 0.0000000001 '1.0e-10 T = 0: Tend = 10: Y(0) = 1: Y(1) = 2 Info = 0 Do Tout = T + 1 IRev = 0 Do Call Dverk_r(N, T, Y(), Tout, Tol, Info, TT, YY(), YYp(), IRev) If IRev = 1 Then Call F1(N, TT, YY(), YYp()) Loop While IRev <> 0 If Info <> 1 Then Debug.Print "Error in Dverk_r: Info =", Info Exit Sub End If Debug.Print T, Y(0), Y(1) Loop While Tout < Tend End Sub 実行結果 1 0.367879441171258 0.908181747039973 2 0.135335283236816 -0.28081155331094 3 4.97870683680296E-02 -0.940205428232924 4 1.83156388888794E-02 -0.635327981975177 5 6.73794699886025E-03 0.290400132462795 6 2.47875217650124E-03 0.962649038827387 7 9.11881965405485E-04 0.754814136309172 8 3.35462627940704E-04 -0.145164571180804 9 1.23409804260653E-04 -0.911006852080961 10 4.53999299169892E-05 -0.839026129147014 使用例 (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_Dverk_r_2() Const N = 2 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double Dim Tol As Double, Info As Long Dim Y1(N - 1) As Double, Int2 As Long, Cont As LongPtr Dim TT As Double, YY(N - 1) As Double, YYp(N - 1) As Double, IRev As Long Tol = 0.0000000001 '1.0e-10 T = 0: Tend = 10: Y(0) = 1: Y(1) = 2 Int2 = 2 Tout = T + 1 Info = 0 Do IRev = 0 Do Call Dverk_r(N, T, Y(), Tend, Tol, Info, TT, YY(), YYp(), IRev, Int2:=Int2, Cont:=Cont) If IRev = 1 Then Call F1(N, TT, YY(), YYp()) Loop While IRev <> 0 While (Info = 1 Or Info = 3) And T >= Tout Call DverkInt(N, Tout, Y1(), Cont) Debug.Print Tout, Y1(0), Y1(1) Tout = Tout + 1 Wend Loop While Tout <= Tend End Sub 実行結果 1 0.367879441761895 0.9081817458561 2 0.135335283109193 -0.280811553054101 3 4.97870683242647E-02 -0.940205428144945 4 0.018315638476038 -0.635327981148217 5 6.73794741978772E-03 0.290400131618177 6 2.4787522740674E-03 0.962649038631452 7 9.11882292284201E-04 0.754814135654083 8 3.35462206406771E-04 -0.14516457033392 9 1.23409537546899E-04 -0.911006851545974 10 4.53999299134007E-05 -0.839026129147008