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() Sub Retard ( N As  Long, F As  LongPtr, T As  Double, Y() As  Double, Tout As  Double, RTol() As  Double, ATol() As  Double, Ngrid As  Long, TGrid() As  Double, Info As  Long, Optional Solout As  LongPtr = NullPtr, Optional Neval As  Long, Optional Nstep As  Long, Optional Naccept As  Long, Optional Nreject As  Long, Optional Mxst As  Long = 0, 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)次ドルマン・プリンス法に基づく陽的ルンゲ・クッタ法のプログラムである. ステップ制御アルゴリズムと密出力機能を持っている. 詳細については下記参考文献を参照せよ. 引数 [in] N 微分方程式の数. (N >= 1) [in] F 微分方程式の値を求めるユーザー定義サブルーチンで, 次のように定義すること. Sub F(N As Long, T As Double, Y() As Double, Yp() As Double, RCont As Double, ICont As Long) Yp(i) = TおよびY()における微分係数の計算値 (i = 0 〜 N-1) End Sub ただし, Nは方程式の数である. Yp()には与えられたTおよびY()における微分係数を設定する. すなわち, Yp(i) = dyi/dt = fi(t, y(t), y(t-τ), ... ) (i = 0 〜 N-1) である. Yp()以外の変数を変更してはならない. y(t-τ)は次の関数呼び出しにより求めることができる. Ylag(i, t-τ, AddressOf Phi, RCont, ICont) [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 と自動的に設定される. [in,out] Info [in] = 0: 初期化して計算を開始する(新たな問題を解く). = 1: Toutだけを変えて計算を続ける(前回呼び出しの計算を再開する). Ngrid = 1 のときにのみ有効. [out] = -1: パラメータ N の誤り. (N < 1) = -4: パラメータ Y() の誤り. = -6: パラメータ RTol() の誤り. (RTol(i) < 0, RTol(i) = 0 かつ ATol(i) = 0) = -7: パラメータ ATol() の誤り. (ATol(i) < 0) = -9: パラメータ Tgrid() の誤り. = -10: パラメータ Info の誤り. (Info <> 0 かつ Info <> 1) = 1: 正常終了. = 2: Soloutによる中断 (正常終了). = 11: 最大ステップ数を超えた. = 12: ステップサイズが小さくなり過ぎた. = 13: スティフな問題の可能性がある (中断された). = 14: Ylagによる中断. (RTol, ATolが小さすぎるか, Mxstが小さい可能性がある) [in] Solout (省略可) 中間結果を出力するユーザー定義サブルーチンで, ステップが成功するごとに呼び出される. 次のように定義すること. (省略時 = NullPtr) Sub Solout(Nr As Long, Told As Double, T As Double, Y() As Double, N As Long, RCont As Double, ICont As Long, Irtrn As Long) Nr番目のステップTでのY()の値が渡されるのでこれを出力する. Toldは前回のTの値, Nは方程式の次数である. Irtrnの値は, 初回, 中間, または, 最終呼び出しのとき, それぞれ 0, 1 または 2 である. Irtrnを使用して積分を中断することができる. SoloutでIrtrnを負の値に設定して戻ると積分を中断して Info = 2 で終了する. Soloutの中で解を変更した場合には Irtrn = 3 と設定して戻ること. 制御情報 RCont および ICont を使用して密出力を行うことができる. 区間[Told, T]内の任意の点T2での解Y(i) (0 <= i <= N-1)を次の関数呼び出しにより得ることができる. Y(i) = Ylag(i, T2, AddressOf Phi, RCont, ICont) End Sub これを省略した場合(Solout = NullPtrの場合), 中間結果出力を行わない. [out] Neval (省略可) 関数評価回数. [out] Nstep (省略可) ステップ数. [out] Naccept (省略可) 採用されたステップ数. [out] Nreject (省略可) 不採用だったステップ数. (最初のステップはカウントされない) [in] Mxst (省略可) 記録する(Ylagによりさかのぼることができる)ステップ数. (省略時 = 500) (Mxst <= 0 であれば省略時の既定値とみなす) [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 は不連続点) Function Phi(I As Long, T As Double) As Double If I = 1 Then Phi = 0.1 End Function Sub F4(N As Long, T As Double, Y() As Double, Yp() As Double, RCont As Double, ICont As Long) Dim Y2L1 As Double, Y2L10 As Double Y2L1 = Ylag(1, T - 1, AddressOf Phi, RCont, ICont) Y2L10 = Ylag(1, T - 10, AddressOf Phi, RCont, ICont) Yp(0) = -Y(0) * Y2L1 + Y2L10 Yp(1) = Y(0) * Y2L1 - Y(1) Yp(2) = Y(1) - Y2L10 End Sub Sub Ex_Retard() Const N = 3, NTout = 12, Ngrid = 0 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, Info As Long, 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 Call Retard(N, AddressOf F4, T, Y(), Tout(I), RTol(), ATol(), Ngrid, Tgrid(), Info) If Info <> 1 Then Debug.Print "Error in Retard: Info =", Info Exit For End If Debug.Print T, Y(0), Y(1), Y(2) Next End Sub 実行結果 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 は不連続点) Function Phi(I As Long, T As Double) As Double If I = 1 Then Phi = 0.1 End Function Sub F4(N As Long, T As Double, Y() As Double, Yp() As Double, RCont As Double, ICont As Long) Dim Y2L1 As Double, Y2L10 As Double Y2L1 = Ylag(1, T - 1, AddressOf Phi, RCont, ICont) Y2L10 = Ylag(1, T - 10, AddressOf Phi, RCont, ICont) Yp(0) = -Y(0) * Y2L1 + Y2L10 Yp(1) = Y(0) * Y2L1 - Y(1) Yp(2) = Y(1) - Y2L10 End Sub Sub Ex_Retard_2() Const N = 3, Ngrid = 12 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tgrid(Ngrid - 1) As Double Dim RTol(0) As Double, ATol(0) As Double, Info As Long, 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 Tgrid(I) = I + 1 Next Tgrid(10) = 20 Tgrid(11) = 40 Tend = 40 Info = 0 Call Retard(N, AddressOf F4, T, Y(), Tend, RTol(), ATol(), Ngrid, Tgrid(), Info, AddressOf SoloutRt) If Info <> 1 Then Debug.Print "Error in Retard: Info =", Info End Sub Sub SoloutRt(Nr As Long, Told As Double, T As Double, Y() As Double, N As Long, RCont As Double, ICont As Long, Irtrn As Long) Dim Y0 As Double, Y1 As Double, Y2 As Double Static Tout As Double If Nr = 1 Then Tout = 5 While T >= Tout Y0 = Ylag(0, Tout, AddressOf Phi, RCont, ICont) Y1 = Ylag(1, Tout, AddressOf Phi, RCont, ICont) Y2 = Ylag(2, Tout, AddressOf Phi, RCont, ICont) Debug.Print Tout, Y0, Y1, Y2 Tout = Tout + 5 Wend End Sub 実行結果 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