Contd5 Contd5_r Contd8 Contd8_r Contx1 Contx1_r Contx2 Contx2_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

Odex_r() Sub Odex_r ( N As  Long, T As  Double, Y() As  Double, Tout As  Double, RTol() As  Double, ATol() 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 Hinit As  Double = 0, Optional Hmax As  Double = 0, Optional MaxIter As  Long = 0, Optional Km As  Long = 0, Optional Nsequ As  Long = 0, Optional Mstab As  Long = 0, Optional Jstab As  Long = 0, Optional Iderr As  Long = 0, Optional Mudif As  Long = 0, Optional Fac1 As  Double = 0, Optional Fac2 As  Double = 0, Optional Fac3 As  Double = 0, Optional Fac4 As  Double = 0, Optional Safe1 As  Double = 0, Optional Safe2 As  Double = 0, Optional Safe3 As  Double = 0, Optional Cnt As  Long = 0  ) 常微分方程式の初期値問題 (補外法 (GBSアルゴリズム)) (リバースコミュニケーション版) 注 - 本プログラムは次バージョンで廃止予定です. 目的 本ルーチンは1階の常微分方程式の初期値問題 dy/dt = f(t, y), ただし t = t0 において y = y0 の解を求める. ただし, t0およびy0は既知でそれぞれtおよびyの初期値である. 上の方程式が連立微分方程式であれば, yはベクトルで表される. Odexは陽的中点則にもとづく補外法(GBSアルゴリズム)のプログラムである. ステップ制御, 次数選択および密出力機能を持っている. 詳細については下記参考文献を参照せよ. Odex_rはOdexのリバースコミュニケーション版である. 引数 [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,out] RCont() 配列 RCont(LRCont - 1) (LRCont >= (2*Km + 5)*N) 密出力のための制御情報. (Iout = 0 のときには参照されない) [in,out] ICont() 配列 ICont(LICont - 1) (LICont >= N) 密出力のための制御情報. (Iout = 0 のときには参照されない) [in,out] Info [in] = 0: 初期化して計算を開始する(新たな問題を解く). = 1: Toutだけを変えて計算を続ける(前回呼び出しの計算を再開する). [out] = -1: パラメータ N の誤り. (N < 1) = -3: パラメータ Y() の誤り. = -5: パラメータ RTol() の誤り. (RTol(i) < 0, RTol(i) = 0 かつ ATol(i) = 0) = -6: パラメータ ATol() の誤り. (ATol(i) < 0) = -7: パラメータ RCont() の誤り. = -8: パラメータ ICont() の誤り. = -9: パラメータ Info の誤り. (Info <> 0 かつ Info <> 1) = -11: パラメータ YY() の誤り. = -12: パラメータ YYp() の誤り. = -14: パラメータ IRev の誤り. (IRev != 0, 1, 5) = 1: 正常終了. = 2: Irtrn による中断 (正常終了). = 11: 最大ステップ数を超えた. [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 で終了する. また, 解を変更した場合には Irtrn = 3 と設定する必要がある. [out] IRev = 5: 初回, 中間, 最後にIRev = 5 で戻るとき, それぞれ 0, 1, 2 を返す. [in,out] IRev リバースコミュニケーションの制御変数 [in] 最初の呼び出し時に 0 に設定しておくこと. 2回目以降の呼び出し時には値を変更してはならない. [out] 0 以外の場合, 下記処理を行いIRevを変更せずに再び本ルーチンを呼び出すこと. = 0: 処理終了. 正常終了かどうかはInfoをチェックすること. = 1: TTおよびYY()における微分係数の計算値をYYp()に設定する. YYp()以外の変数を変更してはならない. = 5: 途中結果を出力する. 通常はどの変数も変更してはならない(Ioutを参照せよ). [in] Iout (省略可) 計算の途中結果を出力するかどうか指定する. (省略時 = 0) = 0: 出力しない. (IRev = 5では戻らない) = 1: 途中結果出力のために成功したステップごとに IRev = 5 で戻る. リバースコミュニケーション版でないOdexのSoloutの呼び出しに相当する. 対応する情報は次のとおりである. Nr = Naccept + 1, Told = 前回のT, T = 現在のT, Y() = 現在のY(). RCont()およびICont()はSoloutと同様に密出力のために使用する. Y(i) = Contx1_r(i, T2, RCont(), ICont()) (Ioutが上記以外の値であれば省略時の既定値とみなす) [out] Neval (省略可) 関数評価回数. [out] Nstep (省略可) ステップ数. [out] Naccept (省略可) 採用されたステップ数. [out] Nreject (省略可) 不採用だったステップ数. (最初のステップはカウントされない) [in] Hinit (省略可) ステップ幅の初期値. (省略時 = 0.0001) H = 1/||f'|| であり, 通常 0.1〜0.001 が適切とされる. (|Hinit| < 0.0001 であれば Hinit = 0.0001 とみなす) [in] Hmax (省略可) ステップ幅の最大値. (省略時 = Tout - T) (Hmax = 0 であれば省略時の既定値とみなす) [in] MaxIter (省略可) 許されるステップ数の最大値. (省略時 = 100000) (MaxIter <= 0 であれば省略時の既定値とみなす) [in] Km (省略可) 補外表の最大列数. (Km >= 3) (省略時 = 9) (Km < 3 であれば省略時の既定値とみなす) [in] Nsequ (省略可) ステップ数列の選択. (省略時 = 1 (Solout不使用時), 4 (Solout使用時)) = 1: 2, 4, 6, 8, 10, 12, 14, 16, ... = 2: 2, 4, 8, 12, 16, 20, 24, 28, ... = 3: 2, 4, 6, 8, 12, 16, 24, 32, ... = 4: 2, 6, 10, 14, 18, 22, 26, 30, ... = 5: 4, 8, 12, 16, 20, 24, 28, 32, ... Solout使用時には1〜3は指定できない. (上記以外の値, あるいは, Solout使用時に1〜3であれば省略時の既定値とみなす) [in] Mstab (省略可) 補外表の1行につき最大でMstab回の安定性チェックを行う. (省略時 = 1) (Mstab <= 0 であれば省略時の既定値とみなす) [in] Jstab (省略可) 補外表の安定性のチェックを行1〜行Jstabについてのみ行う. (省略時 = 2) (Jstab <= 0 であれば省略時の既定値とみなす) [in] Iderr (省略可) 密出力の式における誤差推定. (省略時 = 0) = 0: 行う = 1: 行わない Soloutを使用しない場合, Iderrの指定にかかわらず誤差推定は行わない. (上記以外の値であれば Iderr = 1 とみなす) [in] Mudif (省略可) 補外式の次数を決めるパラメータ. (1 <= Mudif <= 6) (省略時 = 4) (Mudif < 1 または Mudif > 6 であれば省略時の既定値とみなす) [in] Fac1 (省略可) [in] Fac2 (省略可) ステップ幅選択パラメータ. (省略時: Fac1 = 0.02, Fac2 = 4) j番目の要素が Facmin/Fac2 <= j番目のHnew/Hold <= 1/Facmin となるように新しいステップ幅が選ばれる. ただし, Facmin = Fac1^(1/(2*j - 1)) である. (Fac1 = 0 あるいは Fac2 = 0 であればそれぞれ省略時の既定値とみなす) [in] Fac3 (省略可) [in] Fac4 (省略可) 次数選択パラメータ. (省略時: Fac3 = 0.8, Fac4 = 0.9)   W(k-1) <= W(k)*Fac3 : 次数を減らす   W(k) <= W(k-1)*Fac4 : 次数を増やす W(k) = A(k)/H(k)は単位ステップあたりの仕事量である. (A(k)は関数値計算回数, H(k)はステップ幅を表す) (Fac3 = 0 あるいは Fac4 = 0 であればそれぞれ省略時の既定値とみなす) [in] Safe1 (省略可) [in] Safe2 (省略可) ステップ幅制御の安全係数. (省略時: Safe1 = 0.65, Safe2 = 0.94) (Safe1 = 0 あるいは Safe2 = 0 であればそれぞれ省略時の既定値とみなす) [in] Safe3 (省略可) ステップ幅を減らす際の係数. (省略時 = 0.5) (Safe3 <= 0 であれば省略時の既定値とみなす. Safe3 > 1 であれば 1 とみなす) [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_Odex_r() 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 Dim RCont() As Double, ICont() As Long, TT As Double Dim YY(N - 1) As Double, YYp(N - 1) As Double, Irtrn As Long, IRev 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 IRev = 0 Do Call Odex_r(N, T, Y(), Tout, RTol(), ATol(), RCont(), ICont(), Info, TT, YY(), YYp(), Irtrn, IRev) If IRev = 1 Then Call F1(N, TT, YY(), YYp()) Loop While IRev <> 0 If Info <> 1 Then Debug.Print "Error in Odex_r: Info =", Info Exit Do End If Debug.Print T, Y(0), Y(1) Loop While Tout < Tend End Sub 実行結果 1 0.367879441180513 0.90818174702144 2 0.135335283238428 -0.280811553314161 3 4.97870683602588E-02 -0.940205428217376 4 1.83156388764711E-02 -0.635327981950352 5 6.73794699120742E-03 0.290400132478066 6 2.47875218511396E-03 0.962649038810137 7 9.11881969599684E-04 0.754814136300769 8 3.35462632178912E-04 -0.145164571189262 9 1.23409800123739E-04 -0.911006852072663 10 4.53999253959451E-05 -0.839026129137953 使用例 (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_Odex_r_2() Const N = 2, KM = 9 Dim T As Double, Y(N - 1) As Double, Tend As Double Dim RTol(0) As Double, ATol(0) As Double, Info As Long Dim RCont((2 * KM + 5) * N - 1) As Double, ICont(N - 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, Tout As Double, Y0 As Double, Y1 As Double RTol(0) = 0.0000000001 '1.0e-10 ATol(0) = RTol(0) Iout = 1 T = 0: Tend = 10: Y(0) = 1: Y(1) = 2 Tout = 1 Info = 0 IRev = 0 Do Call Odex_r(N, T, Y(), Tend, RTol(), ATol(), RCont(), ICont(), Info, TT, YY(), YYp(), Irtrn, IRev, Iout) If IRev = 1 Then Call F1(N, TT, YY(), YYp()) ElseIf IRev = 5 Then While T >= Tout Y0 = Contx1_r(0, Tout, RCont(), ICont()) Y1 = Contx1_r(1, Tout, RCont(), ICont()) Debug.Print Tout, Y0, Y1 Tout = Tout + 1 Wend End If Loop While IRev <> 0 If Info <> 1 Then Debug.Print "Error in Odex_r: Info =", Info End Sub 実行結果 1 0.367879441157644 0.908181747067181 2 0.135335283236881 -0.280811553311066 3 4.97870683669239E-02 -0.940205428230704 4 1.83156389169758E-02 -0.635327982031361 5 6.73794699664751E-03 0.290400132467187 6 2.47875227629557E-03 0.962649038627776 7 9.11881939436941E-04 0.754814136361094 8 3.35462629214848E-04 -0.145164571183337 9 1.23409825870176E-04 -0.911006852124154 10 4.53998495624352E-05 -0.839026128986288