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

Dopri5() Sub Dopri5 ( N As  Long, F As  LongPtr, T As  Double, Y() As  Double, Tout As  Double, RTol() As  Double, ATol() 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 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 = t0 において y = y0 の解を求める. ただし, t0およびy0は既知でそれぞれtおよびyの初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. Dopri5は, 5(4)次 ドルマン・プリンス法に基づく陽的ルンゲ・クッタ法のプログラムである. ステップ制御アルゴリズムと密出力機能を持っており, 出力点数が非常に多くなったとしても効率的である. 詳細については下記参考文献を参照せよ. 引数 [in] N 微分方程式の数. (N >= 1) [in] F 微分方程式の値を求めるユーザー定義サブルーチンで, 次のように定義すること. Sub F(N As Long, T As Double, Y() As Double, Yp() As Double) Yp(i) = TおよびY()における微分係数の計算値 (i = 0 〜 N-1) End Sub ただし, Nは方程式の数, Yp()は与えられたTおよびY()における微分係数の計算値, すなわち, Yp(i) = dyi/dt = fi(T, Y(0), ..., Y(N-1)) (i = 0 〜 N-1). Yp()以外の変数を変更してはならない. [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,out] Info [in] = 0: 初期化して計算を開始する(新たな問題を解く). = 1: Toutだけを変えて計算を続ける(前回呼び出しの計算を再開する). [out] = -1: パラメータ N の誤り. (N < 1) = -4: パラメータ Y() の誤り. = -6: パラメータ RTol() の誤り. (RTol(i) < 0, RTol(i) = 0 かつ ATol(i) = 0) = -7: パラメータ ATol() の誤り. (ATol(i) < 0) = -8: パラメータ Info の誤り. (Info <> 0 かつ Info <> 1) = 1: 正常終了. = 2: Soloutによる中断 (正常終了). = 11: 最大ステップ数を超えた. = 12: ステップサイズが小さくなり過ぎた. = 13: スティフな問題の可能性がある (中断された). [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) = Contd5(i, T2, RCont, ICont); End Sub これを省略した場合(Solout = NullPtrの場合), 中間結果出力を行わない. [out] Neval (省略可) 関数評価回数. [out] Nstep (省略可) ステップ数. [out] Naccept (省略可) 採用されたステップ数. [out] Nreject (省略可) 不採用だったステップ数. (最初のステップはカウントされない) [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) 次の常微分方程式の初期値問題を解く. 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_Dopri5() Const N = 2 Dim T As Double, Y(N - 1) As Double, Tend As Double, Tout 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: Tend = 10: Y(0) = 1: Y(1) = 2 Info = 0 Do Tout = T + 1 Call Dopri5(N, AddressOf F1, T, Y(), Tout, RTol(), ATol(), Info) If Info <> 1 Then Debug.Print "Error in Dopri5: Info =", Info Exit Do End If Debug.Print T, Y(0), Y(1) Loop While Tout < Tend End Sub 実行結果 1 0.367879441193688 0.908181746995853 2 0.13533528322734 -0.280811553289576 3 4.97870683500648E-02 -0.940205428195432 4 1.83156388723754E-02 -0.635327981941068 5 6.73794700675574E-03 0.290400132447405 6 2.4787521939019E-03 0.962649038792233 7 9.1188198233568E-04 0.75481413627468 8 3.35462627424841E-04 -0.145164571180089 9 1.2340978747767E-04 -0.911006852047013 10 4.53999120696984E-05 -0.839026129110666 使用例 (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_Dopri5_2() Const N = 2 Dim T As Double, Y(N - 1) As Double, Tend 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: Tend = 10: Y(0) = 1: Y(1) = 2 Info = 0 Call Dopri5(N, AddressOf F1, T, Y(), Tend, RTol(), ATol(), Info, AddressOf SoloutD5) If Info <> 1 Then Debug.Print "Error in Dopri5: Info =", Info End Sub Sub SoloutD5(Nr As Long, Told As Double, T As Double, Y() As Double, N As Long, RCont As Double, ICont As Long, rtrn As Long) Dim Y0 As Double, Y1 As Double Static Tout As Double If Nr = 1 Then Tout = 1 While T >= Tout Y0 = Contd5(0, Tout, RCont, ICont) Y1 = Contd5(1, Tout, RCont, ICont) Debug.Print Tout, Y0, Y1 Tout = Tout + 1 Wend End Sub 実行結果 1 0.367879441187881 0.908181747006775 2 0.135335283237205 -0.280811553309717 3 4.97870683667704E-02 -0.940205428229821 4 1.83156388900869E-02 -0.635327981977462 5 6.73794699879586E-03 0.290400132463192 6 2.47875219132359E-03 0.962649038797553 7 9.11881968405946E-04 0.754814136303386 8 3.3546265009956E-04 -0.145164571227014 9 1.23409786683808E-04 -0.911006852045392 10 4.53999126645588E-05 -0.83902612911184