Deabm Deabm_r Derkfa Derkfa_r Dop853a Dop853a_r Dopn1210 Dopn1210_r Dopn43 Dopn43_r Dopn64 Dopn64_r Dopn86 Dopn86_r Dopri5a Dopri5a_r Dverka Dverka_r Odex2a Odex2a_r Odexa Odexa_r Retarda Retarda_r Ylaga

Dverka_r() Sub Dverka_r ( N As  Long, T As  Double, Y() As  Double, Tout As  Double, Tend As  Double, Tol As  Double, Mode As  Long, Info As  Long, TT As  Double, YY() As  Double, YYp() As  Double, IRev As  Long, Optional Neval As  Long, Optional Nstep As  Long, Optional Naccept As  Long, Optional Nreject As  Long, Optional MaxIter As  Long = 0, Optional ErrCntl As  Long = 0, Optional Cnt As  Long = 0, Optional Hinit As  Double = 0, Optional Hmax As  Double = 0, Optional Hmin As  Double = 0, Optional Scal As  Double = 0, Optional Efloor As  Double = 0  ) 常微分方程式の初期値問題 (6(5)次 ルンゲ・クッタ・ヴァーナー法) (リバースコミュニケーション版) 目的 本プログラムは1階の常微分方程式の初期値問題 dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, y は要素数 n のベクトルで表され, 方程式は n 本の連立微分方程式である. また, t0 および y0 はそれぞれ t および y の既知の初期値である. 本プログラムは, 6(5)次 ルンゲ・クッタ・ヴァーナー法プログラム DVERK (文献 (1)) をベースに Enrightら (文献 (2)) の補間式を使用した密出力機能を追加して書き直したものである. 本プログラムは Dverka のリバースコミュニケーション版である. 引数 [in] N 微分方程式の数. (N >= 1) [in,out] T 独立変数 t を表す. 本プログラムは t の初期値から Tend までの積分を行う. Mode の設定により, 途中結果を返すために Tout またはステップごとに戻ることができる. [in] t の初期値 t0. [out] 積分終了時には T = Tend, 途中結果を返すために戻ったときには Mode の設定により T = Tout または T = 直前のステップの終点の値 を返す. [in,out] Y() 配列 Y(LY - 1) (LY >= N) 従属変数 y を表す. [in] t の初期値 t0 における y の初期値 y0. [out] t = T における y の値 (数値解). [in] Tout Mode = 2 または 3 の場合に, 途中結果の確認／出力を行う t を表す. Mode = 0 または 1 では Tout は参照されない. Mode = 2 または 3 では t = Tout における y を求め, それぞれ Info = 2 または 3 として戻る. 続いて新たな Tout における解を求めるために積分を継続したい場合には, 新たな Tout に変更して(Info を含む)他の変数を変更せずに再度呼び出しを行うことができる. T < Tout <= Tend でなければならない (ただし, 後退方向に積分を行う場合には Tend <= Tout < T でなければならない). 先に Tend に達した場合にはその時点で積分を終了し T = Tend, Info = 0 として戻る. [in] Tend 積分を終了する点を表す. Tend に達すると積分を終了し, T = Tend, Info = 0 として戻る. 積分を行うのは前進方向 (Tend > T) でも後退方向 (Tend < T) でもよい. [in] Tol 誤差許容値を表す. (Tol > 0) 大域誤差が Tol に比例するように局所誤差のノルムが制御される. ノルムは, ErrCntl パラメータで指定される誤差制御方式による重み付き最大ノルムが使用される. デフォルトの重みは 1/max(1, abs(y(i))) で, 絶対誤差制御および相対誤差制御を混合したものである. [in] Mode 動作モード. 本プログラムは t0 から Tend までの積分を行うが, 中間結果の確認／出力の方法により4つの動作モードが提供される. どのモードでも T = Tend に達したときには積分を終了し, Info = 0 を返す. = 0: Tend まで戻らない. Tout は無視される. = 1: 成功したステップごとに戻る (Info = 1 を返す). Tout は無視される. T にはそのステップの終了点を返す. 最終ステップでは T = Tend となるようにステップ幅が調整される. = 2: 積分途中で Tout における中間結果を返すために戻る (T = Tout, Info = 2 を返す). 続けて次の Tout における解を求めるために積分を継続するためは, Tout を新たな値に変えて再度呼び出しを行うことができる. Tout 直前のステップでは T = Tout となるようにステップ幅が調節される. = 3: 積分途中で Tout における中間結果を返すために戻る (T = Tout, Info = 3 を返す). 続けて次の Tout における解を求めるために積分を継続するためは, Tout を新たな値に変えて再度呼び出しを行うことができる. Mode = 2 と異なり Tout における Y() の値は補間により求められ, Tout 直前のステップでのステップ幅の調節は行わず Mode = 1 と同じステップを踏む. [in,out] Info [in] 制御コード. = 0: 最初の呼び出し時(新たに問題を開始する場合)には Info = 0 と設定する. 全ての変数の初期化を行ってから計算を開始する. = 1, 2, 3: Info = 1, 2 または 3 で戻った場合, 計算を継続するために Tout だけを変更し, Info の値をそのままにして再呼び出しすることができる. [out] リターンコード. = 0: 正常終了. Tend までの積分が完了した. < 0: (-Info)番目の入力パラメータの誤り. = 1: Mode = 1 の中間結果出力のために戻った. 再呼び出しすることより次のステップまで進むことができる. = 2: Mode = 2 の中間結果出力のために T = Tout において戻った. Tout を再設定して再呼び出しすることができる. = 3: Mode = 3 の中間結果出力のために T = Tout において戻った. Tout を再設定して再呼び出しすることができる. = 11: (エラー) 最大ステップ数を超えた. = 12: (エラー) ステップ幅の選択に失敗した (Tol が小さすぎる可能性がある). = 13: (エラー) ステップ幅が小さくなりすぎた (Tol が小さすぎる可能性がある). [out] TT IRev = 1: YYp() 参照. [out] YY() 配列 YY(LYY - 1) (LYY >= N) IRev = 1: YYp() 参照. [in] YYp() 配列 YYp(LYYp - 1) (LYYp >= N) IRev = 1: T = TT および Y = YY() における微分値 dy/dt = f(t, y) を計算し, 再呼び出し時に YYp() に設定する. [in,out] IRev リバースコミュニケーションの制御変数. [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の場合, 下記処理を行い IRev を変更せずに再び本プログラムを呼び出すこと. = 0: 処理終了. 正常終了かどうかは Info をチェックすること. = 1: T = TT および Y = YY() における微分の計算値を YYp() に設定する. YYp() 以外の変数を変更してはならない. [out] Neval (省略可) 関数評価回数. [out] Nstep (省略可) 全ステップ数. [out] Naccept (省略可) 採用されたステップ数. [out] Nreject (省略可) 不採用だったステップ数. [in] Maxiter (省略可) 許される最大ステップ数. (デフォルト値 = 100000) [in] ErrCntl (省略可) エラー制御方法. 誤差制御のための推定誤差ベクトルの重みを指定する. (デフォルト値 = 0) = 0: 重み = 1/max(1, abs(y(i))) (絶対および相対誤差制御を混合したもの) = 1: 重み = 1 (絶対誤差制御) = 2: 重み = 1/abs(y(i)) (相対誤差制御) = 3: 重み = 1/max(abs(Efloor), abs(y(i))) (abs(y(i)) が下限値 abs(Efloor) より小さくない限り相対誤差制御) [in] Cnt (省略可) カウンタ (Neval, Nstep, Naccept および Nreject) のリセット方法を指定する. (デフォルト値 = 0) = 0: Info = 0 のときにリセットする. = 1: Info = 0 であってもリセットしない. [in] Hinit (省略可) ステップ幅の初期値. (デフォルト値 = プログラムが自動推定する) [in] Hmax (省略可) 最大ステップ幅 = min(abs(Hmax), 2/abs(Scal)) とする. (ただし, Scal = 0 のときは 最大ステップ幅 = abs(Hmax)) (デフォルト値: 最大ステップ幅 = 2/abs(Scal) (ただし, Scal = 0 のときは 最大ステップ幅 = 2)) [in] Hmin (省略可) 最小ステップ幅. (デフォルト値 = 10*max(dwarf, rreb*max((y の重み付きノルム)/Tol, abs(T))). ただし, dwarf は非常に小さな正数(マシン精度) (= 1.0e-50), rreb は相対丸め誤差限界値である.) [in] Scal (省略可) 問題のスケール尺度. (デフォルト値 = 1) 本パラメータは, 最大ステップ幅の決定および判定基準の変更のために使われる. この値が大きければより高精度の結果を求めようとする. [in] Efloor (省略可) ErrCntl = 3 の場合の下限値. 文献 (1) netlib/ode (2) W. H. Enright, et al.: "Interpolants for Runge-Kutta Formulas", ACM Transactions on Mathematocal Software, Vol.12, No.3. (1986) 使用例 次の常微分方程式の初期値問題を解く. 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_Dverka_r() Const N = 2 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double Dim Tol As Double, Mode As Long, Neval As Long, 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 Mode = 2 T = 0: Y(0) = 1: Y(1) = 2 Tend = 10 Info = 0 Do Tout = T + 1 IRev = 0 Do Call Dverka_r(N, T, Y(), Tout, Tend, Tol, Mode, Info, TT, YY(), YYp(), IRev, Neval) If IRev = 1 Then Call F1(N, TT, YY(), YYp()) Loop While IRev <> 0 Debug.Print T, Y(0), Y(1) Loop While Info >= 1 And Info <= 3 Debug.Print Neval, Info End Sub 注 - 同じプログラムで Mode を 2 の代わりに 3 に変えれば密出力(補間)を使用する. 実行結果 1 0.367879441171258 0.908181747039973 2 0.135335283236816 -0.28081155331094 3 4.97870683680297E-02 -0.940205428232924 4 1.83156388888794E-02 -0.635327981975176 5 6.73794699886029E-03 0.290400132462794 6 2.47875217650127E-03 0.962649038827388 7 9.11881965405588E-04 0.754814136309173 8 3.35462627940427E-04 -0.145164571180805 9 1.23409804260466E-04 -0.911006852080962 10 4.53999299169156E-05 -0.839026129147014 2539 0