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

Retard_r() Sub Retard_r ( N As  Long, T As  Double, Y() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Ngrid As  Long, TGrid() As  Double, RCont() As  Double, ICont() As  Long, Info As  Long, TT As  Double, YY() As  Double, YYp() As  Double, Irtrn As  Long, IRev As  Long, Optional Iout As  Long = 0, Optional Neval As  Long, Optional Nstep As  Long, Optional Naccept As  Long, Optional Nreject As  Long, Optional Nrdens As  Long = 0, Optional Hinit As  Double = 0, Optional Hmax As  Double = 0, Optional MaxIter As  Long = 0, Optional Nstiff As  Long = 0, Optional Safe As  Double = 0, Optional Fac1 As  Double = 0, Optional Fac2 As  Double = 0, Optional Beta As  Double = 0, Optional Cnt As  Long = 0  ) 遅延微分方程式の初期値問題 (5(4)次 ドルマン・プリンス法) (リバースコミュニケーション版) 目的 本ルーチンは1階の遅延微分方程式の初期値問題 dy/dt = f(t, y(t), y(t-τ), ...), ただし t = t0 において y = y0 の解を求める. ただし, t0およびy0は既知でそれぞれtおよびyの初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. Retardは5(4)次ドルマン・プリンス法に基づく陽的ルンゲ・クッタ法のプログラムである. ステップ制御アルゴリズムと密出力機能を持っている. 詳細については下記参考文献を参照せよ. Retard_rはRetardのリバースコミュニケーション版である. 引数 [in] N 微分方程式の数. (N >= 1) [in,out] T 本ルーチンはTからToutまでの積分を行う. 積分を開始する点を与え, 最終ステップの最後の点が返される. [in] 独立変数 t の初期値. [out] ステップ終了時の最終値(通常Toutに等しい). 解がこの点まで正常に求められたことを示す. 新たなToutで Info = 1 として再呼び出しすることにより, 新たな点まで積分を継続することができる. [in,out] Y() 配列 Y(LY - 1) (LY >= N) [in] Tの初期値における従属変数Y()の初期値. [out] 最終のT(通常Toutに等しい)において求められた解(の近似値). [in] Tout 解を求めたい点を設定する. 積分を行うのはtについて前進方向(Tout > T)でも後退方向(Tout < T)でもよい. 要求精度を満たすように自動的に選ばれたステップ幅を用いてTからToutに向かって解が求められる. [in] RTol() 配列 RTol(LRTol - 1) (LRTol >= 1) (RTol()の全要素 >= 0) 求める解の精度を指定する相対誤差許容値. 本パラメータはスカラー(LRTol = 1)またはベクトル(LRTol = N)のどちらでもよい. LRTol = 2, ... または N-1の場合, LRTol = 1 とみなす. LRTol > Nの場合, LRTol = N とみなす. LRTol = N であっても LATol = 1 であれば LRTol = 1 とみなす. 許容値は各ステップにおける局所的な誤差テストに使われ, Y(i)の各要素がおおむね次式を満たすようにする (i = 0 〜 LRTol-1).   abs(Y(i)の局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i) RTol(i) = 0 と設定するとその要素については純粋に絶対誤差テストとなる. RTol(i)とATol(i) (i = 0 〜 LRTol-1) が同時に0になってはいけない. [in] ATol() 配列 ATol(LATol - 1) (LATol >= 1) (ATol()の全要素 >= 0) 求める解の精度を指定する絶対誤差許容値. 本パラメータはスカラー(LATol = 1)またはベクトル(LATol = N)のどちらでもよい. LATol = 2, ... または N-1の場合, LATol = 1 とみなす. LATol > Nの場合, LATol = N とみなす. LATol = N であっても LRTol = 1 であれば LATol = 1 とみなす. 許容値は各ステップにおけるローカルな誤差テストに使われ, Y()の各要素がおおむね次式を満たすようにする (i = 0 〜 LATol-1).   abs(Y(i)の局所誤差) <= RTol(i)*abs(Y(i)) + ATol(i) ATol(i) = 0 と設定するとその要素については純粋に相対誤差テストとなる. RTol(i)とATol(i) (i = 0 〜 LRTol-1) が同時に0になってはいけない. [in] Ngrid Tgrid()で指定するグリッド点数. (Ngrid >= 1) 積分区間中に解またはその導関数の不連続点(グリッド点)がある場合, それを指定しておくと計算精度および効率を改善することができる. Ngridは解を求める最終点を含めて数える. すなわち, 不連続点を全く指定しない場合, Ngrid = 1である. (Ngrid <= 0 であれば Ngrid = 1 とみなす) [in] Tgrid() 配列 Tgrid(LTgrid - 1) (LTgrid >= Ngrid) グリッド点のtの値を Tgrid(0), Tgrid(1), ..., Tgrid(Ngrid - 1) に入力する. 最終要素 Tgrid(Ngrid - 1)には解を求める点 (= Tout) を設定すること. 不連続点を全く指定しない場合(Ngrid = 1の場合), Tgrid(0) = Tout と自動的に設定される. [out] RCont() 配列 RCont(LRCont - 1) (LRCont >= Mxst*(5*Nrdens) 解の後方値のための制御情報. Mxstは記録する(Ylagによりさかのぼることができる)ステップ数である. RTol, ATolが特に小さくなければ Mxst = 500 程度になるようにLRContの大きさを決めるとよい. [out] ICont() 配列 ICont(LICont - 1) (LICont >= Nrdens) 解の後方値のための制御情報. [in,out] Info [in] = 0: 初期化して計算を開始する(新たな問題を解く). = 1: Toutだけを変えて計算を続ける(前回呼び出しの計算を再開する). Ngrid = 1 のときにのみ有効. [out] = -1: パラメータ N の誤り. (N < 1) = -3: パラメータ Y() の誤り. = -5: パラメータ RTol() の誤り. (RTol(i) < 0, RTol(i) = 0 かつ ATol(i) = 0) = -6: パラメータ ATol() の誤り. (ATol(i) < 0) = -8: パラメータ Tgrid() の誤り. = -9: パラメータ RCont() の誤り. = -10: パラメータ ICont() の誤り. = -11: パラメータ Info の誤り. (Info <> 0 かつ Info <> 1) = -16: パラメータ IRev の誤り. (IRev != 0, 1, 5) = 1: 正常終了. = 2: Irtrn による中断 (正常終了). = 11: 最大ステップ数を超えた. = 12: ステップサイズが小さくなり過ぎた. = 13: スティフな問題の可能性がある (中断された). = 14: Ylagによる中断. (RTol, ATolが小さすぎるか, Mxstが小さい可能性がある) [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] Irtrn [in] IRev = 5: 積分を中断したい場合以外変更してはならない. Irtrnを負の値に設定すると, 積分を中断して Info = 2 で終了する. [out] IRev = 5: 初回, 中間, 最後にIRev = 5 で戻るとき, それぞれ 0, 1, 2 を返す. [in,out] IRev リバースコミュニケーションの制御変数 [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の場合, 下記処理を行いIRevを変更せずに再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかはInfoをチェックすること. = 1: YYp()に与えられたTTおよびYY()における微分係数を設定する. すなわち, YYp(i) = dyi/dt = fi(t, y(t), y(t-τ), ... ) (i = 0 〜 N-1) である. YYp()以外の変数を変更してはならない. y(t - τ) は次のように Ylag_r を使って求めることができる. IRev = 0 Do Call Ylag_r(I, T, RCont(), ICont(), YY, IRev) If IRev = 1 Then YY = TにおけるY(I)の値 Loop While IRev <> 0 ただし, I は要素の番号(0 〜 N - 1), また T = t - τ である. RCont() および ICont() は変更せずに渡すこと. IRev = 0 のときに YY に結果が入る. t0 - τ <= T <= t0 (ただし, t0 は t の初期値) であれば IRev = 1 で戻るので, T における Y(I) の値を YY に与える. = 5: 途中結果を出力する. 通常はどの変数も変更してはならない(Ioutを参照せよ). [in] Iout (省略可) 計算の途中結果を出力するかどうか指定する (省略時 = 0) = 0: 出力しない. (IRev = 5では戻らない) = 1: 成功したステップごとに出力する(IRev = 5で戻る). リバースコミュニケーション版でないRetardのSoloutの呼び出しに相当する. 対応する情報は次のとおりである. Nr = Naccept + 1, Told = 前回のT, T = 現在のT, Y() = 現在のY(). RCont() および ICont() は密出力のために Ylag_r を呼び出すときに使用する. (Ioutが上記以外の値であれば省略時の既定値とみなす) [out] Neval (省略可) 関数評価回数. [out] Nstep (省略可) ステップ数. [out] Naccept (省略可) 採用されたステップ数. [out] Nreject (省略可) 不採用だったステップ数. (最初のステップはカウントされない) [in] Nrdens (省略可) 記録する(Ylagによりさかのぼることができる)要素数. (省略時 = N) 要素 0〜Nrdens-1 がYlagのために記録される. (Nrdens <= 0 あるいは Nrdens > N であれば省略時の既定値とみなす) [in] Hinit (省略可) ステップ幅の初期値. (省略時 = 初期関数値より自動推定) (Hinit = 0 であれば省略時の既定値とみなす) [in] Hmax (省略可) ステップ幅の最大値. (省略時 = Tout - T) (Hmax = 0 であれば省略時の既定値とみなす) [in] MaxIter (省略可) 許されるステップ数の最大値. (省略時 = 100000) (MaxIter <= 0 であれば省略時の既定値とみなす) [in] Nstiff (省略可) Nstiff回のステップごとにスティフ性のテストを行う. (省略時 = 1000) Nstiff < 0 であればテストは行わない. (Nstiff = 0 であれば省略時の既定値とみなす) [in] Safe (省略可) ステップ幅推定時の安全係数. (0.0001 < Safe < 1) (省略時 = 0.9) (Safe <= 0.0001 または Safe >= 1 であれば省略時の既定値とみなす) [in] Fac1 (省略可) [in] Fac2 (省略可) ステップ幅選択パラメータ. (省略時: Fac1 = 0.2, Fac2 = 10) Fac1 <= Hnew/Hold <= Fac2 となるようにステップ幅が選ばれる. (Fac1=0 あるいは Fac2=0 であればそれぞれ省略時の既定値とみなす) [in] Beta (省略可) ステップ幅制御の安定化パラメータ. (Beta <= 0.2) (省略時 = 0.04) (Beta < 0 または Beta > 0.2 であれば 0 とみなす) [in] Cnt (省略可) Neval, Nstep, Naccept, Nrejectをリセットするタイミングを指定する. (省略時 = 0) = 0: 本ルーチンが呼び出されるたびにリセットする. <> 0: Info = 0 で呼び出されたときにだけリセットする. 出典 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) 使用例 (1) 次の遅延微分方程式を解く. 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 は不連続点) Sub Ex_Retard_r() Const N = 3, NTout = 12, Ngrid = 0 Const Mxst = 500, LRCont = Mxst * (5 * N + 2), LICont = N Dim T As Double, Y(N - 1) As Double Dim Tout(NTout - 1) As Double, Tgrid(0) As Double Dim RTol(0) As Double, ATol(0) As Double, ITol As Long, Info As Long Dim RCont(LRCont - 1) As Double, ICont(LICont - 1) As Long, TT As Double Dim YY(N - 1) As Double, YYp(N - 1) As Double, Irtrn As Long, IRev As Long Dim Y2L1 As Double, Y2L10 As Double, I As Long RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) T = 0: Y(0) = 5: Y(1) = 0.1: Y(2) = 1 For I = 0 To 9 Tout(I) = I + 1 Next Tout(10) = 20 Tout(11) = 40 Info = 0 For I = 0 To NTout - 1 IRev = 0 Do Call Retard_r(N, T, Y(), Tout(I), RTol(), ATol(), Ngrid, Tgrid(), RCont(), ICont(), Info, TT, YY(), YYp(), Irtrn, IRev) If IRev = 1 Then Y2L1 = YlagR(1, TT - 1, RCont(), ICont()) Y2L10 = YlagR(1, TT - 10, RCont(), ICont()) YYp(0) = -YY(0) * Y2L1 + Y2L10 YYp(1) = YY(0) * Y2L1 - YY(1) YYp(2) = YY(1) - Y2L10 End If Loop While IRev <> 0 If Info <> 1 Then Debug.Print "Error in Retard_r: Info =", Info Exit For End If Debug.Print T, Y(0), Y(1), Y(2) Next End Sub Function YlagR(II As Long, T As Double, RCont() As Double, ICont() As Long) As Double Dim YY As Double, IRev As Long IRev = 0 Do Call Ylag_r(II, T, RCont(), ICont(), YY, IRev) If IRev = 1 Then YY = 0.1 Loop While IRev <> 0 YlagR = YY End Function 実行結果 1 4.61934967214384 0.338647989701885 1.14200233815428 2 3.71355732102849 0.798576426134181 1.58786625283732 3 2.22539076042276 1.33022795319548 2.54438128638175 4 0.832713255497743 1.39799979205327 3.86928695244899 5 0.253384515785908 0.904747157403704 4.94186832681039 6 0.139846290897889 0.456681488763738 5.50347222033837 7 0.148082371938012 0.223097926838407 5.72881970122358 8 0.193976398804964 0.115041118545248 5.79098248264979 9 0.258076652158558 6.43200487911403E-02 5.7776032990503 10 0.332855516159152 3.91808850328527E-02 5.72796359880799 20 0.170673966811781 0.864389003690028 5.06493702949818 40 9.12491209293587E-02 0.020299500264818 5.98845137880582 使用例 (2) 次の遅延微分方程式を解く(グリッドと密出力を使用). 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 は不連続点) Sub Ex_Retard_r_2() Const N = 3, Ngrid = 12 Const Mxst = 500, LRCont = Mxst * (5 * N + 2), LICont = N Dim T As Double, Y(N - 1) As Double Dim Tend As Double, Tout As Double, Tgrid(Ngrid - 1) As Double Dim RTol(0) As Double, ATol(0) As Double, ITol As Long, Info As Long Dim RCont(LRCont - 1) As Double, ICont(LICont - 1) As Long, TT As Double Dim YY(N - 1) As Double, YYp(N - 1) As Double, Irtrn As Long, IRev As Long Dim Iout As Long, Y0 As Double, Y1 As Double, Y2 As Double Dim Y2L1 As Double, Y2L10 As Double, I As Long RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) Iout = 1 T = 0: Y(0) = 5: Y(1) = 0.1: Y(2) = 1 For I = 0 To 9 Tgrid(I) = I + 1 Next Tgrid(10) = 20 Tgrid(11) = 40 Tend = 40 Tout = 5 Info = 0 IRev = 0 Do Call Retard_r(N, T, Y(), Tend, RTol(), ATol(), Ngrid, Tgrid(), RCont(), ICont(), Info, TT, YY(), YYp(), Irtrn, IRev, Iout) If IRev = 1 Then Y2L1 = YlagR(1, TT - 1, RCont(), ICont()) Y2L10 = YlagR(1, TT - 10, RCont(), ICont()) YYp(0) = -YY(0) * Y2L1 + Y2L10 YYp(1) = YY(0) * Y2L1 - YY(1) YYp(2) = YY(1) - Y2L10 ElseIf IRev = 5 Then While T >= Tout Y0 = YlagR(0, Tout, RCont(), ICont()) Y1 = YlagR(1, Tout, RCont(), ICont()) Y2 = YlagR(2, Tout, RCont(), ICont()) Debug.Print Tout, Y0, Y1, Y2 Tout = Tout + 5 Wend End If Loop While IRev <> 0 If Info <> 1 Then Debug.Print "Error in Retard_r: Info =", Info End Sub Function YlagR(II As Long, T As Double, RCont() As Double, ICont() As Long) As Double Dim YY As Double, IRev As Long IRev = 0 Do Call Ylag_r(II, T, RCont(), ICont(), YY, IRev) If IRev = 1 Then YY = 0.1 Loop While IRev <> 0 YlagR = YY End Function 実行結果 5 0.253384515785908 0.904747157403704 4.94186832681039 10 0.332855516159152 3.91808850328527E-02 5.72796359880799 15 4.40303382758445 0.163106426937532 1.53385974547801 20 0.170673966811781 0.864389003690028 5.06493702949818 25 0.214955222878298 1.64093322939858E-02 5.86863544482771 30 4.87247653187963 7.33384924459654E-02 1.1541849756744 35 0.422930643367461 1.3593087991098 4.31776055752274 40 9.12491209293587E-02 0.020299500264818 5.98845137880582