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

Retarda() Sub Retarda ( N As  Long, F As  LongPtr, T As  Double, Y() As  Double, Tout As  Double, Tend As  Double, RTol() As  Double, ATol() As  Double, Mode As  Long, Grid() As  Double, Cont() As  Double, ICont() As  Long, Info As  Long, Optional Neval As  Long, Optional Nstep As  Long, Optional Naccept As  Long, Optional Nreject As  Long, Optional MaxIter As  Long = 0, Optional Nstiff As  Long = 0, Optional Ngrid As  Long = 0, Optional Mxst As  Long = 0, Optional Cnt As  Long = 0, Optional Hinit As  Double = 0, Optional Hmax As  Double = 0, Optional Fac1 As  Double = 0, Optional Fac2 As  Double = 0, Optional Safe As  Double = 0, Optional Beta As  Double = 0  ) 遅延微分方程式の数値解 (5(4)次 ドルマン・プリンス法) 目的 本プログラムは1階の遅延微分方程式の初期値問題 dy/dt = f(t, y(t), y(t - τ), y(t - τ2), ...), ただし t = t0 において y = y0 の解を求める. ただし, y は要素数 n のベクトルで表され, 方程式は n 本の連立微分方程式である. また, t0 および y0 はそれぞれ t および y の既知の初期値である. f は過去の値 y(t - τ), y(t - τ2), ... にも依存する. 本プログラムは, 遅延微分方程式用の 5(4)次 ドルマン・プリンス法プログラム RETARD (文献 (1)) を書き直したものである. 引数 [in] N 微分方程式の数. (N >= 1) [in] F 微分方程式の関数値を求めるユーザー定義サブルーチンで, 次のように定義すること. Sub F(N As Long, T As Double, Y() As Double, Yp() As Double) dy/dt を求め Yp() に入れる. End Sub ただし, N は方程式の数, dy/dt (= f(t, y(t), y(t - τ), y(t - τ2), ...)) は与えられた T および Y() における微分の計算値である. y(t - τ), y(t - τ2), ... は Ylaga を使用して求めることができる. [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 に達した場合にはその時点で積分を終了し T = Tend, Info = 0 として戻る. [in] Tend 積分を終了する点を表す. Tend に達すると積分を終了し, T = Tend, Info = 0 として戻る. 積分を行うのは前進方向 (Tend > T) でなければならない. [in] Rtol() 配列 Rtol(LRtol - 1) (LRtol >= 1) (Rtol()の全要素 >= 0) 求める解の精度を指定する相対誤差許容値を表す. 本パラメータは LRtol および LAtol の値によりスカラーまたは配列を選択できる. (LRtol < N または LAtol < N の場合: スカラー, LRtol >= N かつ LAtol >= N の場合: 配列) Rtol は Atol と共に各ステップにおける局所誤差テストに用いられ, Y() の各要素が次式を満たすように自動的に選ばれたステップ幅を用いて積分が行われる. スカラーの場合 (i = 0 〜 N - 1): (Y(i) の局所誤差) <= Rtol(0)*Abs(Y(i)) + Atol(0) ただし, Rtol(0) と Atol(0) が同時に 0 であってはならない. 配列の場合 (i = 0 〜 N - 1): (Y(i) の局所誤差) <= Rtol(i)*Abs(Y(i)) + Atol(i) ただし, Rtol(i) と Atol(i) が同時に 0 であってはならない. [in] ATol() 配列 Atol(LAtol - 1) (LAtol >= 1) (Atol()の全要素 >= 0) 求める解の精度を指定する絶対誤差許容値を表す. 本パラメータは LRtol および LAtol の値によりスカラーまたは配列を選択できる. (LRtol < N または LAtol < N の場合: スカラー, LRtol >= N かつ LAtol >= N の場合: 配列) Atol は Rtol と共に各ステップにおける局所誤差テストに用いられる (上記 Rtol 参照). [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] Grid() 配列 Grid(LGrid - 1) (LGrid >= Ngrid) 積分区間内に解の導関数の不連続点(グリッド点)がある場合, それを指定して計算精度および効率を改善することができる. Ngrid > 0 の場合には, グリッド点を, Grid(0) 〜 Grid(Ngrid - 1) に昇順に格納しておかなければならない (積分区間の始点と終点を除く). Ngrid = 0 の場合には Grid() は参照されない. [in,out] Cont() 配列 Cont(LCont - 1) 制御情報領域. 内容を変更してはならない. Ylaga と情報をを共有するためにも使用される. 本プログラムを Info = -1 として呼び出すと, 他の処理は行わず, 必要な LCont の値を計算し Info に返す. [in,out] ICont() 配列 Icont(LICont - 1) (LICont >= 50) 整数制御情報領域. 内容を変更してはならない. Ylaga と情報を共有するためにも使用される. [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: (エラー) ステップ幅が小さくなりすぎたため計算が継続できない. = 13: (エラー) スティフな問題の可能性がある. [out] Neval (省略可) 関数評価回数. [out] Nstep (省略可) 全ステップ数. [out] Naccept (省略可) 採用されたステップ数. [out] Nreject (省略可) 不採用だったステップ数. [in] Maxiter (省略可) 許される最大ステップ数. (デフォルト値 = 100000) [in] Nstiff (省略可) Nstiff回のステップごとにスティフ性のテストを行う. (デフォルト値 = 1000) Nstiff < 0 であればテストは行わない. [in] Ngrid (省略可) Grid() に設定したグリッド点の数. (デフォルト値 = 0) [in] Mxst (省略可) 過去何ステップ分の情報を保持しておくかを指定する. Ylaga はこれを使用して過去の値 y(t - τ) を求める. (デフォルト値 = 1000) [in] Cnt (省略可) カウンタ (Neval, Nstep, Naccept および Nreject) のリセット方法を指定する. (デフォルト値 = 0) = 0: Info = 0 のときにリセットする. = 1: Info = 0 であってもリセットしない. [in] Hinit (省略可) ステップ幅の初期値. (デフォルト値 = プログラムが自動推定する) [in] Hmax (省略可) ステップ幅の最大値. (デフォルト値 = Abs(Tend - T)) (Hmax = 0 であれば省略時の既定値とみなす) [in] Fac1,Fac2 (省略可) ステップ幅選択パラメータ. (デフォルト値: Fac1 = 0.2, Fac2 = 10) Fac1 <= Hnew/Hold <= Fac2 となるようにステップ幅が選ばれる. [in] Safe (省略可) ステップ幅推定時の安全係数. (0.0001 < Safe < 1) (デフォルト値 = 0.9) [in] Beta (省略可) ステップ幅制御の安定化パラメータ (Beta <= 0.2). (デフォルト値 = 0.04, ただし Beta < 0 の場合 0) 文献 (1) E. Hairer, S.P. Norsett and G. Wanner, "Solving Ordinary Differential Equations. Nonstiff Problems. 2nd edition", Springer Series in Computational Mathematics, Springer-Verlag (1993) (邦訳:「常微分方程式の数値解法Ⅰ 基礎編」スプリンガージャパン (2007)) 使用例 次の遅延微分方程式を解く. y'1(t) = -y1(t)*y2(t - 1) + y2(t - 10) y'2(t) = y1(t)*y2(t - 1) - y2(t) y'3(t) = y2(t) - y2(t - 10) (t <= 0 において y1(t) = 5, y2(t) = 0.1, y3(t) = 1) (t = 1, 2, ..., 10, 20, 40 は不連続点) Dim Cont() As Double, ICont() As Long Sub Phi(I As Long, N As Long, T As Double, Y() As Double, Info As Long) If I = 1 Then Y(0) = 5 If I = 2 Then Y(0) = 0.1 If I = 3 Then Y(0) = 1 End Sub Sub F4(N As Long, T As Double, Y() As Double, Yp() As Double) Dim Y2L1(0) As Double, Y2L10(0) As Double Call Ylaga(2, N, T - 1, Y2L1(), AddressOf Phi, Cont(), ICont()) Call Ylaga(2, N, T - 10, Y2L10(), AddressOf Phi, Cont(), ICont()) Yp(0) = -Y(0) * Y2L1(0) + Y2L10(0) Yp(1) = Y(0) * Y2L1(0) - Y(1) Yp(2) = Y(1) - Y2L10(0) End Sub Sub Ex_Retarda() Const N = 3, Ngrid = 11 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout As Double Dim RTol(0) As Double, ATol(0) As Double, Grid(Ngrid) As Double Dim Mode As Long, Neval As Long, Info As Long, I As Long RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) For I = 0 To 9 Grid(I) = I + 1 Next Grid(10) = 20 Mode = 2 Info = -1 Call Retarda(N, AddressOf F4, T, Y(), Tout, Tend, RTol(), ATol(), Mode, Grid(), Cont(), ICont(), Info, Ngrid:=Ngrid) ReDim Cont(Info - 1), ICont(50 - 1) T = 0: Y(0) = 5: Y(1) = 0.1: Y(2) = 1 Tend = 40 Info = 0 Do Tout = T + 5 Call Retarda(N, AddressOf F4, T, Y(), Tout, Tend, RTol(), ATol(), Mode, Grid(), Cont(), ICont(), Info, Neval, Ngrid:=Ngrid) Debug.Print T, Y(0), Y(1), Y(2) Loop While Info >= 1 And Info <= 3 Debug.Print Neval, Info End Sub 注 - 同じプログラムで Mode を 2 の代わりに 3 に変えれば密出力(補間)を使用する. 実行結果 5 0.253384515787142 0.904747157404 4.94186832680886 10 0.332855516159117 3.91808850328078E-02 5.72796359880807 15 4.40303383018292 0.163106427890271 1.5338597419268 20 0.170673965282798 0.864389000478564 5.06493703423863 25 0.214955223950597 1.64093322434754E-02 5.86863544380592 30 4.87247654002283 7.33384930344983E-02 1.15418496694268 35 0.42293063684405 1.35930879446553 4.31776056869042 40 9.12491214440587E-02 2.02995001131839E-02 5.98845137844276 4438 0